Annotation of OpenXM/src/k097/debug/ahg.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: 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;
! 189: rz = Ringp(f);
! 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