Annotation of OpenXM/src/k097/debug/toric0.k, Revision 1.1
1.1 ! maekawa 1:
! 2: /* SSWork/yacc/ip-k/ex3.ccc, 1996, 8/12. This is original. ---> */
! 3: /* debug/toric0.k */
! 4: /* graver basis を求める関数.
! 5: toric の generator を求める関数.
! 6: A-hypergeometric の Fourier 変換の indicial を求める関数.
! 7: */
! 8:
! 9:
! 10: def toric0_toMonom(aa,i,offset, ring)
! 11: {
! 12: local j,ans,m;
! 13: m = Length(aa);
! 14: ans = PolyR("1",ring);
! 15: for (j=0; j<m; j++) {
! 16: ans = ans*(z[offset+j]^(aa[j,i]));
! 17: }
! 18: return(ans);
! 19: }
! 20:
! 21: def testgraver() {
! 22: local a,ans;
! 23: a = [[1,1,1,1,1,1,1],
! 24: [2,3,2,2,2,0,1],
! 25: [2,2,3,2,2,0,1],
! 26: [2,2,2,3,2,0,1],
! 27: [1,1,1,1,2,0,1]];
! 28: ans = Graver(a);
! 29: return(ans); /* degree 33. About 1 hour? cf. this/Mirror/wp22211.ml */
! 30: }
! 31:
! 32: def testgraver2() {
! 33: local a,ans;
! 34: a = [[1,1,1,1,1,1],
! 35: [0,0,0,1,1,1],
! 36: [0,1,0,0,1,0],
! 37: [0,0,1,0,0,1]];
! 38: ans = Graver(a);
! 39: return(ans);
! 40: }
! 41:
! 42:
! 43:
! 44: HelpAdd(["Graver",
! 45: ["Graver(a) (matrix a ) returns a graver basis of the affine toric ideal",
! 46: "defined by the matrix << a >>.",
! 47: "It defines a ring of variables z[0], z[1], ..., z[n-1], ..... and ",
! 48: "the output belongs to this ring.",
! 49: "Example: Graver([[1,1,1,1],[0,1,2,3]]):",
! 50: "[ -z[1]^2+z[2]*z[0] , -z[2]^2+z[3]*z[1] , -z[2]*z[1]+z[3]*z[0] ,",
! 51: " -z[1]^3+z[3]*z[0]^2 , -z[2]^3+z[3]^2*z[0] ] "
! 52: ]]);
! 53:
! 54: def Graver(a) {
! 55: local aa,i,j,rz,n,d,ideal,ans,univ,rule,nn,weight,elim;
! 56: d = Length(a); n = Length(a[0]);
! 57: aa = NewMatrix(d+n,2*n);
! 58: for (i=0; i<d; i++) {
! 59: for (j=0; j<n; j++) {
! 60: aa[i,j] = a[i,j];
! 61: }
! 62: }
! 63: for (i=0; i<n; i++) {
! 64: aa[d+i,i] = 1;
! 65: aa[d+i,n+i] = 1;
! 66: }
! 67: Println(aa);
! 68:
! 69: weight = [ ]; elim = [ ];
! 70: for (i= 2*n; i< 2*n + n+d; i++) {
! 71: weight = Join(weight,[Indexed("z",i), 1]);
! 72: elim = Append(elim, Indexed("z",i));
! 73: }
! 74: Println(weight);
! 75: Println(elim);
! 76:
! 77: rz = RingPonIndexedVariables("z",2*n+d+n, [weight]);
! 78: /* z[0], ..., z[2n-1], ... */
! 79: ideal = [ ];
! 80: for (i=0; i< 2*n; i++) {
! 81: ideal = Append(ideal, z[i] - toric0_toMonom(aa,i,2*n,rz));
! 82: }
! 83: Println(" --------- input ideal -------------");
! 84: Print(" z["); Print( 2*n ); Print( "] --- z["); Print( 2*n+n+d-1);
! 85: Println("] should be eliminated.");
! 86: Println(ideal);
! 87:
! 88: ans = Groebner(ideal);
! 89: Println(" -------------- gb is ----------------- "); Println(ans);
! 90: ans = Eliminatev(ans,elim);
! 91: Println(" ------------ eliminated -------------- "); Println(ans);
! 92:
! 93: rule = [[h, PolyR("1",rz) ] ];
! 94: nn = Length(ans); univ = [ ];
! 95: for (i=0; i<nn; i++) {
! 96: univ = Append(univ,Replace(ans[i],rule));
! 97: }
! 98: ans = ReducedBase(univ);
! 99: Println(" ----------- removed redundant elements ----------- ");
! 100: Println(ans);
! 101:
! 102: rule = [ ];
! 103: for (i= n; i<2*n; i++) {
! 104: rule = Append(rule, [ z[i], PolyR("1",rz)]);
! 105: }
! 106: Println(" ------------ a rule for replacement ------------ ");
! 107: Println(rule);
! 108:
! 109: univ = [ ];
! 110: nn = Length( ans );
! 111: for (i=0; i<nn; i++) {
! 112: univ = Append(univ,Replace(ans[i],rule));
! 113: }
! 114: Println(" ------------ graver basis is ------------ ");
! 115: Println(univ);
! 116: Println(" ");
! 117: return(univ);
! 118: }
! 119:
! 120:
! 121: def toric(aa) {
! 122: local i,j,rz,n,d,ideal,ans,univ,rule,nn,weight,elim;
! 123: d = Length(aa); n = Length(aa[0]);
! 124: Println(aa);
! 125:
! 126: weight = [ ]; elim = [ ];
! 127: for (i= n; i< n+d; i++) {
! 128: weight = Join(weight,[Indexed("z",i), 1]);
! 129: elim = Append(elim, Indexed("z",i));
! 130: }
! 131: Println(weight);
! 132: Println(elim);
! 133:
! 134: rz = RingPonIndexedVariables("z",n+d, [weight]);
! 135: /* z[0], ..., z[n-1], ... , z[n+d]*/
! 136: ideal = [ ];
! 137: for (i=0; i< n; i++) {
! 138: ideal = Append(ideal, z[i] - toric0_toMonom(aa,i,n,rz));
! 139: }
! 140: Println(" --------- input ideal -------------");
! 141: Print(" z["); Print( n ); Print( "] --- z["); Print( n+d-1);
! 142: Println("] should be eliminated.");
! 143: Println(ideal);
! 144:
! 145: ans = Groebner(ideal);
! 146: Println(" -------------- gb is ----------------- "); Println(ans);
! 147: ans = Eliminatev(ans,elim);
! 148: Println(" ------------ eliminated -------------- "); Println(ans);
! 149:
! 150: rule = [[h, PolyR("1",rz) ] ];
! 151: nn = Length(ans); univ = [ ];
! 152: for (i=0; i<nn; i++) {
! 153: univ = Append(univ,Replace(ans[i],rule));
! 154: }
! 155: ans = ReducedBase(univ);
! 156: Println(" ----------- removed redundant elements ----------- ");
! 157: Println(" ---------- generators of the toric ideal are ----- ");
! 158: Println(ans);
! 159: Println(" ");
! 160: return(ans);
! 161: }
! 162:
! 163: def zindicial0(input,n,m) {
! 164: local rz,weight, ww,i,rule,zinverse,m,d,ans,elim;
! 165: ww = [ ]; elim = [ ];
! 166: /* z[n+1], ..., z[n+m] がパラメータ変数 */
! 167: /* Dz[0] --- Dz[n-2], z[0] --- z[n-2] を消去する. */
! 168: for (i=0; i<n-1; i++) {
! 169: ww = Join(ww,[Indexed("Dz",i), 1]);
! 170: ww = Join(ww,[Indexed("z",i), 1]);
! 171: if (i != n-1) {
! 172: elim = Append(elim,Indexed("Dz",i));
! 173: elim = Append(elim,Indexed("z",i));
! 174: }
! 175: }
! 176: weight = [[Indexed("z",n),1] , ww];
! 177: Print("-------- weight ---------: "); Println(weight);
! 178: rz = RingDonIndexedVariables("z",n+1+m, weight);
! 179: input = Mapto(input,rz);
! 180: Println("------------ input ------------"); Println(input);
! 181:
! 182: /* F-homogenization. z[0], ..., z[n-1],
! 183: z[n] is the homogenization variable*/
! 184: /* z[n]^(-1) とは書けないのはつらい. 1 を戻すので bug ともいえる. */
! 185: zinverse = PolyR(AddString([Indexed("z",n),"^(-1)"]),rz);
! 186: rule = [[Dz[n-1], Dz[n-1]*z[n]], [z[n-1],z[n-1]*zinverse]];
! 187: input = Replace(input,rule);
! 188: m = Length(input);
! 189: for (i=0; i<m; i++) {
! 190: d = -Degree(Replace(input[i],[[z[n],zinverse]]),z[n]);
! 191: if (d < 0) {
! 192: input[i] = z[n]^(-d)*input[i];
! 193: }
! 194: }
! 195: Print("------ input : "); Println(input);
! 196: ans = Groebner(input);
! 197: m = Length(ans);
! 198: for (i=0; i<m; i++) {
! 199: /* FW principal part をとる. */
! 200: ans[i] = Init_w(ans[i],[z[n]],[1]);
! 201: }
! 202: Print("--------FW principal parts : ");Println(ans);
! 203: ans = Eliminatev(ans,elim);
! 204: m = Length(ans);
! 205: /* h, z[n] を 1 にする. */
! 206: for (i=0; i<m; i++) {
! 207: ans[i] = Replace(ans[i],[[h,PolyR("1",rz)],[z[n],PolyR("1",rz)]]);
! 208: }
! 209: Println(" "); Println(" ");
! 210: return(ans);
! 211: }
! 212:
! 213: def zrho(f,n) {
! 214: local ans,i,top,w,rz;
! 215: ans = 0;
! 216: rz = Ringp(f);
! 217: while(true) {
! 218: if ( f == Poly("0")) sm1(" exit ");
! 219: top = Init(f);
! 220: f = f-top;
! 221: w = Exponent(top,[Dz[n-1]]);
! 222: top = Replace(top,[[Dz[n-1],PolyR("1",rz)]])*zipoch(z[n],w[0]);
! 223: ans = ans + top;
! 224: }
! 225: return(ans);
! 226: }
! 227:
! 228: def zipoch(f,w) {
! 229: local ans,i;
! 230: ans = 1;
! 231: for (i=0; i<w; i++) {
! 232: ans = ans*(f-i);
! 233: }
! 234: return(ans);
! 235: }
! 236:
! 237:
! 238: fff = ["z[0]*Dz[0]+z[1]*Dz[1]+z[2]*Dz[2]+z[3]*Dz[3]-z[5]",
! 239: " z[1]*Dz[1] +z[3]*Dz[3]-z[6]",
! 240: " z[2]*Dz[2]+z[3]*Dz[3]+z[7]",
! 241: "z[0]*z[3]-z[1]*z[2]"];
! 242: /* z[4] はつかったらだめよ. */
! 243: /* Println(" type in zindicial0(fff,4,3); "); */
! 244:
! 245: def foo2() {
! 246: local ans,n,i,m,r,tmp;
! 247: ans = zindicial0(fff,4,3); n = 4;
! 248: Println(ans);
! 249: m = Length(ans);
! 250: r = [ ];
! 251: for (i=0; i<m; i++) {
! 252: tmp = ans[i];
! 253: tmp = Replace(tmp,[[z[n-1],Poly("1")]]);
! 254: tmp = zrho(tmp,n);
! 255: Print(i); Print(" : ");Println(tmp);
! 256: r = Append(r,tmp);
! 257: }
! 258: return(r);
! 259: }
! 260:
! 261: /* Println("Type in foo2()."); */
! 262:
! 263: def zindicial(fff,n,mm) {
! 264: local ans,n,i,m,r,tmp;
! 265: ans = zindicial0(fff,n,mm);
! 266: Println(ans);
! 267: m = Length(ans);
! 268: r = [ ];
! 269: Println(AddString(["------ The indicial polynomial along z[",
! 270: ToString(n-1),
! 271: "] = 0 is the minimal degree polynomial of the following",
! 272: "polynomials."]));
! 273: Println(AddString(["z[",ToString(n),"] is equal to s."]));
! 274: for (i=0; i<m; i++) {
! 275: tmp = ans[i];
! 276: tmp = Replace(tmp,[[z[n-1],Poly("1")]]);
! 277: tmp = zrho(tmp,n);
! 278: Print(i); Print(" : ");Println(tmp);
! 279: r = Append(r,tmp);
! 280: }
! 281: Println(" ");
! 282: return(r);
! 283: }
! 284:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to toric0.k CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / k097 / debug