[BACK]Return to ahg.k CVS log [TXT][DIR] Up to [local] / OpenXM / src / k097 / debug

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>