[BACK]Return to complex.k CVS log [TXT][DIR] Up to [local] / OpenXM / src / k097 / lib / restriction

Annotation of OpenXM/src/k097/lib/restriction/complex.k, Revision 1.2

1.2     ! takayama    1: /* $OpenXM: OpenXM/src/k097/lib/restriction/complex.k,v 1.1 2001/01/04 12:29:33 takayama Exp $ */
        !             2: /* Document of this module is at k097/Doc/complex.texi */
1.1       takayama    3:
                      4: load["lib/restriction/restriction.k"];;
                      5:
                      6: def Res_solv(m,d,rng) {
                      7:   local r,rr,ans,ac;
                      8:   ac = Length(Arglist);
                      9:   r = GetRing(Poly("1")); /* Save the current ring. */
                     10:   if (ac < 3) {
                     11:     rng = null;
                     12:     rr = GetRing(m);
                     13:     if (Tag(rr) == 0) rr = GetRing(d);
                     14:     if (Tag(rr) != 0) SetRing(rr);
                     15:   }else{
                     16:     SetRing(rng);
                     17:   }
                     18:   m=DC(m,"polynomial"); d = DC(d,"polynomial");
                     19:
                     20:   if (IsRing(rng)) {
                     21:     sm1(" [ m d rng] res-solv /ans set ");
                     22:   }else{
                     23:     sm1(" [ m d] res-solv /ans set ");
                     24:   }
                     25:
                     26:   SetRing(r);
                     27:   return(ans);
                     28: }
                     29:
                     30: /* m : D^p ---> D^q/jj,   u m = d mod jj */
                     31: def Res_solv2(m,d,jj,rng) {
                     32:   local r,rr,ans,ac,pp,qq,kk;
                     33:
                     34:   ac = Length(Arglist);
                     35:   r = GetRing(Poly("1")); /* Save the current ring. */
                     36:   if (ac < 4) {
                     37:     rng = null;
                     38:     rng = GetRing(m);
                     39:     if (Tag(rng) == 0) rng = GetRing(d);
                     40:   }
                     41:
                     42:   pp = Length(m);
                     43:   if (!IsArray(m[0])) {
                     44:      sm1(" m { [ 2 1 roll ] } map /m set ");
                     45:   }
                     46:   qq = Length(m[0]);
                     47:   if (Length(jj) > 0) {
                     48:     if (!IsArray(jj[0])) {
                     49:        sm1(" jj { [ 2 1 roll ] } map /jj set ");
                     50:     }
                     51:     if (qq != Length(jj[0])) {
                     52:       Error(" Matrix size mismatch in m and jj of Kernel2(m,jj,r).");
                     53:     }
                     54:   }
                     55:   m = Join(m,jj);
                     56:
                     57:   ans = Res_solv(m,d,rng);
                     58:   /* Println(ans); */
                     59:   SetRing(r);
                     60:   return([Firstn(ans[0],pp),ans[1]]);
                     61: }
                     62: /* Res_solv2([x,y],[x^2+y^2],[x]):*/
                     63:
                     64: def Res_solv_h(m,d,rng) {
                     65:   local r,rr,ans,ac;
                     66:   ac = Length(Arglist);
                     67:   r = GetRing(Poly("1")); /* Save the current ring. */
                     68:   if (ac < 3) {
                     69:     rng = null;
                     70:     rr = GetRing(m);
                     71:     if (Tag(rr) == 0) rr = GetRing(d);
                     72:     if (Tag(rr) != 0) SetRing(rr);
                     73:   }else{
                     74:     SetRing(rng);
                     75:   }
                     76:   m=DC(m,"polynomial"); d = DC(d,"polynomial");
                     77:
                     78:   if (IsRing(rng)) {
                     79:     sm1(" [ m d rng] res-solv-h /ans set ");
                     80:   }else{
                     81:     sm1(" [ m d] res-solv-h /ans set ");
                     82:   }
                     83:
                     84:   SetRing(r);
                     85:   return(ans);
                     86: }
                     87:
                     88: /* m : D^p ---> D^q/jj,   u m = d mod jj */
                     89: def Res_solv2_h(m,d,jj,rng) {
                     90:   local r,rr,ans,ac,pp,qq,kk;
                     91:
                     92:   ac = Length(Arglist);
                     93:   r = GetRing(Poly("1")); /* Save the current ring. */
                     94:   if (ac < 4) {
                     95:     rng = null;
                     96:     rng = GetRing(m);
                     97:     if (Tag(rng) == 0) rng = GetRing(d);
                     98:   }
                     99:
                    100:   pp = Length(m);
                    101:   if (!IsArray(m[0])) {
                    102:      sm1(" m { [ 2 1 roll ] } map /m set ");
                    103:   }
                    104:   qq = Length(m[0]);
                    105:   if (Length(jj) > 0) {
                    106:     if (!IsArray(jj[0])) {
                    107:        sm1(" jj { [ 2 1 roll ] } map /jj set ");
                    108:     }
                    109:     if (qq != Length(jj[0])) {
                    110:       Error(" Matrix size mismatch in m and jj of Kernel2(m,jj,r).");
                    111:     }
                    112:   }
                    113:   m = Join(m,jj);
                    114:
                    115:   ans = Res_solv_h(m,d,rng);
                    116:   /* Println(ans); */
                    117:   SetRing(r);
                    118:   return([Firstn(ans[0],pp),ans[1]]);
                    119: }
                    120: /* Res_solv2_h([x,y],[x^2+y^2],[x]):*/
                    121:
                    122: def Getxvars() {
                    123:   local v,n,i,ans,ans2;
                    124:   sm1(" getvNames /v set ");
                    125:   sm1(" [(NN)] system_variable (universalNumber) dc /n set ");
                    126:   ans = [];
                    127:   for (i=1; i<n; i++) {
                    128:     ans = Append(ans,v[i]);
                    129:   }
                    130:   sm1(" ans from_records /ans2 set ");
                    131:   return([ans,ans2]);
                    132: }
                    133:
                    134: /* This works only for D cf. Getxvars(). */
                    135: def Intersection(i1,i2,rng) {
                    136:   local r,rr,ans,n,tt,vv,ac;
                    137:   ac = Length(Arglist);
                    138:   r = GetRing(Poly("1")); /* Save the current ring */
                    139:   if (ac < 3) {
                    140:     rr = GetRing(i1);
                    141:     if (Tag(rr) == 0) rr = GetRing(i2);
                    142:     if (Tag(rr) != 0) SetRing(rr);
                    143:   }else{
                    144:     SetRing(rng);
                    145:   }
                    146:   /*
                    147:   i1=DC(i1,"polynomial"); i2 = DC(i2,"polynomial");
                    148:   */
                    149:   i1 = ReParse(i1); i2 = ReParse(i2);
                    150:
                    151:   vv = Getxvars();
                    152:   vv = vv[1];
                    153:   if (Length(i1) == 0) {
                    154:     ans = i2;
                    155:   }else if (!IsArray(i1[0]))  {
                    156:     sm1(" [i1 i2 vv] intersection /ans set ");
                    157:   }else {
                    158:     n = Length(i1[0]);
                    159:     sm1(" i1  fromVectors  /i1 set ");
                    160:     sm1(" i2  fromVectors  /i2 set ");
                    161:     sm1(" [i1 i2 vv] intersection /tt set ");
                    162:     sm1(" [n (integer) dc tt] toVectors /ans set ");
                    163:   }
                    164:
                    165:   SetRing(r);
                    166:   return(ans);
                    167: }
                    168:
                    169: def Firstn(a,n) {
                    170:   local r,i,j,ans;
                    171:   if (Length(a) == 0) {
                    172:     return([ ]);
                    173:   }
                    174:   if (!IsArray(a[0])) {
                    175:     r = NewArray(n);
                    176:     for (i=0; i<n; i++) {
                    177:       r[i] = a[i];
                    178:     }
                    179:     return(r);
                    180:   }else{
                    181:     ans = [ ];
                    182:     for (j=0; j < Length(a); j++) {
                    183:       r = NewArray(n);
                    184:       for (i=0; i<n; i++) {
                    185:         r[i] = a[j,i];
                    186:       }
                    187:       ans = Append(ans,r);
                    188:     }
                    189:     return(ans);
                    190:   }
                    191: }
                    192:
                    193: /* Kernel is also defined in lib/minimal/minimal.k */
                    194: def Kernel(f,v) {
                    195:   local ans;
                    196:   /* v :  string or ring */
                    197:   if (Length(Arglist) < 2) {
                    198:     sm1(" [f] syz /ans set ");
                    199:   }else{
                    200:     sm1(" [f v] syz /ans set ");
                    201:   }
                    202:   return(ans);
                    203: }
                    204:
                    205: def Kernel_h(f,v) {
                    206:   local ans;
                    207:   /* v :  string or ring */
                    208:   if (Length(Arglist) < 2) {
                    209:     sm1(" [f] syz_h /ans set ");
                    210:   }else{
                    211:     sm1(" [f v] syz_h /ans set ");
                    212:   }
                    213:   return(ans);
                    214: }
                    215:
                    216: /*  Kernel of (D^p --- m ---> D^q/jj)  */
                    217: def Kernel2(m,jj,r)
                    218: {
                    219:   local crng,ac,pp,qq,kk;
                    220:   ac = Length(Arglist);
                    221:   crng = GetRing(Poly("1"));
                    222:   if (ac < 3) {
                    223:     r = GetRing(m);
                    224:   }
                    225:   pp = Length(m);
                    226:   if (!IsArray(m[0])) {
                    227:      sm1(" m { [ 2 1 roll ] } map /m set ");
                    228:   }
                    229:   qq = Length(m[0]);
                    230:   if (Length(jj) > 0) {
                    231:     if (!IsArray(jj[0])) {
                    232:        sm1(" jj { [ 2 1 roll ] } map /jj set ");
                    233:     }
                    234:     if (qq != Length(jj[0])) {
                    235:       Error(" Matrix size mismatch in m and jj of Kernel2(m,jj,r).");
                    236:     }
                    237:   }
                    238:   m = Join(m,jj);
                    239:   kk = Kernel(m,r);
                    240:   SetRing(crng);
                    241:   return(Firstn(kk[0],pp));
                    242: }
                    243:
                    244: /*  Kernel of (D^p --- m ---> D^q/jj)  */
                    245: def Kernel2_h(m,jj,r)
                    246: {
                    247:   local crng,ac,pp,qq,kk;
                    248:   ac = Length(Arglist);
                    249:   crng = GetRing(Poly("1"));
                    250:   if (ac < 3) {
                    251:     r = GetRing(m);
                    252:   }
                    253:   pp = Length(m);
                    254:   if (!IsArray(m[0])) {
                    255:      sm1(" m { [ 2 1 roll ] } map /m set ");
                    256:   }
                    257:   qq = Length(m[0]);
                    258:   if (Length(jj) > 0) {
                    259:     if (!IsArray(jj[0])) {
                    260:        sm1(" jj { [ 2 1 roll ] } map /jj set ");
                    261:     }
                    262:     if (qq != Length(jj[0])) {
                    263:       Error(" Matrix size mismatch in m and jj of Kernel2(m,jj,r).");
                    264:     }
                    265:   }
                    266:   m = Join(m,jj);
                    267:   kk = Kernel_h(m,r);
                    268:   SetRing(crng);
                    269:   return(Firstn(kk[0],pp));
                    270: }
                    271:
                    272: /* From lib/minimal/minimal.k */
                    273: def Gb(m,rng) {
                    274:   local r,rr,ans,ac;
                    275:   ac = Length(Arglist);
                    276:   r = GetRing(Poly("1")); /* Save the current ring. */
                    277:   if (ac < 2) {
                    278:     rng = null;
                    279:     rr = GetRing(m);
                    280:     if (Tag(rr) != 0) SetRing(rr);
                    281:   }else{
                    282:     rr = rng;
                    283:     SetRing(rr);
                    284:   }
                    285:   /* m=DC(m,"polynomial");  */
                    286:   m = ReParse(m);
                    287:   sm1(" [m rr] gb /ans set ");
                    288:   SetRing(r);
                    289:   return(ans);
                    290: }
                    291:
                    292: def Gb_h(m,rng) {
                    293:   local r,rr,ans,ac;
                    294:   ac = Length(Arglist);
                    295:   r = GetRing(Poly("1")); /* Save the current ring. */
                    296:   if (ac < 2) {
                    297:     rng = null;
                    298:     rr = GetRing(m);
                    299:     if (Tag(rr) != 0) SetRing(rr);
                    300:   }else{
                    301:     rr = rng;
                    302:     SetRing(rr);
                    303:   }
                    304:   /* m=DC(m,"polynomial");  */
                    305:   m = ReParse(m);
                    306:   sm1(" [m rr] gb_h /ans set ");
                    307:   SetRing(r);
                    308:   return(ans);
                    309: }
                    310:
                    311: def Res_shiftMatrix(m,v,rng) {
                    312:   local n,ans,r,ac,i,j,b1,b2;
                    313:   sm1(" 40 (string) dc /b1 set ");
                    314:   sm1(" 41 (string) dc /b2 set ");
                    315:   ac = Length(Arglist);
                    316:   r = GetRing(Poly("1")); /* Save the current ring. */
                    317:   if (ac < 3) {
                    318:   }else{
                    319:     SetRing(rng);
                    320:   }
                    321:   n = Length(m);
                    322:   ans = NewVector(n);
                    323:   for (i=0; i<n; i++) {
                    324:     ans[i] = NewVector(n);
                    325:     for (j=0; j<n; j++) {
                    326:       ans[i,j] = Poly("0");
                    327:     }
                    328:     ans[i,i] = Poly(AddString([
                    329:         DC(v,"string"),"^",b1,DC(m[i],"string"),b2]));
                    330:   }
                    331:   SetRing(r);
                    332:   return(ans);
                    333: }
                    334:
