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