Annotation of OpenXM/src/k097/debug/toric0.k, Revision 1.1.1.1
1.1 maekawa 1:
2: /* SSWork/yacc/ip-k/ex3.ccc, 1996, 8/12. This is original. ---> */
3: /* debug/toric0.k */
4: /* graver basis を求める関数.
5: toric の generator を求める関数.
6: A-hypergeometric の Fourier 変換の indicial を求める関数.
7: */
8:
9:
10: def toric0_toMonom(aa,i,offset, ring)
11: {
12: local j,ans,m;
13: m = Length(aa);
14: ans = PolyR("1",ring);
15: for (j=0; j<m; j++) {
16: ans = ans*(z[offset+j]^(aa[j,i]));
17: }
18: return(ans);
19: }
20:
21: def testgraver() {
22: local a,ans;
23: a = [[1,1,1,1,1,1,1],
24: [2,3,2,2,2,0,1],
25: [2,2,3,2,2,0,1],
26: [2,2,2,3,2,0,1],
27: [1,1,1,1,2,0,1]];
28: ans = Graver(a);
29: return(ans); /* degree 33. About 1 hour? cf. this/Mirror/wp22211.ml */
30: }
31:
32: def testgraver2() {
33: local a,ans;
34: a = [[1,1,1,1,1,1],
35: [0,0,0,1,1,1],
36: [0,1,0,0,1,0],
37: [0,0,1,0,0,1]];
38: ans = Graver(a);
39: return(ans);
40: }
41:
42:
43:
44: HelpAdd(["Graver",
45: ["Graver(a) (matrix a ) returns a graver basis of the affine toric ideal",
46: "defined by the matrix << a >>.",
47: "It defines a ring of variables z[0], z[1], ..., z[n-1], ..... and ",
48: "the output belongs to this ring.",
49: "Example: Graver([[1,1,1,1],[0,1,2,3]]):",
50: "[ -z[1]^2+z[2]*z[0] , -z[2]^2+z[3]*z[1] , -z[2]*z[1]+z[3]*z[0] ,",
51: " -z[1]^3+z[3]*z[0]^2 , -z[2]^3+z[3]^2*z[0] ] "
52: ]]);
53:
54: def Graver(a) {
55: local aa,i,j,rz,n,d,ideal,ans,univ,rule,nn,weight,elim;
56: d = Length(a); n = Length(a[0]);
57: aa = NewMatrix(d+n,2*n);
58: for (i=0; i<d; i++) {
59: for (j=0; j<n; j++) {
60: aa[i,j] = a[i,j];
61: }
62: }
63: for (i=0; i<n; i++) {
64: aa[d+i,i] = 1;
65: aa[d+i,n+i] = 1;
66: }
67: Println(aa);
68:
69: weight = [ ]; elim = [ ];
70: for (i= 2*n; i< 2*n + n+d; i++) {
71: weight = Join(weight,[Indexed("z",i), 1]);
72: elim = Append(elim, Indexed("z",i));
73: }
74: Println(weight);
75: Println(elim);
76:
77: rz = RingPonIndexedVariables("z",2*n+d+n, [weight]);
78: /* z[0], ..., z[2n-1], ... */
79: ideal = [ ];
80: for (i=0; i< 2*n; i++) {
81: ideal = Append(ideal, z[i] - toric0_toMonom(aa,i,2*n,rz));
82: }
83: Println(" --------- input ideal -------------");
84: Print(" z["); Print( 2*n ); Print( "] --- z["); Print( 2*n+n+d-1);
85: Println("] should be eliminated.");
86: Println(ideal);
87:
88: ans = Groebner(ideal);
89: Println(" -------------- gb is ----------------- "); Println(ans);
90: ans = Eliminatev(ans,elim);
91: Println(" ------------ eliminated -------------- "); Println(ans);
92:
93: rule = [[h, PolyR("1",rz) ] ];
94: nn = Length(ans); univ = [ ];
95: for (i=0; i<nn; i++) {
96: univ = Append(univ,Replace(ans[i],rule));
97: }
98: ans = ReducedBase(univ);
99: Println(" ----------- removed redundant elements ----------- ");
100: Println(ans);
101:
102: rule = [ ];
103: for (i= n; i<2*n; i++) {
104: rule = Append(rule, [ z[i], PolyR("1",rz)]);
105: }
106: Println(" ------------ a rule for replacement ------------ ");
107: Println(rule);
108:
109: univ = [ ];
110: nn = Length( ans );
111: for (i=0; i<nn; i++) {
112: univ = Append(univ,Replace(ans[i],rule));
113: }
114: Println(" ------------ graver basis is ------------ ");
115: Println(univ);
116: Println(" ");
117: return(univ);
118: }
119:
120:
121: def toric(aa) {
122: local i,j,rz,n,d,ideal,ans,univ,rule,nn,weight,elim;
123: d = Length(aa); n = Length(aa[0]);
124: Println(aa);
125:
126: weight = [ ]; elim = [ ];
127: for (i= n; i< n+d; i++) {
128: weight = Join(weight,[Indexed("z",i), 1]);
129: elim = Append(elim, Indexed("z",i));
130: }
131: Println(weight);
132: Println(elim);
133:
134: rz = RingPonIndexedVariables("z",n+d, [weight]);
135: /* z[0], ..., z[n-1], ... , z[n+d]*/
136: ideal = [ ];
137: for (i=0; i< n; i++) {
138: ideal = Append(ideal, z[i] - toric0_toMonom(aa,i,n,rz));
139: }
140: Println(" --------- input ideal -------------");
141: Print(" z["); Print( n ); Print( "] --- z["); Print( n+d-1);
142: Println("] should be eliminated.");
143: Println(ideal);
144:
145: ans = Groebner(ideal);
146: Println(" -------------- gb is ----------------- "); Println(ans);
147: ans = Eliminatev(ans,elim);
148: Println(" ------------ eliminated -------------- "); Println(ans);
149:
150: rule = [[h, PolyR("1",rz) ] ];
151: nn = Length(ans); univ = [ ];
152: for (i=0; i<nn; i++) {
153: univ = Append(univ,Replace(ans[i],rule));
154: }
155: ans = ReducedBase(univ);
156: Println(" ----------- removed redundant elements ----------- ");
157: Println(" ---------- generators of the toric ideal are ----- ");
158: Println(ans);
159: Println(" ");
160: return(ans);
161: }
162:
163: def zindicial0(input,n,m) {
164: local rz,weight, ww,i,rule,zinverse,m,d,ans,elim;
165: ww = [ ]; elim = [ ];
166: /* z[n+1], ..., z[n+m] がパラメータ変数 */
167: /* Dz[0] --- Dz[n-2], z[0] --- z[n-2] を消去する. */
168: for (i=0; i<n-1; i++) {
169: ww = Join(ww,[Indexed("Dz",i), 1]);
170: ww = Join(ww,[Indexed("z",i), 1]);
171: if (i != n-1) {
172: elim = Append(elim,Indexed("Dz",i));
173: elim = Append(elim,Indexed("z",i));
174: }
175: }
176: weight = [[Indexed("z",n),1] , ww];
177: Print("-------- weight ---------: "); Println(weight);
178: rz = RingDonIndexedVariables("z",n+1+m, weight);
179: input = Mapto(input,rz);
180: Println("------------ input ------------"); Println(input);
181:
182: /* F-homogenization. z[0], ..., z[n-1],
183: z[n] is the homogenization variable*/
184: /* z[n]^(-1) とは書けないのはつらい. 1 を戻すので bug ともいえる. */
185: zinverse = PolyR(AddString([Indexed("z",n),"^(-1)"]),rz);
186: rule = [[Dz[n-1], Dz[n-1]*z[n]], [z[n-1],z[n-1]*zinverse]];
187: input = Replace(input,rule);
188: m = Length(input);
189: for (i=0; i<m; i++) {
190: d = -Degree(Replace(input[i],[[z[n],zinverse]]),z[n]);
191: if (d < 0) {
192: input[i] = z[n]^(-d)*input[i];
193: }
194: }
195: Print("------ input : "); Println(input);
196: ans = Groebner(input);
197: m = Length(ans);
198: for (i=0; i<m; i++) {
199: /* FW principal part をとる. */
200: ans[i] = Init_w(ans[i],[z[n]],[1]);
201: }
202: Print("--------FW principal parts : ");Println(ans);
203: ans = Eliminatev(ans,elim);
204: m = Length(ans);
205: /* h, z[n] を 1 にする. */
206: for (i=0; i<m; i++) {
207: ans[i] = Replace(ans[i],[[h,PolyR("1",rz)],[z[n],PolyR("1",rz)]]);
208: }
209: Println(" "); Println(" ");
210: return(ans);
211: }
212:
213: def zrho(f,n) {
214: local ans,i,top,w,rz;
215: ans = 0;
216: rz = Ringp(f);
217: while(true) {
218: if ( f == Poly("0")) sm1(" exit ");
219: top = Init(f);
220: f = f-top;
221: w = Exponent(top,[Dz[n-1]]);
222: top = Replace(top,[[Dz[n-1],PolyR("1",rz)]])*zipoch(z[n],w[0]);
223: ans = ans + top;
224: }
225: return(ans);
226: }
227:
228: def zipoch(f,w) {
229: local ans,i;
230: ans = 1;
231: for (i=0; i<w; i++) {
232: ans = ans*(f-i);
233: }
234: return(ans);
235: }
236:
237:
238: fff = ["z[0]*Dz[0]+z[1]*Dz[1]+z[2]*Dz[2]+z[3]*Dz[3]-z[5]",
239: " z[1]*Dz[1] +z[3]*Dz[3]-z[6]",
240: " z[2]*Dz[2]+z[3]*Dz[3]+z[7]",
241: "z[0]*z[3]-z[1]*z[2]"];
242: /* z[4] はつかったらだめよ. */
243: /* Println(" type in zindicial0(fff,4,3); "); */
244:
245: def foo2() {
246: local ans,n,i,m,r,tmp;
247: ans = zindicial0(fff,4,3); n = 4;
248: Println(ans);
249: m = Length(ans);
250: r = [ ];
251: for (i=0; i<m; i++) {
252: tmp = ans[i];
253: tmp = Replace(tmp,[[z[n-1],Poly("1")]]);
254: tmp = zrho(tmp,n);
255: Print(i); Print(" : ");Println(tmp);
256: r = Append(r,tmp);
257: }
258: return(r);
259: }
260:
261: /* Println("Type in foo2()."); */
262:
263: def zindicial(fff,n,mm) {
264: local ans,n,i,m,r,tmp;
265: ans = zindicial0(fff,n,mm);
266: Println(ans);
267: m = Length(ans);
268: r = [ ];
269: Println(AddString(["------ The indicial polynomial along z[",
270: ToString(n-1),
271: "] = 0 is the minimal degree polynomial of the following",
272: "polynomials."]));
273: Println(AddString(["z[",ToString(n),"] is equal to s."]));
274: for (i=0; i<m; i++) {
275: tmp = ans[i];
276: tmp = Replace(tmp,[[z[n-1],Poly("1")]]);
277: tmp = zrho(tmp,n);
278: Print(i); Print(" : ");Println(tmp);
279: r = Append(r,tmp);
280: }
281: Println(" ");
282: return(r);
283: }
284:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to toric0.k CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / k097 / debug