Annotation of OpenXM/src/k097/debug/ahg2.k, Revision 1.2
1.1 maekawa 1:
2: /* SSWork/yacc/debug/ahg.k
3: /* cf. debug/toric0.k */
4: /* toric の generator を求める関数.
5: A-hypergeometric の indicial ideal を求める関数.
6: */
1.2 ! takayama 7: load("indexed.k");
! 8:
1.1 maekawa 9: ShimomuraSpecial = true ; OnePath= true;
10: Vvv = false;
11: SetRingVariables_Verbose = false;
12: def void QuietKan() {
13: sm1(" [(KanGBmessage) 0] system_variable ");
14: }
15:
16: def testhg1() {
17: a = [[1,1,1,1,1,1],
18: [0,0,0,1,1,1],
19: [0,1,0,0,1,0],
20: [0,0,1,0,0,1]];
21: return(idhg(a));
22: }
23:
24: def testhg2() {
25: a = [[1,1,1,1,1],
26: [0,2,3,4,3],
27: [0,1,1,0,2]];
28: return(idhg(a));
29: }
30:
31: def idhg(a) {
32: local a,ans,rd,i,ans2,ans3,n,ff,d,zlist;
33: ans = toric(a);
34: if (ShimomuraSpecial) {
35: /* 先に, toric の initial part をとってしまう. */
36: /* 本当は, FW の initial part をとるべきなのかも? */
37: if (Vvv) {Println("-------- S-special ---------");}
38: ans = Map(ans,"Init");
39: }
40: ans = Map(ans,"ToString");
41: if (Vvv) {Println(ans);}
42:
43: rd = RingDonIndexedVariables("z",Length(a[0])+1+Length(a));
44: /* 4 秒程度かかる. */
45: ans = Map(ans,"Poly");
46: n = Length(a[0]); d = Length(a);
47: ans2 = NewArray(Length(ans)); /* ans2 には, toric */
48: PSfor (i=0; i< Length(ans); i++) {
49: ans2[i] = ztoDz(ans[i],n);
50: }
51: if (Vvv) {Println(ans2);}
52: ans3 = atolin(a); /* ans3 には, 一次式 */
53: if (Vvv) {Println(ans3);}
54: ff = Map(Join(ans2,ans3),"ToString");
55: ans = zindicial(ff,n,d);
56: zlist = [ ];
57: PSfor(i= n; i<n+d+1; i++) {
58: zlist = Append(zlist,Indexed("z",i));
59: }
60: return([ans,zlist]);
61: }
62:
63: def toric0_toMonom(aa,i,offset, ring)
64: {
65: local j,ans,m;
66: m = Length(aa);
67: ans = PolyR("1",ring);
68: for (j=0; j<m; j++) {
69: ans = ans*(z[offset+j]^(aa[j,i]));
70: }
71: return(ans);
72: }
73:
74:
75: def toric(aa) {
76: local i,j,rz,n,d,ideal,ans,univ,rule,nn,weight,elim;
77: d = Length(aa); n = Length(aa[0]);
78: if (Vvv) {Println(aa);}
79:
80: weight = [ ]; elim = [ ];
81: PSfor (i= n; i< n+d; i++) {
82: weight = Join(weight,[Indexed("z",i), 1]);
83: elim = Append(elim, Indexed("z",i));
84: }
85: weight = Append([weight],[Indexed("z",n-1),1]);
86: if (Vvv) {Println(weight); Println(elim);}
87:
88: rz = RingPonIndexedVariables("z",n+d, weight); /* 4 秒程度かかる. */
89: /* z[0], ..., z[n-1], ... , z[n+d]*/
90: ideal = [ ];
91: PSfor (i=0; i< n; i++) {
92: ideal = Append(ideal, z[i] - toric0_toMonom(aa,i,n,rz));
93: }
94: if (Vvv) {
95: Println(" --------- input ideal -------------");
96: Print(" z["); Print( n ); Print( "] --- z["); Print( n+d-1);
97: Println("] should be eliminated.");
98: Println(ideal);
99: }
100:
101: ans = Groebner(ideal);
102: if (Vvv) {Println(" -------------- gb is ----------------- "); Println(ans);}
103: ans = Eliminatev(ans,elim);
104: if (Vvv) {Println(" ------------ eliminated -------------- "); Println(ans);}
105:
106: rule = [[h, PolyR("1",rz) ] ];
107: nn = Length(ans); univ = [ ];
108: PSfor (i=0; i<nn; i++) {
109: univ = Append(univ,Replace(ans[i],rule));
110: }
111: ans = ReducedBase(univ);
112: if (Vvv) {
113: Println(" ----------- removed redundant elements ----------- ");
114: Println(" ---------- generators of the toric ideal are ----- ");
115: Println(ans);
116: Println(" ");
117: }
118: return(ans);
119: }
120:
121: def zindicial0(input,n,m) {
122: local rz,weight, ww,i,rule,zinverse,m,d,ans,elim,tmp;
123: if (!OnePath) {
124: ww = [ ]; elim = [ ];
125: weight = [[Indexed("z",n),1]];
126: if (Vvv) {Print("-------- weight ---------: "); Println(weight);}
127: rz = RingDonIndexedVariables("z",n+1+m, weight);
128: z = NewArray(n+1+m); Dz = NewArray(n+1+m);
129: PSfor(i=0; i<n+1+m; i++) {
130: z[i] = PolyR(Indexed("z",i),rz);
131: Dz[i] = PolyR(Indexed("Dz",i),rz);
132: }
133: input = Mapto(input,rz);
134: if (Vvv) {Println("------------ input ------------"); Println(input);}
135:
136: /* F-homogenization. z[0], ..., z[n-1],
137: z[n] is the homogenization variable*/
138: /* z[n]^(-1) とは書けないのはつらい. 1 を戻すので bug ともいえる. */
139: zinverse = PolyR(AddString([Indexed("z",n),"^(-1)"]),rz);
140: rule = [[Dz[n-1], Dz[n-1]*z[n]], [z[n-1],z[n-1]*zinverse]];
141: input = Replace(input,rule);
142: m = Length(input);
143: PSfor (i=0; i<m; i++) {
144: d = -Degree(Replace(input[i],[[z[n],zinverse]]),z[n]);
145: if (d < 0) {
146: input[i] = z[n]^(-d)*input[i];
147: }
148: }
149: if (Vvv) {Print("------ input : "); Println(input);}
150: ans = GroebnerTime(input);
151: m = Length(ans);
152: PSfor (i=0; i<m; i++) {
153: /* FW principal part をとる. */
154: /* ans[i] = Init_w(ans[i],[z[n]],[1]); この関数は遅い. */
155: tmp = Coefficients(ans[i],z[n]);
156: ans[i] = tmp[1,0];
157: }
158: if (Vvv) {Print("--------FW principal parts : ");Println(ans);}
159: /* もう一回, GB の計算. */
160: input = Map(ans,"ToString");
161: }
162:
163: /* もう一回, GB の計算. */
164: ww = [ ]; elim = [ ];
165: /* z[n+1], ..., z[n+m] がパラメータ変数 */
166: /* Dz[0] --- Dz[n-2], z[0] --- z[n-2] を消去する. */
167: PSfor (i=0; i<n-1; i++) {
168: ww = Join(ww,[Indexed("Dz",i), 1]);
169: /* ww = Join(ww,[Indexed("z",i), 1]); */
170: if (i != n-1) {
171: elim = Append(elim,Indexed("Dz",i));
172: /* elim = Append(elim,Indexed("z",i)); */
173: }
174: }
175: weight = [[Indexed("z",n),1] , ww];
176: if (Vvv) {Print("-------- weight ---------: "); Println(weight);}
177: rz = RingDonIndexedVariables("z",n+1+m, weight);
178: z = NewArray(n+1+m); Dz = NewArray(n+1+m);
179: PSfor(i=0; i<n+1+m; i++) { /* これをもう一度やらないと, mklib で sm1 に
180: したときのみ z が undefined
181: になる. どうしてか? */
182: z[i] = PolyR(Indexed("z",i),rz);
183: Dz[i] = PolyR(Indexed("Dz",i),rz);
184: }
185: input = Mapto(input,rz);
186: if (OnePath) {
187: /* F-homogenization. z[0], ..., z[n-1],
188: z[n] is the homogenization variable*/
189: /* z[n]^(-1) とは書けないのはつらい. 1 を戻すので bug ともいえる. */
190: zinverse = PolyR(AddString([Indexed("z",n),"^(-1)"]),rz);
191: rule = [[Dz[n-1], Dz[n-1]*z[n]], [z[n-1],z[n-1]*zinverse]];
192: input = Replace(input,rule);
193: m = Length(input);
194: PSfor (i=0; i<m; i++) {
195: d = -Degree(Replace(input[i],[[z[n],zinverse]]),z[n]);
196: if (d < 0) {
197: input[i] = z[n]^(-d)*input[i];
198: }
199: }
200: }
201: if (Vvv) {Print("------ input : "); Println(input);}
202: ans = GroebnerTime(input);
203: m = Length(ans);
204: PSfor (i=0; i<m; i++) {
205: /* FW principal part をとる. */
206: /* ans[i] = Init_w(ans[i],[PolyR(Indexed("z",n),rz)],[1]); */
207: tmp = Coefficients(ans[i],z[n]);
208: ans[i] = tmp[1,0];
209: }
210: if (Vvv) {Print("--------FW principal parts : ");Println(ans);}
211:
212:
213:
214: ans = Eliminatev(ans,elim);
215: m = Length(ans);
216: /* h, z[n] を 1 にする. */
217: for (i=0; i<m; i++) {
218: ans[i] = Replace(ans[i],[[h,PolyR("1",rz)],[PolyR(Indexed("z",n),rz),PolyR("1",rz)]]);
219: }
220: if (Vvv) {Println(" "); Println(" ");}
221: return(ans);
222: }
223:
224: def zrho(f,n) {
225: local ans,i,top,w,rz;
226: ans = 0;
227: rz = Ringp(f);
228: while(true) {
229: if ( f == Poly("0")) sm1(" exit ");
230: top = Init(f);
231: f = f-top;
232: w = Exponent(top,[PolyR(Indexed("Dz",n-1),rz)]);
233: top = Replace(top,[[PolyR(Indexed("Dz",n-1),rz),PolyR("1",rz)]])*zipoch(z[n],w[0]);
234: ans = ans + top;
235: }
236: return(ans);
237: }
238:
239: def zipoch(f,w) {
240: local ans,i;
241: ans = 1;
242: PSfor (i=0; i<w; i++) {
243: ans = ans*(f-i);
244: }
245: return(ans);
246: }
247:
248:
249:
250:
251: def zindicial(fff,n,mm) {
252: local ans,n,i,m,r,tmp;
253: ans = zindicial0(fff,n,mm);
254: if (Vvv) {Println(ans);}
255: m = Length(ans);
256: r = [ ];
257: if (Vvv) {
258: Println(AddString(["------ The generic indicial polynomial along z[",
259: ToString(n-1),
260: "] = 0 is the minimal degree polynomial of the following",
261: "polynomials."]));
262: Println(AddString(["z[",ToString(n),"] is equal to s."]));
263: }
264: PSfor (i=0; i<m; i++) {
265: tmp = ans[i];
266: tmp = Replace(tmp,[[Poly(Indexed("z",n-1)),Poly("1")]]);
267: tmp = zrho(tmp,n);
268: if (Vvv) {Print(i); Print(" : ");Println(tmp);}
269: r = Append(r,tmp);
270: }
271: if (Vvv) {Println(" ");}
272: return(r);
273: }
274:
275: def ztoDz(f,n) {
276: local rule,i;
277: rule = NewArray(n);
278: PSfor(i=0; i<n; i++) {
279: rule[i] = [z[i],Dz[i]];
280: }
281: return(Replace(f,rule));
282: }
283:
284: /* 行列より A-HG の線形方程式をだす. */
285: def atolin(a) {
286: local d,n,eqs,ans,i,j;
287: d = Length(a);
288: n = Length(a[0]);
289: eqs = NewArray(d);
290: PSfor (i=0; i<d; i++) {
291: ans = 0;
292: PSfor (j=0; j<n; j++) {
293: ans = ans + a[i,j]*z[j]*Dz[j];
294: }
295: ans = ans - z[n+1+i];
296: eqs[i] = ans;
297: }
298: return(eqs);
299: }
300:
301:
302:
303:
304: sm1(" ; ");
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ahg2.k CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / k097 / debug