Annotation of OpenXM/src/k097/debug/ahg2.k, Revision 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