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