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