1.2     ! takayama  335: def ChangeRing(f) {
        !           336:   local r;
        !           337:   r = GetRing(f);
        !           338:   if (Tag(r) == 14) {
        !           339:     SetRing(r);
        !           340:     return(true);
1.1       takayama  341:   }else{
1.2     ! takayama  342:     return(false);
1.1       takayama  343:   }
                    344: }
                    345:
                    346:
                    347: def test2() {
                    348:   RingD("x,y,z");
                    349:   /* Step 1. J */
                    350:   Println(" ---------- J --------------");
                    351:   mm = [[1],
                    352:         [x*Dx],
                    353:         [y*Dy],
                    354:         [z*Dz]];
                    355:   b1 = Res_solv(mm,[Dz]);
                    356:   Println(b1);
                    357:   b2 = Res_solv(mm,[Dy]);
                    358:   Println(b2);
                    359:   /* Step 2. K */
                    360:   Println(" --------- K -------------");
                    361:   mm2 = [[Dz],
                    362:          [Dy],
                    363:         [x*Dx],
                    364:         [y*Dy],
                    365:         [z*Dz]];
                    366:    k1 = Kernel(mm2);
                    367:    Pmat(Firstn(k1[0],2));
                    368:    aa=Kernel2([[Dz],[Dy]],[x*Dx,y*Dy,z*Dz]);
                    369:    Pmat(aa);
                    370: }
                    371:
                    372: def test3() {
                    373:    RingD("x,y,z");
                    374:    mm = [[-Dz],[Dy],[x*Dx],[0]];
                    375:    kk = Kernel(mm);
                    376:    Pmat(kk[0]);
                    377:    rrr= RingD("x,y,z,t",[["t",1,"x",-1,"y",-1,"z",-1,
                    378:                                 "Dx",1,"Dy",1,"Dz",1]]);
                    379:    kk0 = Reparse(kk[0]);
                    380:    gg = Gb(kk0*Res_shiftMatrix([1,1,0,2],t));
                    381:    Pmat(gg[0]); Pmat(gg[1]);
                    382:    gg2= Substitute(Substitute(gg[0],t,1),h,1);
                    383:    Pmat(gg2);
                    384:    Println("----------------------------");
                    385:    mm2 = [[0,0],
                    386:           [0,0],
                    387:           [0,0],
                    388:           [-Dy,-Dz]];
                    389:    jj2 = [[x*Dx,0],
                    390:           [y*Dy,0],
                    391:           [z,0],
                    392:           [0,x*Dx],
                    393:           [0,y],
                    394:           [0,z*Dz]];
                    395:    kk2 = Kernel2(mm2,jj2);
                    396:    Pmat(kk2);
                    397:    Println("-----------------------");
                    398:    ii = Intersection(gg2,kk2);
                    399:    Pmat(ii);
                    400:    SetRing(rrr);
                    401:    ii = Reparse(ii);
                    402:    gg3 = Gb(ii*Res_shiftMatrix([1,1,0,2],t));
                    403:    Pmat(gg3[0]);
                    404:    gg4= Substitute(Substitute(gg3[0],t,1),h,1);
                    405:    Println("---- page 20, observation 4 -----");
                    406:    Pmat(gg4);
                    407: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>