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

Annotation of OpenXM/src/k097/lib/minimal/minimal.k, Revision 1.16

1.16    ! takayama    1: /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.15 2000/06/14 07:44:05 takayama Exp $ */
1.1       takayama    2: #define DEBUG 1
                      3: /* #define ORDINARY 1 */
1.4       takayama    4: /* If you run this program on openxm version 1.1.2 (FreeBSD),
                      5:    make a symbolic link by the command
                      6:    ln -s /usr/bin/cpp /lib/cpp
                      7: */
1.6       takayama    8: #define OFFSET 0
1.14      takayama    9: #define TOTAL_STRATEGY 1
1.6       takayama   10: /* #define OFFSET 20*/
1.1       takayama   11: /* Test sequences.
                     12:    Use load["minimal.k"];;
                     13:
                     14:    a=Sminimal(v);
                     15:    b=a[0];
                     16:    b[1]*b[0]:
                     17:    b[2]*b[1]:
                     18:
                     19:    a = test0();
                     20:    b = a[0];
                     21:    b[1]*b[0]:
                     22:    b[2]*b[1]:
                     23:    a = Sminimal(b[0]);
                     24:
                     25:    a = test1();
                     26:    b=a[0];
                     27:    b[1]*b[0]:
                     28:    b[2]*b[1]:
                     29:
                     30: */
                     31:
                     32:
                     33: load("cohom.k");
                     34: def load_tower() {
                     35:   if (Boundp("k0-tower.sm1.loaded")) {
                     36:   }else{
                     37:     sm1(" [(parse) (k0-tower.sm1) pushfile ] extension ");
                     38:     sm1(" /k0-tower.sm1.loaded 1 def ");
                     39:   }
1.7       takayama   40:   sm1(" oxNoX ");
1.1       takayama   41: }
                     42: load_tower();
                     43: SonAutoReduce = true;
                     44: def Factor(f) {
                     45:    sm1(f, " fctr /FunctionValue set");
                     46: }
                     47: def Reverse(f) {
                     48:    sm1(f," reverse /FunctionValue set");
                     49: }
                     50: def Sgroebner(f) {
                     51:    sm1(" [f] groebner /FunctionValue set");
                     52: }
                     53: def test0() {
                     54:   local f;
                     55:   Sweyl("x,y,z");
                     56:   f = [x^2+y^2+z^2, x*y+x*z+y*z, x*z^2+y*z^2, y^3-x^2*z - x*y*z+y*z^2,
                     57:        -y^2*z^2 + x*z^3 + y*z^3, -z^4];
                     58:   frame=SresolutionFrame(f);
                     59:   Println(frame);
                     60:   /* return(frame); */
                     61:   return(SlaScala(f));
                     62: }
                     63: def test1() {
                     64:   local f;
                     65:   Sweyl("x,y,z");
                     66:   f = [x^2+y^2+z^2, x*y+x*z+y*z, x*z^2+y*z^2, y^3-x^2*z - x*y*z+y*z^2,
                     67:        -y^2*z^2 + x*z^3 + y*z^3, -z^4];
                     68:   return(Sminimal(f));
                     69: }
                     70:
                     71:
                     72:
                     73: def Sweyl(v,w) {
                     74:   /* extern WeightOfSweyl ; */
                     75:   local ww,i,n;
                     76:   if(Length(Arglist) == 1) {
                     77:     sm1(" [v s_ring_of_differential_operators 0 [(schreyer) 1]] define_ring ");
                     78:     sm1(" define_ring_variables ");
                     79:
                     80:     sm1(" [ v to_records pop ] /ww set ");
                     81:     n = Length(ww);
                     82:     WeightOfSweyl = NewArray(n*4);
                     83:     for (i=0; i< n; i++) {
                     84:       WeightOfSweyl[2*i] = ww[i];
                     85:       WeightOfSweyl[2*i+1] = 1;
                     86:     }
                     87:     for (i=0; i< n; i++) {
                     88:       WeightOfSweyl[2*n+2*i] = AddString(["D",ww[i]]);
                     89:       WeightOfSweyl[2*n+2*i+1] = 1;
                     90:     }
                     91:
                     92:   }else{
                     93:     sm1(" [v s_ring_of_differential_operators w s_weight_vector 0 [(schreyer) 1]] define_ring ");
                     94:     sm1(" define_ring_variables ");
                     95:     WeightOfSweyl = w[0];
                     96:   }
                     97: }
                     98:
                     99:
                    100: def Spoly(f) {
                    101:   sm1(f, " toString tparse /FunctionValue set ");
                    102: }
                    103:
                    104: def SreplaceZeroByZeroPoly(f) {
                    105:   if (IsArray(f)) {
                    106:      return(Map(f,"SreplaceZeroByZeroPoly"));
                    107:   }else{
                    108:      if (IsInteger(f)) {
                    109:        return(Poly(ToString(f)));
                    110:      }else{
                    111:        return(f);
                    112:      }
                    113:   }
                    114: }
                    115: def Shomogenize(f) {
                    116:   f = SreplaceZeroByZeroPoly(f);
                    117:   if (IsArray(f)) {
                    118:     sm1(f," sHomogenize2  /FunctionValue set ");
                    119:     /* sm1(f," {sHomogenize2} map  /FunctionValue set ");  */
                    120:     /* Is it correct? Double check.*/
                    121:   }else{
                    122:     sm1(f, " sHomogenize /FunctionValue set ");
                    123:   }
                    124: }
                    125:
                    126: def StoTower() {
                    127:   sm1("  [(AvoidTheSameRing)] pushEnv [ [(AvoidTheSameRing) 0] system_variable (mmLarger) (tower) switch_function ] pop popEnv ");
                    128: }
                    129:
                    130: def SsetTower(tower) {
                    131: sm1(" [(AvoidTheSameRing)] pushEnv
                    132:       [ [(AvoidTheSameRing) 0] system_variable
                    133:         [(gbListTower) tower (list) dc] system_variable
                    134:       ] pop popEnv ");
1.14      takayama  135:       /* sm1("(hoge) message show_ring "); */
1.1       takayama  136: }
                    137:
                    138: def SresolutionFrameWithTower(g,opt) {
                    139:   local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
1.15      takayama  140:         gbasis, nohomog;
                    141:   nohomog = false;
                    142:   count = -1;
1.1       takayama  143:   if (Length(Arglist) >= 2) {
1.15      takayama  144:     if (IsInteger(opt)) {
                    145:       count = opt;
                    146:     }else if (IsString(opt)) {
                    147:       if (opt == "homogenized") {
                    148:          nohomog = true;
                    149:       }else{
                    150:          Println("Warning: unknown option");
                    151:          Println(opt);
                    152:       }
                    153:     }
1.1       takayama  154:   }else{
                    155:     count = -1;
                    156:   }
                    157:
                    158:   sm1(" setupEnvForResolution ");
                    159:   /* If I do not put this macro, homogenization
                    160:      make a strange behavior. For example,
                    161:      [(2*x*Dx + 3*y*Dy+6) (0)] homogenize returns
                    162:      [(2*x*Dx*h + 3*y*Dy*h+6*h^3) (0)].
                    163:      4/19, 2000.
                    164:   */
                    165:
                    166:   sm1(" (mmLarger) (matrix) switch_function ");
1.15      takayama  167:   if (! nohomog) {
                    168:     Println("Automatic homogenization.");
                    169:     g = Map(g,"Shomogenize");
                    170:   }else{
                    171:     Println("No automatic homogenization.");
                    172:   }
1.1       takayama  173:   if (SonAutoReduce) {
                    174:     sm1("[ (AutoReduce) ] system_variable /autof set ");
                    175:     sm1("[ (AutoReduce) 1 ] system_variable ");
                    176:   }
                    177:   gbasis = Sgroebner(g);
                    178:   g = gbasis[0];
                    179:   if (SonAutoReduce) {
                    180:     sm1("[ (AutoReduce) autof] system_variable  ");
                    181:   }
                    182:
                    183:   g = Init(g);
                    184:
                    185: /*  sm1(" setupEnvForResolution-sugar "); */
                    186:   /* -sugar is fine? */
                    187:   sm1(" setupEnvForResolution ");
                    188:
                    189:   Println(g);
                    190:   startingGB = g;
                    191:   /* ans = [ SzeroMap(g) ];  It has not been implemented. see resol1.withZeroMap */
                    192:   ans = [ ];
                    193:   gbTower = [ ];
                    194:   skelton = [ ];
                    195:   while (true) {
                    196:     /* sm1(g," res0Frame /ff set "); */
                    197:     withSkel = Sres0FrameWithSkelton(g);
                    198:     ff = withSkel[0];
                    199:     ans = Append(ans, ff[0]);
                    200:     gbTower = Join([ ff[1] ], gbTower);
                    201:     skelton = Join([ withSkel[1] ], skelton);
                    202:     g = ff[0];
                    203:     if (Length(g) == 0) break;
                    204:     SsetTower( gbTower );
                    205:     if (count == 0) break;
                    206:     count = count - 1;
                    207:   }
                    208:   return([ans,Reverse(gbTower),Join([ [ ] ], Reverse(skelton)),gbasis]);
                    209: }
                    210: HelpAdd(["SresolutionFrameWithTower",
                    211: ["It returs [resolution of the initial, gbTower, skelton, gbasis]",
1.15      takayama  212:  "option: \"homogenized\" (no automatic homogenization) ",
1.1       takayama  213:  "Example: Sweyl(\"x,y\");",
                    214:  "         a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
                    215:
                    216: def SresolutionFrame(f,opt) {
                    217:   local ans;
1.15      takayama  218:   ans = SresolutionFrameWithTower(f,opt);
1.1       takayama  219:   return(ans[0]);
                    220: }
                    221: /* ---------------------------- */
                    222: def ToGradedPolySet(g) {
                    223:   sm1(g," (gradedPolySet) dc /FunctionValue set ");
                    224: }
                    225:
                    226: def NewPolynomialVector(size) {
                    227:   sm1(size," (integer) dc newPolyVector /FunctionValue set ");
                    228: }
                    229:
                    230: def  SturnOffHomogenization() {
                    231:   sm1("
                    232:     [(Homogenize)] system_variable 1 eq
                    233:     { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
                    234:       [(Homogenize) 0] system_variable
                    235:       [(ReduceLowerTerms) 0] system_variable
                    236:     } {  } ifelse
                    237:   ");
                    238: }
                    239: def  SturnOnHomogenization() {
                    240:   sm1("
                    241:     [(Homogenize)] system_variable 0 eq
                    242:     { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message
                    243:       [(Homogenize) 1] system_variable
                    244:       [(ReduceLowerTerms) 1] system_variable
                    245:     } {  } ifelse
                    246:   ");
                    247: }
                    248:
                    249: def SschreyerSkelton(g) {
                    250:   sm1(" [(schreyerSkelton) g] gbext /FunctionValue set ");
                    251: }
                    252: def Stoes(g) {
                    253:   if (IsArray(g)) {
                    254:     sm1(g," {toes} map /FunctionValue set ");
                    255:   }else{
                    256:     sm1(g," toes /FunctionValue set ");
                    257:   }
                    258: }
                    259: def Stoes_vec(g) {
                    260:     sm1(g," toes /FunctionValue set ");
                    261: }
                    262:
                    263: def Sres0Frame(g) {
                    264:   local ans;
                    265:   ans = Sres0FrameWithSkelton(g);
                    266:   return(ans[0]);
                    267: }
                    268: def Sres0FrameWithSkelton(g) {
                    269:   local t_syz, nexttower, m, t_gb, skel, betti,
                    270:         gg, k, i, j, pair, tmp, si, sj, grG, syzAll, gLength;
                    271:
                    272:   SturnOffHomogenization();
                    273:
                    274:   g = Stoes(g);
                    275:   skel = SschreyerSkelton(g);
                    276:   /* Print("Skelton is ");
                    277:   sm1_pmat(skel); */
                    278:   betti = Length(skel);
                    279:
                    280:   gLength = Length(g);
                    281:   grG = ToGradedPolySet(g);
                    282:   syzAll = NewPolynomialVector(betti);
                    283:   for (k=0; k<betti; k++) {
                    284:     pair = skel[k];
                    285:     i = pair[0,0];
                    286:     j = pair[0,1];
                    287:     si = pair[1,0];
                    288:     sj = pair[1,1];
                    289:     /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
                    290:     Print(".");
                    291:
                    292:     t_syz = NewPolynomialVector(gLength);
                    293:     t_syz[i] = si;
                    294:     t_syz[j] = sj;
                    295:     syzAll[k] = t_syz;
                    296:   }
                    297:   t_syz = syzAll;
                    298:   Print("Done. betti="); Println(betti);
                    299:   /* Println(g);  g is in a format such as
                    300:     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
                    301:     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
                    302:     [y-es*x , 3*es^4*y*Dy-es^5*x , 3*es^5*y*Dy-es^6*x , ...]
                    303:     [3*es^3*y*Dy-es^5*x ]
                    304:   */
                    305:   nexttower = Init(g);
                    306:   SturnOnHomogenization();
                    307:   return([[t_syz, nexttower],skel]);
                    308: }
                    309:
                    310:
                    311: def StotalDegree(f) {
1.14      takayama  312:   local d0;
                    313:   sm1(" [(grade) f] gbext (universalNumber) dc /d0 set ");
                    314:   /* Print("degree of "); Print(f); Print(" is "); Println(d0); */
                    315:   return(d0);
1.1       takayama  316: }
                    317:
                    318: /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
                    319: def Sord_w(f,w) {
                    320:   local neww,i,n;
                    321:   n = Length(w);
                    322:   neww = NewArray(n);
                    323:   for (i=0; i<n; i=i+2) {
                    324:     neww[i] = ToString(w[i]);
                    325:   }
                    326:   for (i=1; i<n; i=i+2) {
                    327:     neww[i] = IntegerToSm1Integer(w[i]);
                    328:   }
                    329:   sm1(" f neww ord_w (universalNumber) dc /FunctionValue set ");
                    330: }
                    331:
                    332:
                    333: /* This is not satisfactory. */
                    334: def SinitOfArray(f) {
                    335:   local p,pos,top;
                    336:   if (IsArray(f)) {
                    337:      sm1(f," toes init /p set ");
                    338:      sm1(p," (es). degree (universalNumber) dc /pos set ");
                    339:      return([Init(f[pos]),pos]);
                    340:   } else {
                    341:      return(Init(f));
                    342:   }
                    343: }
                    344:
                    345: def test_SinitOfArray() {
                    346:   local f, frame,p,tower,i,j,k;
                    347:   Sweyl("x,y,z");
                    348:   f = [x^2+y^2+z^2, x*y+x*z+y*z, x*z^2+y*z^2, y^3-x^2*z - x*y*z+y*z^2,
                    349:        -y^2*z^2 + x*z^3 + y*z^3, -z^4];
                    350:   p=SresolutionFrameWithTower(f);
                    351:   sm1_pmat(p);
                    352:   sm1_pmat(SgenerateTable(p[1]));
                    353:   return(p);
                    354:   frame = p[0];
                    355:   sm1_pmat(p[1]);
                    356:   sm1_pmat(frame);
                    357:   sm1_pmat(Map(frame[0],"SinitOfArray"));
                    358:   sm1_pmat(Map(frame[1],"SinitOfArray"));
                    359:   return(p);
                    360: }
                    361:
                    362: /* f is assumed to be a monomial with toes. */
                    363: def Sdegree(f,tower,level) {
1.6       takayama  364:   local i,ww, wd;
                    365:   /* extern WeightOfSweyl; */
                    366:   ww = WeightOfSweyl;
1.5       takayama  367:   f = Init(f);
1.1       takayama  368:   if (level <= 1) return(StotalDegree(f));
                    369:   i = Degree(f,es);
1.6       takayama  370: #ifdef TOTAL_STRATEGY
                    371:   return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
                    372: #endif
                    373:   /* Strategy must be compatible with ordering.  */
                    374:   /* Weight vector must be non-negative, too.  */
                    375:   /* See Sdegree, SgenerateTable, reductionTable. */
                    376:   wd = Sord_w(f,ww);
                    377:   return(wd+Sdegree(tower[level-2,i],tower,level-1));
                    378:
1.1       takayama  379: }
                    380:
                    381: def SgenerateTable(tower) {
                    382:   local height, n,i,j, ans, ans_at_each_floor;
1.16    ! takayama  383:
        !           384:   /*
        !           385:   Print("SgenerateTable: tower=");Println(tower);
        !           386:   sm1(" print_switch_status "); */
1.1       takayama  387:   height = Length(tower);
                    388:   ans = NewArray(height);
                    389:   for (i=0; i<height; i++) {
                    390:     n = Length(tower[i]);
                    391:     ans_at_each_floor=NewArray(n);
                    392:     for (j=0; j<n; j++) {
1.6       takayama  393:       ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1)
                    394:                             + OFFSET;
1.1       takayama  395:       /* Println([i,j,ans_at_each_floor[j]]); */
                    396:     }
                    397:     ans[i] = ans_at_each_floor;
                    398:   }
                    399:   return(ans);
                    400: }
                    401: Sweyl("x,y,z");
                    402: v=[[2*x*Dx + 3*y*Dy+6, 0],
                    403:    [3*x^2*Dy + 2*y*Dx, 0],
                    404:    [0,  x^2+y^2],
                    405:    [0,  x*y]];
                    406: /*  SresolutionFrameWithTower(v); */
                    407:
                    408: def SnewArrayOfFormat(p) {
                    409:   if (IsArray(p)) {
                    410:      return(Map(p,"SnewArrayOfFormat"));
                    411:   }else{
                    412:      return(null);
                    413:   }
                    414: }
1.4       takayama  415: def ScopyArray(a) {
                    416:   local n, i,ans;
                    417:   n = Length(a);
                    418:   ans = NewArray(n);
                    419:   for (i=0; i<n; i++) {
                    420:     ans[i] = a[i];
                    421:   }
                    422:   return(ans);
                    423: }
1.1       takayama  424: def SminOfStrategy(a) {
                    425:   local n,i,ans,tt;
                    426:   ans = 100000; /* very big number */
                    427:   if (IsArray(a)) {
                    428:     n = Length(a);
                    429:     for (i=0; i<n; i++) {
                    430:       if (IsArray(a[i])) {
                    431:         tt = SminOfStrategy(a[i]);
                    432:         if (tt < ans) ans = tt;
                    433:       }else{
                    434:         if (a[i] < ans) ans = a[i];
                    435:       }
                    436:     }
                    437:   }else{
                    438:      if (a < ans) ans = a;
                    439:   }
                    440:   return(ans);
                    441: }
                    442: def SmaxOfStrategy(a) {
                    443:   local n,i,ans,tt;
                    444:   ans = -100000; /* very small number */
                    445:   if (IsArray(a)) {
                    446:     n = Length(a);
                    447:     for (i=0; i<n; i++) {
                    448:       if (IsArray(a[i])) {
                    449:         tt = SmaxOfStrategy(a[i]);
                    450:         if (tt > ans) ans = tt;
                    451:       }else{
                    452:         if (a[i] > ans) ans = a[i];
                    453:       }
                    454:     }
                    455:   }else{
                    456:      if (a > ans) ans = a;
                    457:   }
                    458:   return(ans);
                    459: }
                    460:
                    461:
1.15      takayama  462: def SlaScala(g,opt) {
1.1       takayama  463:   local rf, tower, reductionTable, skel, redundantTable, bases,
                    464:         strategy, maxOfStrategy, height, level, n, i,
                    465:         freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
1.4       takayama  466:         redundantTable_ordinary, redundant_seq_ordinary,
                    467:         reductionTable_tmp;
1.1       takayama  468:   /* extern WeightOfSweyl; */
                    469:   ww = WeightOfSweyl;
1.6       takayama  470:   Print("WeightOfSweyl="); Println(WeightOfSweyl);
1.15      takayama  471:   rf = SresolutionFrameWithTower(g,opt);
1.14      takayama  472:   Print("rf="); sm1_pmat(rf);
1.1       takayama  473:   redundant_seq = 1;   redundant_seq_ordinary = 1;
                    474:   tower = rf[1];
1.16    ! takayama  475:
        !           476:   Println("Generating reduction table which gives an order of reduction.");
        !           477:   Print("WeghtOfSweyl="); Println(WeightOfSweyl);
        !           478:   Print("tower"); Println(tower);
1.1       takayama  479:   reductionTable = SgenerateTable(tower);
1.16    ! takayama  480:   Print("reductionTable="); sm1_pmat(reductionTable);
        !           481:
1.1       takayama  482:   skel = rf[2];
                    483:   redundantTable = SnewArrayOfFormat(rf[1]);
                    484:   redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
                    485:   reducer = SnewArrayOfFormat(rf[1]);
                    486:   freeRes = SnewArrayOfFormat(rf[1]);
                    487:   bettiTable = SsetBettiTable(rf[1],g);
                    488:
                    489:   strategy = SminOfStrategy( reductionTable );
                    490:   maxOfStrategy = SmaxOfStrategy( reductionTable );
                    491:   height = Length(reductionTable);
                    492:   while (strategy <= maxOfStrategy) {
                    493:     for (level = 0; level < height; level++) {
                    494:       n = Length(reductionTable[level]);
1.4       takayama  495:       reductionTable_tmp = ScopyArray(reductionTable[level]);
                    496:       while (SthereIs(reductionTable_tmp,strategy)) {
                    497:         i = SnextI(reductionTable_tmp,strategy,redundantTable,
                    498:                    skel,level,freeRes);
                    499:         Println([level,i]);
                    500:         reductionTable_tmp[i] = -200000;
1.1       takayama  501:         if (reductionTable[level,i] == strategy) {
1.16    ! takayama  502:            Print("Processing [level,i]= "); Print([level,i]);
1.1       takayama  503:            Print("   Strategy = "); Println(strategy);
                    504:            if (level == 0) {
                    505:              if (IsNull(redundantTable[level,i])) {
                    506:                bases = freeRes[level];
                    507:                /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
                    508:                pos = SwhereInGB(tower[0,i],rf[3,0]);
                    509:                bases[i] = rf[3,0,pos];
                    510:                redundantTable[level,i] = 0;
                    511:                redundantTable_ordinary[level,i] = 0;
                    512:                freeRes[level] = bases;
                    513:                /* Println(["GB=",i,bases,tower[0,i]]); */
                    514:              }
                    515:            }else{ /* level >= 1 */
                    516:              if (IsNull(redundantTable[level,i])) {
                    517:                bases = freeRes[level];
                    518:                f = SpairAndReduction(skel,level,i,freeRes,tower,ww);
                    519:                if (f[0] != Poly("0")) {
                    520:                   place = f[3];
                    521:                   /* (level-1, place) is the place for f[0],
                    522:                      which is a newly obtained  GB. */
                    523: #ifdef ORDINARY
                    524:                   redundantTable[level-1,place] = redundant_seq;
                    525:                   redundant_seq++;
                    526: #else
                    527:                   if (f[4] > f[5]) {
                    528:                     /* Zero in the gr-module */
                    529:                     Print("v-degree of [org,remainder] = ");
                    530:                     Println([f[4],f[5]]);
                    531:                     Print("[level,i] = "); Println([level,i]);
                    532:                     redundantTable[level-1,place] = 0;
                    533:                   }else{
                    534:                     redundantTable[level-1,place] = redundant_seq;
                    535:                     redundant_seq++;
                    536:                   }
                    537: #endif
                    538:                   redundantTable_ordinary[level-1,place]
                    539:                      =redundant_seq_ordinary;
                    540:                   redundant_seq_ordinary++;
                    541:                   bases[i] = SunitOfFormat(place,f[1])-f[1];  /* syzygy */
                    542:                   redundantTable[level,i] = 0;
                    543:                   redundantTable_ordinary[level,i] = 0;
                    544:                   /* i must be equal to f[2], I think. Double check. */
                    545:                   freeRes[level] = bases;
                    546:                   bases = freeRes[level-1];
                    547:                   bases[place] = f[0];
                    548:                   freeRes[level-1] = bases;
                    549:                   reducer[level-1,place] = f[1];
                    550:                }else{
                    551:                   redundantTable[level,i] = 0;
                    552:                   bases = freeRes[level];
                    553:                   bases[i] = f[1];  /* Put the syzygy. */
                    554:                   freeRes[level] = bases;
                    555:                }
                    556:              }
                    557:            } /* end of level >= 1 */
                    558:         }
                    559:       }
                    560:     }
                    561:     strategy++;
                    562:   }
                    563:   n = Length(freeRes);
                    564:   freeResV = SnewArrayOfFormat(freeRes);
                    565:   for (i=0; i<n; i++) {
                    566:     bases = freeRes[i];
                    567:     bases = Sbases_to_vec(bases,bettiTable[i]);
                    568:     freeResV[i] = bases;
                    569:   }
                    570:   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
                    571: }
1.4       takayama  572:
                    573: def SthereIs(reductionTable_tmp,strategy) {
                    574:   local n,i;
                    575:   n = Length(reductionTable_tmp);
                    576:   for (i=0; i<n; i++) {
                    577:     if (reductionTable_tmp[i] == strategy) {
                    578:       return(true);
                    579:     }
                    580:   }
                    581:   return(false);
                    582: }
                    583:
                    584: def SnextI(reductionTable_tmp,strategy,redundantTable,
                    585:                                   skel,level,freeRes)
                    586: {
                    587:    local ii,n,p,myindex,i,j,bases;
                    588:    n = Length(reductionTable_tmp);
                    589:    if (level == 0) {
                    590:      for (ii=0; ii<n; ii++) {
                    591:        if (reductionTable_tmp[ii] == strategy) {
                    592:           return(ii);
                    593:         }
                    594:       }
                    595:    }else{
                    596:      for (ii=0; ii<n; ii++) {
                    597:        if (reductionTable_tmp[ii] == strategy) {
                    598:          p = skel[level,ii];
                    599:          myindex = p[0];
                    600:          i = myindex[0]; j = myindex[1];
                    601:          bases = freeRes[level-1];
                    602:          if (IsNull(bases[i]) || IsNull(bases[j])) {
                    603:
                    604:          }else{
                    605:            return(ii);
                    606:          }
                    607:        }
                    608:      }
                    609:    }
1.5       takayama  610:    Print("reductionTable_tmp=");
1.4       takayama  611:    Println(reductionTable_tmp);
1.5       takayama  612:    Println("See also reductionTable, strategy, level,i");
1.4       takayama  613:    Error("SnextI: bases[i] or bases[j] is null for all combinations.");
                    614: }
                    615:
                    616:
1.1       takayama  617:
                    618: def SsetBettiTable(freeRes,g) {
                    619:   local level,i, n,bases,ans;
                    620:   ans = NewArray(Length(freeRes)+1);
                    621:   n = Length(freeRes);
                    622:   if (IsArray(g[0])) {
                    623:     ans[0] = Length(g[0]);
                    624:   }else{
                    625:     ans[0] = 1;
                    626:   }
                    627:   for (level=0; level<n; level++) {
                    628:     bases = freeRes[level];
                    629:     if (IsArray(bases)) {
                    630:       ans[level+1] = Length(bases);
                    631:     }else{
                    632:       ans[level+1] = 1;
                    633:     }
                    634:   }
                    635:   return(ans);
                    636: }
                    637:
                    638: def SwhereInGB(f,tower) {
                    639:   local i,n,p,q;
                    640:   n = Length(tower);
                    641:   for (i=0; i<n; i++) {
                    642:     p = MonomialPart(tower[i]);
                    643:     q = MonomialPart(f);
                    644:     if (p == q) return(i);
                    645:   }
                    646:   Println([f,tower]);
                    647:   Error("whereInGB : [f,myset]: f could not be found in the myset.");
                    648: }
                    649: def SunitOfFormat(pos,forms) {
                    650:   local ans,i,n;
                    651:   n = Length(forms);
                    652:   ans = NewArray(n);
                    653:   for (i=0; i<n; i++) {
                    654:     if (i != pos) {
                    655:       ans[i] = Poly("0");
                    656:     }else{
                    657:       ans[i] = Poly("1");
                    658:     }
                    659:   }
                    660:   return(ans);
                    661: }
                    662:
                    663: def Error(s) {
                    664:   sm1(" s error ");
                    665: }
                    666:
                    667: def IsNull(s) {
                    668:   if (Stag(s) == 0) return(true);
                    669:   else return(false);
                    670: }
                    671:
                    672: def StowerOf(tower,level) {
                    673:   local ans,i;
                    674:   ans = [ ];
                    675:   if (level == 0) return([[]]);
                    676:   for (i=0; i<level; i++) {
                    677:     ans = Append(ans,tower[i]);
                    678:   }
                    679:   return(Reverse(ans));
                    680: }
                    681:
                    682: def Sspolynomial(f,g) {
                    683:   if (IsArray(f)) {
                    684:     f = Stoes_vec(f);
                    685:   }
                    686:   if (IsArray(g)) {
                    687:     g = Stoes_vec(g);
                    688:   }
                    689:   sm1("f g spol /FunctionValue set");
                    690: }
                    691:
                    692: def MonomialPart(f) {
                    693:   sm1(" [(lmonom) f] gbext /FunctionValue set ");
                    694: }
                    695:
1.14      takayama  696: /* WARNING:
                    697:   When you use SwhereInTower, you have to change gbList
                    698:   as below. Ofcourse, you should restrore the gbList
                    699:   SsetTower(StowerOf(tower,level));
                    700:   pos = SwhereInTower(syzHead,tower[level]);
                    701: */
1.1       takayama  702: def SwhereInTower(f,tower) {
                    703:   local i,n,p,q;
                    704:   if (f == Poly("0")) return(-1);
                    705:   n = Length(tower);
                    706:   for (i=0; i<n; i++) {
                    707:     p = MonomialPart(tower[i]);
                    708:     q = MonomialPart(f);
                    709:     if (p == q) return(i);
                    710:   }
                    711:   Println([f,tower]);
                    712:   Error("[f,tower]: f could not be found in the tower.");
                    713: }
                    714:
                    715: def Stag(f) {
                    716:   sm1(f," tag (universalNumber) dc /FunctionValue set");
                    717: }
                    718:
                    719: def SpairAndReduction(skel,level,ii,freeRes,tower,ww) {
                    720:   local i, j, myindex, p, bases, tower2, gi, gj,
                    721:        si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
                    722:        vdeg,vdeg_reduced;
                    723:   Println("SpairAndReduction:");
                    724:
                    725:   if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
                    726:   p = skel[level,ii];
                    727:   myindex = p[0];
                    728:   i = myindex[0]; j = myindex[1];
                    729:   bases = freeRes[level-1];
                    730:   Println(["p and bases ",p,bases]);
                    731:   if (IsNull(bases[i]) || IsNull(bases[j])) {
                    732:     Println([level,i,j,bases[i],bases[j]]);
                    733:     Error("level, i, j : bases[i], bases[j]  must not be NULL.");
                    734:   }
                    735:
                    736:   tower2 = StowerOf(tower,level-1);
                    737:   SsetTower(tower2);
1.14      takayama  738:   Println(["level=",level]);
                    739:   Println(["tower2=",tower2]);
1.1       takayama  740:   /** sm1(" show_ring ");   */
                    741:
                    742:   gi = Stoes_vec(bases[i]);
                    743:   gj = Stoes_vec(bases[j]);
                    744:
                    745:   ssp = Sspolynomial(gi,gj);
                    746:   si = ssp[0,0];
                    747:   sj = ssp[0,1];
                    748:   syzHead = si*es^i;
                    749:   /* This will be the head term, I think. But, double check. */
                    750:   Println([si*es^i,sj*es^j]);
                    751:
                    752:   Print("[gi, gj] = "); Println([gi,gj]);
                    753:   sm1(" [(Homogenize)] system_variable message ");
                    754:   Print("Reduce the element "); Println(si*gi+sj*gj);
                    755:   Print("by  "); Println(bases);
                    756:
                    757:   tmp = Sreduction(si*gi+sj*gj, bases);
                    758:
                    759:   Print("result is "); Println(tmp);
                    760:
1.3       takayama  761:   /* This is essential part for V-minimal resolution. */
                    762:   /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
                    763:   vdeg = SvDegree(si*gi,tower,level-1,ww);
1.1       takayama  764:   vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
                    765:   Print("vdegree of the original = "); Println(vdeg);
                    766:   Print("vdegree of the remainder = "); Println(vdeg_reduced);
                    767:
                    768:   t_syz = tmp[2];
                    769:   si = si*tmp[1]+t_syz[i];
                    770:   sj = sj*tmp[1]+t_syz[j];
                    771:   t_syz[i] = si;
                    772:   t_syz[j] = sj;
1.14      takayama  773:
                    774:   SsetTower(StowerOf(tower,level));
1.1       takayama  775:   pos = SwhereInTower(syzHead,tower[level]);
1.14      takayama  776:
                    777:   SsetTower(StowerOf(tower,level-1));
1.1       takayama  778:   pos2 = SwhereInTower(tmp[0],tower[level-1]);
                    779:   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
                    780:   /* pos is the place to put syzygy at level. */
                    781:   /* pos2 is the place to put a new GB at level-1. */
                    782:   Println(ans);
                    783:   return(ans);
                    784: }
                    785:
                    786: def Sreduction(f,myset) {
                    787:   local n, indexTable, set2, i, j, tmp, t_syz;
                    788:   n = Length(myset);
                    789:   indexTable = NewArray(n);
                    790:   set2 = [ ];
                    791:   j = 0;
                    792:   for (i=0; i<n; i++) {
                    793:     if (IsNull(myset[i])) {
                    794:       indexTable[i] = -1;
                    795: /*    }else if (myset[i] == Poly("0")) {
                    796:       indexTable[i] = -1;  */
                    797:     }else{
                    798:       set2 = Append(set2,Stoes_vec(myset[i]));
                    799:       indexTable[i] = j;
                    800:       j++;
                    801:     }
                    802:   }
                    803:   sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
                    804:   t_syz = NewArray(n);
                    805:   for (i=0; i<n; i++) {
                    806:     if (indexTable[i] != -1) {
                    807:       t_syz[i] = tmp[2, indexTable[i]];
                    808:     }else{
                    809:       t_syz[i] = Poly("0");
                    810:     }
                    811:   }
                    812:   return([tmp[0],tmp[1],t_syz]);
                    813: }
                    814:
                    815: def Warning(s) {
                    816:   Print("Warning: ");
                    817:   Println(s);
                    818: }
                    819: def RingOf(f) {
                    820:   local r;
                    821:   if (IsPolynomial(f)) {
                    822:     if (f != Poly("0")) {
                    823:       sm1(f," (ring) dc /r set ");
                    824:     }else{
                    825:       sm1(" [(CurrentRingp)] system_variable /r set ");
                    826:     }
                    827:   }else{
                    828:     Warning("RingOf(f): the argument f must be a polynomial. Return the current ring.");
                    829:     sm1(" [(CurrentRingp)] system_variable /r set ");
                    830:   }
                    831:   return(r);
                    832: }
                    833:
                    834: def Sfrom_es(f,size) {
                    835:   local c,ans, i, d, myes, myee, j,n,r,ans2;
                    836:   if (Length(Arglist) < 2) size = -1;
                    837:   if (IsArray(f)) return(f);
                    838:   r = RingOf(f);
                    839:   myes = PolyR("es",r);
                    840:   myee = PolyR("e_",r);
                    841:   if (Degree(f,myee) > 0 && size == -1) {
                    842:     if (size == -1) {
                    843:        sm1(f," (array) dc /ans set");
                    844:        return(ans);
                    845:     }
                    846:   }
                    847:
                    848: /*
                    849:     Coefficients(x^2-1,x):
                    850:     [    [    2 , 0 ]  , [    1 , -1 ]  ]
                    851: */
                    852:   if (Degree(f,myee) > 0) {
                    853:     c = Coefficients(f,myee);
                    854:   }else{
                    855:     c = Coefficients(f,myes);
                    856:   }
                    857:   if (size < 0) {
                    858:     size = c[0,0]+1;
                    859:   }
                    860:   ans = NewArray(size);
                    861:   for (i=0; i<size; i++) {ans[i] = 0;}
                    862:   n = Length(c[0]);
                    863:   for (j=0; j<n; j++) {
                    864:     d = c[0,j];
                    865:     ans[d] = c[1,j];
                    866:   }
                    867:   return(ans);
                    868: }
                    869:
                    870: def Sbases_to_vec(bases,size) {
                    871:   local n, giveSize, newbases,i;
                    872:   /*  bases = [1+es*x, [1,2,3*x]] */
                    873:   if (Length(Arglist) > 1) {
                    874:     giveSize = true;
                    875:   }else{
                    876:     giveSize = false;
                    877:   }
                    878:   n = Length(bases);
                    879:   newbases = NewArray(n);
                    880:   for (i=0; i<n; i++) {
                    881:      if (giveSize) {
                    882:        newbases[i] = Sfrom_es(bases[i], size);
                    883:      }else{
                    884:        newbases[i] = Sfrom_es(bases[i]);
                    885:      }
                    886:   }
                    887:   return(newbases);
                    888: }
                    889:
1.14      takayama  890: HelpAdd(["Sminimal",
                    891: ["It constructs the V-minimal free resolution by LaScala-Stillman's algorithm",
1.15      takayama  892:  "option: \"homogenized\" (no automatic homogenization ",
1.14      takayama  893:  "Example:  Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
                    894:  "          v=[[2*x*Dx + 3*y*Dy+6, 0],",
                    895:  "             [3*x^2*Dy + 2*y*Dx, 0],",
                    896:  "             [0,  x^2+y^2],",
                    897:  "             [0,  x*y]];",
                    898:  "         a=Sminimal(v);",
                    899:  "         Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
                    900:  "         b = ReParse(a[0]); sm1_pmat(b); ",
                    901:  "         IsExact_h(b,[x,y]):",
                    902:  "Note:  a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);
                    903:
1.15      takayama  904: def Sminimal(g,opt) {
1.1       takayama  905:   local r, freeRes, redundantTable, reducer, maxLevel,
                    906:         minRes, seq, maxSeq, level, betti, q, bases, dr,
1.14      takayama  907:         betti_levelplus, newbases, i, j,qq, tminRes;
1.16    ! takayama  908:   if (Length(Arglist) < 2) {
        !           909:      opt = null;
        !           910:   }
        !           911:   ScheckIfSchreyer("Sminimal:0");
1.15      takayama  912:   r = SlaScala(g,opt);
1.1       takayama  913:   /* Should I turn off the tower?? */
1.16    ! takayama  914:   ScheckIfSchreyer("Sminimal:1");
1.1       takayama  915:   freeRes = r[0];
                    916:   redundantTable = r[1];
                    917:   reducer = r[2];
                    918:   minRes = SnewArrayOfFormat(freeRes);
                    919:   seq = 0;
                    920:   maxSeq = SgetMaxSeq(redundantTable);
                    921:   maxLevel = Length(freeRes);
                    922:   for (level = 0; level < maxLevel; level++) {
                    923:     minRes[level] = freeRes[level];
                    924:   }
                    925:   seq=maxSeq+1;
                    926:   while (seq > 1) {
                    927:     seq--;
                    928:     for (level = 0; level < maxLevel; level++) {
                    929:       betti = Length(freeRes[level]);
                    930:       for (q = 0; q<betti; q++) {
                    931:         if (redundantTable[level,q] == seq) {
                    932:           Print("[seq,level,q]="); Println([seq,level,q]);
                    933:           if (level < maxLevel-1) {
                    934:             bases = freeRes[level+1];
                    935:             dr = reducer[level,q];
                    936:             dr[q] = -1;
                    937:             newbases = SnewArrayOfFormat(bases);
                    938:             betti_levelplus = Length(bases);
                    939:             /*
                    940:                bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
                    941:             */
                    942:             for (i=0; i<betti_levelplus; i++) {
                    943:               newbases[i] = bases[i] + bases[i,q]*dr;
                    944:             }
                    945:             Println(["level, q =", level,q]);
                    946:             Println("bases="); sm1_pmat(bases);
                    947:             Println("dr="); sm1_pmat(dr);
                    948:             Println("newbases="); sm1_pmat(newbases);
                    949:             minRes[level+1] = newbases;
                    950:             freeRes = minRes;
                    951: #ifdef DEBUG
                    952:             for (qq=0; qq<betti; qq++) {
                    953:               if ((redundantTable[level,qq] >= seq) &&
                    954:                   (redundantTable[level,qq] <= maxSeq)) {
                    955:                 for (i=0; i<betti_levelplus; i++) {
                    956:                   if (!IsZero(newbases[i,qq])) {
                    957:                     Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
                    958:                     Print("redundantTable ="); sm1_pmat(redundantTable[level]);
                    959:                     Error("Stop in Sminimal for debugging.");
                    960:                   }
                    961:                 }
                    962:               }
                    963:             }
                    964: #endif
                    965:           }
                    966:         }
                    967:       }
                    968:     }
                    969:    }
1.14      takayama  970:    tminRes = Stetris(minRes,redundantTable);
                    971:    return([SpruneZeroRow(tminRes), tminRes,
1.3       takayama  972:           [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
1.1       takayama  973:   /* r[4] is the redundantTable_ordinary */
1.3       takayama  974:   /* r[0] is the freeResolution */
1.1       takayama  975: }
                    976:
                    977:
                    978: def IsZero(f) {
                    979:   if (IsPolynomial(f)) {
                    980:     return( f == Poly("0"));
                    981:   }else if (IsInteger(f)) {
                    982:     return( f == 0);
                    983:   }else if (IsSm1Integer(f)) {
                    984:     return( f == true );
                    985:   }else if (IsDouble(f)) {
                    986:     return( f == 0.0 );
                    987:   }else if (IsRational(f)) {
                    988:     return(IsZero(Denominator(f)));
                    989:   }else{
                    990:     Error("IsZero: cannot deal with this data type.");
                    991:   }
                    992: }
                    993: def SgetMaxSeq(redundantTable) {
                    994:    local level,i,n,ans, levelMax,bases;
                    995:    levelMax = Length( redundantTable );
                    996:    ans = 0;
                    997:    for (level = 0; level < levelMax; level++) {
                    998:      bases = redundantTable[level];
                    999:      n = Length(bases);
                   1000:      for (i=0; i<n; i++) {
                   1001:        if (IsInteger( bases[i] )) {
                   1002:           if (bases[i] > ans) {
                   1003:              ans = bases[i];
                   1004:           }
                   1005:        }
                   1006:      }
                   1007:    }
                   1008:    return(ans);
                   1009: }
                   1010:
                   1011: def Stetris(freeRes,redundantTable) {
                   1012:   local level, i, j, resLength, minRes,
                   1013:         bases, newbases, newbases2;
                   1014:   minRes = SnewArrayOfFormat(freeRes);
                   1015:   resLength = Length( freeRes );
                   1016:   for (level=0; level<resLength; level++) {
                   1017:     bases = freeRes[level];
                   1018:     newbases = SnewArrayOfFormat(bases);
                   1019:     betti = Length(bases); j = 0;
                   1020:     /* Delete rows */
                   1021:     for (i=0; i<betti; i++) {
                   1022:       if (redundantTable[level,i] < 1) {
                   1023:          newbases[j] = bases[i];
                   1024:          j++;
                   1025:       }
                   1026:     }
                   1027:     bases = SfirstN(newbases,j);
                   1028:     if (level > 0) {
                   1029:       /* Delete columns */
                   1030:       newbases = Transpose(bases);
                   1031:       betti = Length(newbases); j = 0;
                   1032:       newbases2 = SnewArrayOfFormat(newbases);
                   1033:       for (i=0; i<betti; i++) {
                   1034:         if (redundantTable[level-1,i] < 1) {
                   1035:            newbases2[j] = newbases[i];
                   1036:            j++;
                   1037:         }
                   1038:       }
                   1039:       newbases = Transpose(SfirstN(newbases2,j));
                   1040:     }else{
                   1041:       newbases = bases;
                   1042:     }
                   1043:     Println(["level=", level]);
                   1044:     sm1_pmat(bases);
                   1045:     sm1_pmat(newbases);
                   1046:
                   1047:     minRes[level] = newbases;
                   1048:   }
                   1049:   return(minRes);
                   1050: }
                   1051:
                   1052: def SfirstN(bases,k) {
                   1053:    local ans,i;
                   1054:    ans = NewArray(k);
                   1055:    for (i=0; i<k; i++) {
                   1056:      ans[i] = bases[i];
                   1057:    }
                   1058:    return(ans);
                   1059: }
                   1060:
                   1061:
                   1062: /* usage:  tt is tower. ww is weight.
                   1063:     a = SresolutionFrameWithTower(v);
                   1064:     tt = a[1];
                   1065:     ww = [x,1,y,1,Dx,1,Dy,1];
                   1066:     SvDegree(x*es,tt,1,ww):
                   1067:
                   1068: In(17)=tt:
                   1069: [[2*x*Dx , e_*x^2 , e_*x*y , 3*x^2*Dy , e_*y^3 , 9*x*y*Dy^2 , 27*y^2*Dy^3 ]  ,
                   1070:  [es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ]  ,
                   1071:  [3*es^3*y*Dy ]  ]
                   1072: In(18)=SvDegree(x*es,tt,1,ww):
                   1073: 3
                   1074: In(19)=SvDegree(x*es^3,tt,1,ww):
                   1075: 4
                   1076: In(20)=SvDegree(x,tt,2,ww):
                   1077: 4
                   1078:
                   1079: */
                   1080: def SvDegree(f,tower,level,w) {
                   1081:   local i,ans;
                   1082:   if (IsZero(f)) return(null);
1.3       takayama 1083:   f = Init(f);
1.1       takayama 1084:   if (level <= 0) {
                   1085:     return(Sord_w(f,w));
                   1086:   }
                   1087:   i = Degree(f,es);
                   1088:   ans = Sord_w(f,w) +
                   1089:         SvDegree(tower[level-1,i],tower,level-1,w);
                   1090:   return(ans);
                   1091: }
                   1092:
1.2       takayama 1093: def Sannfs(f,v) {
                   1094:   local f2;
                   1095:   f2 = ToString(f);
                   1096:   if (IsArray(v)) {
                   1097:      v = Map(v,"ToString");
                   1098:   }
                   1099:   sm1(" [f2 v] annfs /FunctionValue set ");
                   1100: }
                   1101:
                   1102: /* Sannfs2("x^3-y^2"); */
                   1103: def Sannfs2(f) {
                   1104:   local p,pp;
                   1105:   p = Sannfs(f,"x,y");
1.6       takayama 1106:   sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
1.5       takayama 1107: /*
                   1108:   Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
                   1109:                ["x",-1,"y",-1,"Dx",1,"Dy",1]]); */
1.6       takayama 1110:   /* Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1]]); */
1.10      takayama 1111:
1.6       takayama 1112:   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
                   1113:   pp = Map(p,"Spoly");
                   1114:   return(Sminimal_v(pp));
                   1115:   /* return(Sminimal(pp)); */
                   1116: }
                   1117:
1.10      takayama 1118: HelpAdd(["Sannfs2",
                   1119: ["Sannfs2(f) constructs the V-minimal free resolution for the weight (-1,1)",
                   1120:  "of the Laplace transform of the annihilating ideal of the polynomial f in x,y.",
                   1121:  "See also Sminimal_v, Sannfs3.",
                   1122:  "Example: a=Sannfs2(\"x^3-y^2\");",
                   1123:  "         b=a[0]; sm1_pmat(b);",
                   1124:  "         b[1]*b[0]:",
                   1125:  "Example: a=Sannfs2(\"x*y*(x-y)*(x+y)\");",
                   1126:  "         b=a[0]; sm1_pmat(b);",
                   1127:  "         b[1]*b[0]:"
                   1128: ]]);
                   1129:
1.6       takayama 1130: /* Do not forget to turn on TOTAL_STRATEGY */
                   1131: def Sannfs2_laScala(f) {
                   1132:   local p,pp;
                   1133:   p = Sannfs(f,"x,y");
                   1134:   /*   Do not make laplace transform.
                   1135:     sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
                   1136:     p = [p];
                   1137:   */
                   1138:   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
1.2       takayama 1139:   pp = Map(p[0],"Spoly");
                   1140:   return(Sminimal(pp));
                   1141: }
                   1142:
1.11      takayama 1143: def Sannfs2_laScala2(f) {
                   1144:   local p,pp;
                   1145:   p = Sannfs(f,"x,y");
                   1146:   sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
                   1147:   p = [p];
                   1148:   Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
                   1149:                ["x",-1,"y",-1,"Dx",1,"Dy",1]]);
                   1150:   pp = Map(p[0],"Spoly");
                   1151:   return(Sminimal(pp));
                   1152: }
                   1153:
1.3       takayama 1154: def Sannfs3(f) {
                   1155:   local p,pp;
                   1156:   p = Sannfs(f,"x,y,z");
1.6       takayama 1157:   sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
1.3       takayama 1158:   Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
1.6       takayama 1159:   pp = Map(p,"Spoly");
                   1160:   return(Sminimal_v(pp));
1.3       takayama 1161: }
                   1162:
1.10      takayama 1163: HelpAdd(["Sannfs3",
                   1164: ["Sannfs3(f) constructs the V-minimal free resolution for the weight (-1,1)",
                   1165:  "of the Laplace transform of the annihilating ideal of the polynomial f in x,y,z.",
                   1166:  "See also Sminimal_v, Sannfs2.",
                   1167:  "Example: a=Sannfs3(\"x^3-y^2*z^2\");",
                   1168:  "         b=a[0]; sm1_pmat(b);",
                   1169:  "         b[1]*b[0]: b[2]*b[1]:"]]);
                   1170:
1.2       takayama 1171: /*
                   1172:   The betti numbers of most examples are 2,1. (0-th and 1-th).
                   1173:   a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2.
                   1174:   a=Sannfs2("x^3-y^2-x");    : it causes an error. It should be fixed.
1.3       takayama 1175:   a=Sannfs2("x*y*(x-y)");    : it causes an error. It should be fixed.
1.2       takayama 1176:
                   1177: */
                   1178:
1.11      takayama 1179: def Sannfs3_laScala2(f) {
                   1180:   local p,pp;
                   1181:   p = Sannfs(f,"x,y,z");
                   1182:   sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
                   1183:   Sweyl("x,y,z",[["x",1,"y",1,"z",1,"Dx",1,"Dy",1,"Dz",1,"h",1],
                   1184:                  ["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
                   1185:   pp = Map(p,"Spoly");
                   1186:   return(Sminimal(pp));
                   1187: }
1.5       takayama 1188:
                   1189:
1.6       takayama 1190: /*  The below does not use LaScala-Stillman's algorithm. */
1.5       takayama 1191: def Sschreyer(g) {
                   1192:   local rf, tower, reductionTable, skel, redundantTable, bases,
                   1193:         strategy, maxOfStrategy, height, level, n, i,
                   1194:         freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
                   1195:         redundantTable_ordinary, redundant_seq_ordinary,
1.6       takayama 1196:         reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase;
1.5       takayama 1197:   /* extern WeightOfSweyl; */
                   1198:   ww = WeightOfSweyl;
                   1199:   Print("WeghtOfSweyl="); Println(WeightOfSweyl);
                   1200:   rf = SresolutionFrameWithTower(g);
                   1201:   redundant_seq = 1;   redundant_seq_ordinary = 1;
                   1202:   tower = rf[1];
1.16    ! takayama 1203:   Println("Generating reduction table which gives an order of reduction.");
        !          1204:   Println("But, you are in Sschreyer...., you may not use LaScala-Stillman");
        !          1205:   Print("WeghtOfSweyl="); Println(WeightOfSweyl);
        !          1206:   Print("tower"); Println(tower);
1.5       takayama 1207:   reductionTable = SgenerateTable(tower);
                   1208:   skel = rf[2];
                   1209:   redundantTable = SnewArrayOfFormat(rf[1]);
                   1210:   redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
                   1211:   reducer = SnewArrayOfFormat(rf[1]);
                   1212:   freeRes = SnewArrayOfFormat(rf[1]);
                   1213:   bettiTable = SsetBettiTable(rf[1],g);
                   1214:
                   1215:   height = Length(reductionTable);
                   1216:   for (level = 0; level < height; level++) {
                   1217:       n = Length(reductionTable[level]);
                   1218:       for (i=0; i<n; i++) {
                   1219:            Println([level,i]);
                   1220:            Print("Processing "); Print([level,i]);
                   1221:            if (level == 0) {
                   1222:              if (IsNull(redundantTable[level,i])) {
                   1223:                bases = freeRes[level];
                   1224:                /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
                   1225:                pos = SwhereInGB(tower[0,i],rf[3,0]);
                   1226:                bases[i] = rf[3,0,pos];
                   1227:                /* redundantTable[level,i] = 0;
                   1228:                redundantTable_ordinary[level,i] = 0; */
                   1229:                freeRes[level] = bases;
                   1230:                /* Println(["GB=",i,bases,tower[0,i]]); */
                   1231:              }
                   1232:            }else{ /* level >= 1 */
                   1233:              if (IsNull(redundantTable[level,i])) {
                   1234:                bases = freeRes[level];
                   1235:                f = SpairAndReduction2(skel,level,i,freeRes,tower,
                   1236:                                       ww,redundantTable);
                   1237:                if (f[0] != Poly("0")) {
                   1238:                   place = f[3];
                   1239:                   /* (level-1, place) is the place for f[0],
                   1240:                      which is a newly obtained  GB. */
                   1241: #ifdef ORDINARY
                   1242:                   redundantTable[level-1,place] = redundant_seq;
                   1243:                   redundant_seq++;
                   1244: #else
                   1245:                   if (f[4] > f[5]) {
                   1246:                     /* Zero in the gr-module */
                   1247:                     Print("v-degree of [org,remainder] = ");
                   1248:                     Println([f[4],f[5]]);
                   1249:                     Print("[level,i] = "); Println([level,i]);
                   1250:                     redundantTable[level-1,place] = 0;
                   1251:                   }else{
                   1252:                     redundantTable[level-1,place] = redundant_seq;
                   1253:                     redundant_seq++;
                   1254:                   }
                   1255: #endif
                   1256:                   redundantTable_ordinary[level-1,place]
                   1257:                      =redundant_seq_ordinary;
                   1258:                   redundant_seq_ordinary++;
                   1259:                   bases[i] = SunitOfFormat(place,f[1])-f[1];  /* syzygy */
                   1260:                   /* redundantTable[level,i] = 0;
                   1261:                   redundantTable_ordinary[level,i] = 0; */
                   1262:                   /* i must be equal to f[2], I think. Double check. */
                   1263:
                   1264:                   /* Correction Of Constant */
1.9       takayama 1265:                   /* Correction of syzygy */
1.7       takayama 1266:                   c2 = f[6];  /* or -f[6]?  Double check. */
                   1267:                   Print("c2="); Println(c2);
1.5       takayama 1268:                   nn = Length(bases);
                   1269:                   for (ii=0; ii<nn;ii++) {
1.8       takayama 1270:                      if ((ii != i) && (! IsNull(bases[ii]))) {
1.7       takayama 1271:                        m = Length(bases[ii]);
                   1272:                        for (jj=0; jj<m; jj++) {
                   1273:                          if (jj != place) {
                   1274:                            bases[ii,jj] = bases[ii,jj]*c2;
                   1275:                          }
                   1276:                        }
1.5       takayama 1277:                      }
                   1278:                   }
                   1279:
1.7       takayama 1280:                   Print("Old freeRes[level] = "); sm1_pmat(freeRes[level]);
1.5       takayama 1281:                   freeRes[level] = bases;
1.7       takayama 1282:                   Print("New freeRes[level] = "); sm1_pmat(freeRes[level]);
1.6       takayama 1283:
                   1284:                  /* Update the freeRes[level-1] */
1.7       takayama 1285:                   Print("Old freeRes[level-1] = "); sm1_pmat(freeRes[level-1]);
1.6       takayama 1286:                   bases = freeRes[level-1];
                   1287:                   bases[place] = f[0];
                   1288:                   freeRes[level-1] = bases;
1.7       takayama 1289:                   Print("New freeRes[level-1] = "); sm1_pmat(freeRes[level-1]);
1.6       takayama 1290:
1.9       takayama 1291:                   reducer[level-1,place] = f[1]-SunitOfFormat(place,f[1]);
                   1292:                    /* This reducer is different from that of SlaScala(). */
                   1293:
                   1294:                   reducerBasis = reducer[level-1];
                   1295:                   nn = Length(reducerBasis);
                   1296:                   for (ii=0; ii<nn;ii++) {
                   1297:                      if ((ii != place) && (! IsNull(reducerBasis[ii]))) {
                   1298:                        m = Length(reducerBasis[ii]);
                   1299:                        for (jj=0; jj<m; jj++) {
                   1300:                          if (jj != place) {
                   1301:                            reducerBasis[ii,jj] = reducerBasis[ii,jj]*c2;
                   1302:                          }
                   1303:                        }
                   1304:                      }
                   1305:                   }
                   1306:                   reducer[level-1] = reducerBasis;
                   1307:
1.5       takayama 1308:                }else{
                   1309:                   /* redundantTable[level,i] = 0; */
                   1310:                   bases = freeRes[level];
                   1311:                   bases[i] = f[1];  /* Put the syzygy. */
                   1312:                   freeRes[level] = bases;
                   1313:                }
                   1314:              }  /* end of level >= 1 */
                   1315:           }
                   1316:     } /* i loop */
1.6       takayama 1317:
                   1318:     /* Triangulate reducer */
                   1319:     if (level >= 1) {
                   1320:       Println(" ");
                   1321:       Print("Triangulating reducer at level "); Println(level-1);
1.9       takayama 1322:       Println("freeRes[level]="); sm1_pmat(freeRes[level]);
1.6       takayama 1323:       reducerBase = reducer[level-1];
                   1324:       Print("reducerBase=");  Println(reducerBase);
1.9       takayama 1325:       Println("Compare freeRes[level] and reducerBase (put -1)");
1.6       takayama 1326:       m = Length(reducerBase);
                   1327:       for (ii=m-1; ii>=0; ii--) {
                   1328:         if (!IsNull(reducerBase[ii])) {
                   1329:            for (jj=ii-1; jj>=0; jj--) {
                   1330:              if (!IsNull(reducerBase[jj])) {
                   1331:               if (!IsZero(reducerBase[jj,ii])) {
1.9       takayama 1332:                 /* reducerBase[ii,ii] should be always constant. */
                   1333:                 reducerBase[jj] = reducerBase[ii,ii]*reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
1.6       takayama 1334:               }
                   1335:              }
                   1336:            }
                   1337:          }
                   1338:        }
                   1339:        Println("New reducer");
                   1340:        sm1_pmat(reducerBase);
                   1341:        reducer[level-1] = reducerBase;
                   1342:     }
                   1343:
1.5       takayama 1344:   } /* level loop */
                   1345:   n = Length(freeRes);
                   1346:   freeResV = SnewArrayOfFormat(freeRes);
                   1347:   for (i=0; i<n; i++) {
                   1348:     bases = freeRes[i];
                   1349:     bases = Sbases_to_vec(bases,bettiTable[i]);
                   1350:     freeResV[i] = bases;
                   1351:   }
1.6       takayama 1352:
                   1353:   /* Mark the non-redundant elements. */
                   1354:   for (i=0; i<n; i++) {
                   1355:     m = Length(redundantTable[i]);
                   1356:     for (jj=0; jj<m; jj++) {
                   1357:       if (IsNull(redundantTable[i,jj])) {
                   1358:         redundantTable[i,jj] = 0;
                   1359:       }
                   1360:     }
                   1361:   }
                   1362:
                   1363:
1.5       takayama 1364:   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
                   1365: }
                   1366:
                   1367: def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) {
                   1368:   local i, j, myindex, p, bases, tower2, gi, gj,
                   1369:        si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
                   1370:        vdeg,vdeg_reduced,n,c2;
1.6       takayama 1371:   Println("SpairAndReduction2 : -------------------------");
1.5       takayama 1372:
                   1373:   if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
                   1374:   p = skel[level,ii];
                   1375:   myindex = p[0];
                   1376:   i = myindex[0]; j = myindex[1];
                   1377:   bases = freeRes[level-1];
                   1378:   Println(["p and bases ",p,bases]);
                   1379:   if (IsNull(bases[i]) || IsNull(bases[j])) {
                   1380:     Println([level,i,j,bases[i],bases[j]]);
                   1381:     Error("level, i, j : bases[i], bases[j]  must not be NULL.");
                   1382:   }
                   1383:
                   1384:   tower2 = StowerOf(tower,level-1);
                   1385:   SsetTower(tower2);
1.14      takayama 1386:   Println(["level=",level]);
                   1387:   Println(["tower2=",tower2]);
1.5       takayama 1388:   /** sm1(" show_ring ");   */
                   1389:
                   1390:   gi = Stoes_vec(bases[i]);
                   1391:   gj = Stoes_vec(bases[j]);
                   1392:
                   1393:   ssp = Sspolynomial(gi,gj);
                   1394:   si = ssp[0,0];
                   1395:   sj = ssp[0,1];
                   1396:   syzHead = si*es^i;
                   1397:   /* This will be the head term, I think. But, double check. */
                   1398:   Println([si*es^i,sj*es^j]);
                   1399:
                   1400:   Print("[gi, gj] = "); Println([gi,gj]);
                   1401:   sm1(" [(Homogenize)] system_variable message ");
                   1402:   Print("Reduce the element "); Println(si*gi+sj*gj);
                   1403:   Print("by  "); Println(bases);
                   1404:
                   1405:   tmp = Sreduction(si*gi+sj*gj, bases);
                   1406:
                   1407:   Print("result is "); Println(tmp);
1.6       takayama 1408:   if (!IsZero(tmp[0])) {
                   1409:     Print("Error: base = ");
                   1410:     Println(Map(bases,"Stoes_vec"));
                   1411:     Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring.");
                   1412:   }
1.5       takayama 1413:   t_syz = tmp[2];
                   1414:   si = si*tmp[1]+t_syz[i];
                   1415:   sj = sj*tmp[1]+t_syz[j];
                   1416:   t_syz[i] = si;
                   1417:   t_syz[j] = sj;
                   1418:
                   1419:   c2 = null;
                   1420:   /* tmp[0] must be zero */
                   1421:   n = Length(t_syz);
                   1422:   for (i=0; i<n; i++) {
1.6       takayama 1423:      if (IsConstant(t_syz[i])){
                   1424:       if (!IsZero(t_syz[i])) {
1.5       takayama 1425:        if (IsNull(redundantTable[level-1,i])) {
                   1426:          /* i must equal to pos2 below. */
                   1427:          c2 = -t_syz[i];
1.6       takayama 1428:          tmp[0] = c2*Stoes_vec(freeRes[level-1,i]);
1.5       takayama 1429:          t_syz[i] = 0;
1.6       takayama 1430:          /* tmp[0] = t_syz . g */
1.5       takayama 1431:          /* break; does not work. Use */
                   1432:          i = n;
                   1433:        }
1.6       takayama 1434:       }
1.5       takayama 1435:      }
                   1436:   }
                   1437:
                   1438:   /* This is essential part for V-minimal resolution. */
                   1439:   /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
                   1440:   vdeg = SvDegree(si*gi,tower,level-1,ww);
                   1441:   vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
                   1442:   Print("vdegree of the original = "); Println(vdeg);
                   1443:   Print("vdegree of the remainder = "); Println(vdeg_reduced);
                   1444:
1.11      takayama 1445:   if (!IsNull(vdeg_reduced)) {
                   1446:     if (vdeg_reduced < vdeg) {
                   1447:       Println("--- Special in V-minimal!");
                   1448:       Println(tmp[0]);
                   1449:       Println("syzygy="); sm1_pmat(t_syz);
                   1450:       Print("[vdeg, vdeg_reduced] = "); Println([vdeg,vdeg_reduced]);
                   1451:     }
                   1452:   }
                   1453:
1.14      takayama 1454:   SsetTower(StowerOf(tower,level));
                   1455:   pos = SwhereInTower(syzHead,tower[level]);
1.11      takayama 1456:
1.14      takayama 1457:   SsetTower(StowerOf(tower,level-1));
1.5       takayama 1458:   pos2 = SwhereInTower(tmp[0],tower[level-1]);
                   1459:   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
                   1460:   /* pos is the place to put syzygy at level. */
                   1461:   /* pos2 is the place to put a new GB at level-1. */
                   1462:   Println(ans);
1.14      takayama 1463:   Println("--- end of SpairAndReduction2  ");
1.5       takayama 1464:   return(ans);
                   1465: }
1.6       takayama 1466:
1.10      takayama 1467: HelpAdd(["Sminimal_v",
                   1468: ["It constructs the V-minimal free resolution from the Schreyer resolution",
                   1469:  "step by step.",
1.14      takayama 1470:  "This code still contains bugs. It sometimes outputs wrong answer.",
1.10      takayama 1471:  "Example:   Sweyl(\"x,y\",[[\"x\",-1,\"y\",-1,\"Dx\",1,\"Dy\",1]]);",
                   1472:  "          v=[[2*x*Dx + 3*y*Dy+6, 0],",
                   1473:  "             [3*x^2*Dy + 2*y*Dx, 0],",
                   1474:  "             [0,  x^2+y^2],",
                   1475:  "             [0,  x*y]];",
                   1476:  "         a=Sminimal_v(v);",
                   1477:  "         sm1_pmat(a[0]); b=a[0]; b[1]*b[0]:",
                   1478:  "Note:  a[0] is the V-minimal resolution. a[3] is the Schreyer resolution."]]);
                   1479:
1.14      takayama 1480: /* This code still contains bugs. It sometimes outputs wrong answer. */
                   1481: /* See test12() in minimal-test.k.  */
                   1482: /* There may be remaining 1, too */
1.6       takayama 1483: def Sminimal_v(g) {
                   1484:   local r, freeRes, redundantTable, reducer, maxLevel,
                   1485:         minRes, seq, maxSeq, level, betti, q, bases, dr,
1.10      takayama 1486:         betti_levelplus, newbases, i, j,qq,tminRes;
1.6       takayama 1487:   r = Sschreyer(g);
                   1488:   sm1_pmat(r);
                   1489:   Debug_Sminimal_v = r;
                   1490:   Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v");
                   1491:   /* Should I turn off the tower?? */
                   1492:   freeRes = r[0];
                   1493:   redundantTable = r[1];
                   1494:   reducer = r[2];
                   1495:   minRes = SnewArrayOfFormat(freeRes);
                   1496:   seq = 0;
                   1497:   maxSeq = SgetMaxSeq(redundantTable);
                   1498:   maxLevel = Length(freeRes);
                   1499:   for (level = 0; level < maxLevel; level++) {
                   1500:     minRes[level] = freeRes[level];
                   1501:   }
                   1502:   for (level = 0; level < maxLevel; level++) {
                   1503:       betti = Length(freeRes[level]);
                   1504:       for (q = betti-1; q>=0; q--) {
                   1505:         if (redundantTable[level,q] > 0) {
                   1506:           Print("[seq,level,q]="); Println([seq,level,q]);
                   1507:           if (level < maxLevel-1) {
                   1508:             bases = freeRes[level+1];
                   1509:             dr = reducer[level,q];
1.9       takayama 1510:             /* dr[q] = -1;  We do not need this in our reducer format. */
                   1511:             /* dr[q] should be a non-zero constant. */
1.6       takayama 1512:             newbases = SnewArrayOfFormat(bases);
                   1513:             betti_levelplus = Length(bases);
                   1514:             /*
                   1515:                bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
                   1516:             */
                   1517:             for (i=0; i<betti_levelplus; i++) {
1.9       takayama 1518:               newbases[i] = dr[q]*bases[i] - bases[i,q]*dr;
1.6       takayama 1519:             }
                   1520:             Println(["level, q =", level,q]);
                   1521:             Println("bases="); sm1_pmat(bases);
                   1522:             Println("dr="); sm1_pmat(dr);
                   1523:             Println("newbases="); sm1_pmat(newbases);
                   1524:             minRes[level+1] = newbases;
                   1525:             freeRes = minRes;
                   1526: #ifdef DEBUG
1.9       takayama 1527:             for (qq=q; qq<betti; qq++) {
1.6       takayama 1528:                 for (i=0; i<betti_levelplus; i++) {
1.9       takayama 1529:                   if ((!IsZero(newbases[i,qq])) && (redundantTable[level,qq] >0)) {
1.6       takayama 1530:                     Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
                   1531:                     Print("redundantTable ="); sm1_pmat(redundantTable[level]);
                   1532:                     Error("Stop in Sminimal for debugging.");
                   1533:                   }
                   1534:                 }
                   1535:             }
                   1536: #endif
                   1537:           }
                   1538:         }
                   1539:       }
                   1540:    }
1.10      takayama 1541:    tminRes = Stetris(minRes,redundantTable);
                   1542:    return([SpruneZeroRow(tminRes), tminRes,
1.6       takayama 1543:           [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
                   1544:   /* r[4] is the redundantTable_ordinary */
                   1545:   /* r[0] is the freeResolution */
                   1546: }
                   1547:
                   1548: /* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */
1.10      takayama 1549: /* x y (x+y-1)(x-2),  x^3-y^2, x^3 - y^2 z^2,
                   1550:    x y z (x+y+z-1) seems to be interesting, because the first syzygy
                   1551:   contains 1.
                   1552: */
                   1553:
                   1554: def CopyArray(m) {
                   1555:   local ans,i,n;
                   1556:   if (IsArray(m)) {
                   1557:      n = Length(m);
                   1558:      ans = NewArray(n);
                   1559:      for (i=0; i<n; i++) {
                   1560:        ans[i] = CopyArray(m[i]);
                   1561:      }
                   1562:      return(ans);
                   1563:   }else{
                   1564:      return(m);
                   1565:   }
                   1566: }
                   1567: HelpAdd(["CopyArray",
                   1568: ["It duplicates the argument array recursively.",
                   1569:  "Example: m=[1,[2,3]];",
                   1570:  "         a=CopyArray(m); a[1] = \"Hello\";",
                   1571:  "         Println(m); Println(a);"]]);
                   1572:
                   1573: def IsZeroVector(m) {
                   1574:   local n,i;
                   1575:   n = Length(m);
                   1576:   for (i=0; i<n; i++) {
                   1577:     if (!IsZero(m[i])) {
                   1578:       return(false);
                   1579:     }
                   1580:   }
                   1581:   return(true);
                   1582: }
                   1583:
                   1584: def SpruneZeroRow(res) {
                   1585:   local minRes, n,i,j,m, base,base2,newbase,newbase2, newMinRes;
                   1586:
                   1587:   minRes = CopyArray(res);
                   1588:   n = Length(minRes);
                   1589:   for (i=0; i<n; i++) {
                   1590:     base = minRes[i];
                   1591:     m = Length(base);
                   1592:     if (i != n-1) {
                   1593:       base2 = minRes[i+1];
                   1594:       base2 = Transpose(base2);
                   1595:     }
                   1596:     newbase = [ ];
                   1597:     newbase2 = [ ];
                   1598:     for (j=0; j<m; j++) {
                   1599:       if (!IsZeroVector(base[j])) {
                   1600:         newbase = Append(newbase,base[j]);
                   1601:         if (i != n-1) {
                   1602:           newbase2 = Append(newbase2,base2[j]);
                   1603:         }
                   1604:       }
                   1605:     }
                   1606:     minRes[i] = newbase;
                   1607:     if (i != n-1) {
                   1608:       if (newbase2 == [ ]) {
                   1609:         minRes[i+1] = [ ];
                   1610:       }else{
                   1611:         minRes[i+1] = Transpose(newbase2);
                   1612:       }
                   1613:     }
                   1614:   }
                   1615:
                   1616:   newMinRes = [ ];
                   1617:   n = Length(minRes);
                   1618:   i = 0;
                   1619:   while (i < n ) {
                   1620:     base = minRes[i];
                   1621:     if (base == [ ]) {
                   1622:       i = n; /* break; */
                   1623:     }else{
                   1624:       newMinRes = Append(newMinRes,base);
                   1625:     }
                   1626:     i++;
                   1627:   }
                   1628:   return(newMinRes);
                   1629: }
                   1630:
                   1631: def testAnnfs2(f) {
                   1632:   local a,i,n;
                   1633:   a = Sannfs2(f);
                   1634:   b=a[0];
                   1635:   n = Length(b);
                   1636:   Println("------ V-minimal free resolution -----");
                   1637:   sm1_pmat(b);
                   1638:   Println("----- Is it complex?  ---------------");
                   1639:   for (i=0; i<n-1; i++) {
                   1640:     Println(b[i+1]*b[i]);
                   1641:   }
                   1642:   return(a);
                   1643: }
                   1644: def testAnnfs3(f) {
                   1645:   local a,i,n;
                   1646:   a = Sannfs3(f);
                   1647:   b=a[0];
                   1648:   n = Length(b);
                   1649:   Println("------ V-minimal free resolution -----");
                   1650:   sm1_pmat(b);
                   1651:   Println("----- Is it complex?  ---------------");
                   1652:   for (i=0; i<n-1; i++) {
                   1653:     Println(b[i+1]*b[i]);
                   1654:   }
1.11      takayama 1655:   return(a);
                   1656: }
                   1657:
                   1658: def ToString_array(p) {
                   1659:   local ans;
                   1660:   if (IsArray(p)) {
                   1661:     ans = Map(p,"ToString_array");
                   1662:   }else{
                   1663:     ans = ToString(p);
                   1664:   }
                   1665:   return(ans);
                   1666: }
                   1667:
                   1668: /* sm1_res_div([[x],[y]],[[x^2],[x*y],[y^2]],[x,y]): */
                   1669:
                   1670: def sm1_res_div(I,J,V) {
                   1671:   I = ToString_array(I);
                   1672:   J = ToString_array(J);
                   1673:   V = ToString_array(V);
                   1674:   sm1(" [[ I J]  V ] res*div /FunctionValue set ");
                   1675: }
                   1676:
                   1677: /* It has not yet been working */
                   1678: def sm1_res_kernel_image(m,n,v) {
                   1679:   m = ToString_array(m);
                   1680:   n = ToString_array(n);
                   1681:   v = ToString_array(v);
                   1682:   sm1(" [m n v] res-kernel-image /FunctionValue set ");
                   1683: }
                   1684: def Skernel(m,v) {
                   1685:   m = ToString_array(m);
                   1686:   v = ToString_array(v);
                   1687:   sm1(" [ m v ] syz /FunctionValue set ");
                   1688: }
                   1689:
                   1690: def test3() {
                   1691:   local a1,a2,b1,b2;
                   1692:   a1 = Sannfs3("x^3-y^2*z^2");
                   1693:   a1 = a1[0];
                   1694:   a2 = Sannfs3_laScala2("x^3-y^2*z^2");
                   1695:   a2 = a2[0];
                   1696:   b1 = a1[1];
                   1697:   b2 = a2[1];
                   1698:   sm1_pmat(b2);
                   1699:   Println("  OVER ");
                   1700:   sm1_pmat(b1);
                   1701:   return([sm1_res_div(b2,b1,["x","y","z"]),b2,b1,a2,a1]);
                   1702: }
                   1703:
                   1704: def test4() {
                   1705:   local a,b;
                   1706:   a = Sannfs3_laScala2("x^3-y^2*z^2");
                   1707:   b = a[0];
                   1708:   sm1_pmat( sm1_res_kernel_image(b[0],b[1],[x,y,z]));
                   1709:   sm1_pmat( sm1_res_kernel_image(b[1],b[2],[x,y,z]));
                   1710:   return(a);
                   1711: }
                   1712:
                   1713: def sm1_gb(f,v) {
                   1714:   f =ToString_array(f);
                   1715:   v = ToString_array(v);
                   1716:   sm1(" [f v] gb /FunctionValue set ");
1.13      takayama 1717: }
                   1718:
1.11      takayama 1719:
1.12      takayama 1720: def SisComplex(a) {
                   1721:   local n,i,j,k,b,p,q;
                   1722:   n = Length(a);
                   1723:   for (i=0; i<n-1; i++) {
                   1724:     if (Length(a[i+1]) != 0) {
                   1725:       b = a[i+1]*a[i];
                   1726:       p = Length(b); q = Length(b[0]);
                   1727:       for (j=0; j<p; j++) {
                   1728:         for (k=0; k<q; k++) {
                   1729:           if (!IsZero(b[j,k])) {
                   1730:              Print("Is is not complex at ");
                   1731:              Println([i,j,k]);
                   1732:              return(false);
                   1733:           }
                   1734:         }
                   1735:       }
                   1736:     }
                   1737:   }
                   1738:   return(true);
1.14      takayama 1739: }
                   1740:
                   1741: def IsExact_h(c,v) {
                   1742:   local a;
                   1743:   v = ToString_array(v);
                   1744:   a = [c,v];
                   1745:   sm1(a," isExact_h /FunctionValue set ");
                   1746: }
                   1747: HelpAdd(["IsExact_h",
                   1748: ["IsExact_h(complex,var): bool",
                   1749:  "It checks the given complex is exact or not in D<h> (homogenized Weyl algebra)",
                   1750:  "cf. ReParse"
                   1751: ]]);
                   1752:
                   1753: def ReParse(a) {
                   1754:   local c;
                   1755:   if (IsArray(a)) {
                   1756:     c = Map(a,"ReParse");
                   1757:   }else{
                   1758:     sm1(a," toString . /c set");
                   1759:   }
                   1760:   return(c);
                   1761: }
                   1762: HelpAdd(["ReParse",
                   1763: ["Reparse(obj): obj",
                   1764:  "It parses the given object in the current ring.",
                   1765:  "Outputs from SlaScala, Sschreyer may cause a trouble in other functions,",
                   1766:  "because it uses the Schreyer order.",
                   1767:  "In this case, ReParse the outputs from these functions.",
                   1768:  "cf. IsExaxt_h"
                   1769: ]]);
1.16    ! takayama 1770:
        !          1771: def ScheckIfSchreyer(s) {
        !          1772:   local ss;
        !          1773:   sm1(" (report) (grade) switch_function /ss set ");
        !          1774:   if (ss != "module1v") {
        !          1775:      Print("ScheckIfSchreyer: from "); Println(s);
        !          1776:      Error("grade is not module1v");
        !          1777:   }
        !          1778:   /*
        !          1779:   sm1(" (report) (mmLarger) switch_function /ss set ");
        !          1780:   if (ss != "tower") {
        !          1781:      Print("ScheckIfSchreyer: from "); Println(s);
        !          1782:      Error("mmLarger is not tower");
        !          1783:   }
        !          1784:   */
        !          1785:   sm1(" [(Schreyer)] system_variable (universalNumber) dc /ss set ");
        !          1786:   if (ss != 1) {
        !          1787:      Print("ScheckIfSchreyer: from "); Println(s);
        !          1788:      Error("Schreyer order is not set.");
        !          1789:   }
        !          1790:   /* More check will be necessary. */
        !          1791:   return(true);
        !          1792: }

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