[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.7

1.7     ! takayama    1: /* $OpenXM: OpenXM/src/k097/lib/minimal/minimal.k,v 1.6 2000/05/06 07:58:37 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
                      9: #define TOTAL_STRATEGY
                     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 ");
                    135: }
                    136:
                    137: def SresolutionFrameWithTower(g,opt) {
                    138:   local gbTower, ans, ff, count, startingGB, opts, skelton,withSkel, autof,
                    139:         gbasis;
                    140:   if (Length(Arglist) >= 2) {
                    141:     if (IsInteger(opt)) count = opt;
                    142:   }else{
                    143:     count = -1;
                    144:   }
                    145:
                    146:   sm1(" setupEnvForResolution ");
                    147:   /* If I do not put this macro, homogenization
                    148:      make a strange behavior. For example,
                    149:      [(2*x*Dx + 3*y*Dy+6) (0)] homogenize returns
                    150:      [(2*x*Dx*h + 3*y*Dy*h+6*h^3) (0)].
                    151:      4/19, 2000.
                    152:   */
                    153:
                    154:   sm1(" (mmLarger) (matrix) switch_function ");
                    155:   g = Map(g,"Shomogenize");
                    156:   if (SonAutoReduce) {
                    157:     sm1("[ (AutoReduce) ] system_variable /autof set ");
                    158:     sm1("[ (AutoReduce) 1 ] system_variable ");
                    159:   }
                    160:   gbasis = Sgroebner(g);
                    161:   g = gbasis[0];
                    162:   if (SonAutoReduce) {
                    163:     sm1("[ (AutoReduce) autof] system_variable  ");
                    164:   }
                    165:
                    166:   g = Init(g);
                    167:
                    168: /*  sm1(" setupEnvForResolution-sugar "); */
                    169:   /* -sugar is fine? */
                    170:   sm1(" setupEnvForResolution ");
                    171:
                    172:   Println(g);
                    173:   startingGB = g;
                    174:   /* ans = [ SzeroMap(g) ];  It has not been implemented. see resol1.withZeroMap */
                    175:   ans = [ ];
                    176:   gbTower = [ ];
                    177:   skelton = [ ];
                    178:   while (true) {
                    179:     /* sm1(g," res0Frame /ff set "); */
                    180:     withSkel = Sres0FrameWithSkelton(g);
                    181:     ff = withSkel[0];
                    182:     ans = Append(ans, ff[0]);
                    183:     gbTower = Join([ ff[1] ], gbTower);
                    184:     skelton = Join([ withSkel[1] ], skelton);
                    185:     g = ff[0];
                    186:     if (Length(g) == 0) break;
                    187:     SsetTower( gbTower );
                    188:     if (count == 0) break;
                    189:     count = count - 1;
                    190:   }
                    191:   return([ans,Reverse(gbTower),Join([ [ ] ], Reverse(skelton)),gbasis]);
                    192: }
                    193: HelpAdd(["SresolutionFrameWithTower",
                    194: ["It returs [resolution of the initial, gbTower, skelton, gbasis]",
                    195:  "Example: Sweyl(\"x,y\");",
                    196:  "         a=SresolutionFrameWithTower([x^3,x*y,y^3-1]);"]]);
                    197:
                    198: def SresolutionFrame(f,opt) {
                    199:   local ans;
                    200:   ans = SresolutionFrameWithTower(f);
                    201:   return(ans[0]);
                    202: }
                    203: /* ---------------------------- */
                    204: def ToGradedPolySet(g) {
                    205:   sm1(g," (gradedPolySet) dc /FunctionValue set ");
                    206: }
                    207:
                    208: def NewPolynomialVector(size) {
                    209:   sm1(size," (integer) dc newPolyVector /FunctionValue set ");
                    210: }
                    211:
                    212: def  SturnOffHomogenization() {
                    213:   sm1("
                    214:     [(Homogenize)] system_variable 1 eq
                    215:     { (Warning: Homogenization and ReduceLowerTerms options are automatically turned off.) message
                    216:       [(Homogenize) 0] system_variable
                    217:       [(ReduceLowerTerms) 0] system_variable
                    218:     } {  } ifelse
                    219:   ");
                    220: }
                    221: def  SturnOnHomogenization() {
                    222:   sm1("
                    223:     [(Homogenize)] system_variable 0 eq
                    224:     { (Warning: Homogenization and ReduceLowerTerms options are automatically turned ON.) message
                    225:       [(Homogenize) 1] system_variable
                    226:       [(ReduceLowerTerms) 1] system_variable
                    227:     } {  } ifelse
                    228:   ");
                    229: }
                    230:
                    231: def SschreyerSkelton(g) {
                    232:   sm1(" [(schreyerSkelton) g] gbext /FunctionValue set ");
                    233: }
                    234: def Stoes(g) {
                    235:   if (IsArray(g)) {
                    236:     sm1(g," {toes} map /FunctionValue set ");
                    237:   }else{
                    238:     sm1(g," toes /FunctionValue set ");
                    239:   }
                    240: }
                    241: def Stoes_vec(g) {
                    242:     sm1(g," toes /FunctionValue set ");
                    243: }
                    244:
                    245: def Sres0Frame(g) {
                    246:   local ans;
                    247:   ans = Sres0FrameWithSkelton(g);
                    248:   return(ans[0]);
                    249: }
                    250: def Sres0FrameWithSkelton(g) {
                    251:   local t_syz, nexttower, m, t_gb, skel, betti,
                    252:         gg, k, i, j, pair, tmp, si, sj, grG, syzAll, gLength;
                    253:
                    254:   SturnOffHomogenization();
                    255:
                    256:   g = Stoes(g);
                    257:   skel = SschreyerSkelton(g);
                    258:   /* Print("Skelton is ");
                    259:   sm1_pmat(skel); */
                    260:   betti = Length(skel);
                    261:
                    262:   gLength = Length(g);
                    263:   grG = ToGradedPolySet(g);
                    264:   syzAll = NewPolynomialVector(betti);
                    265:   for (k=0; k<betti; k++) {
                    266:     pair = skel[k];
                    267:     i = pair[0,0];
                    268:     j = pair[0,1];
                    269:     si = pair[1,0];
                    270:     sj = pair[1,1];
                    271:     /* si g[i] + sj g[j] + \sum tmp[2][k] g[k] = 0 in res0 */
                    272:     Print(".");
                    273:
                    274:     t_syz = NewPolynomialVector(gLength);
                    275:     t_syz[i] = si;
                    276:     t_syz[j] = sj;
                    277:     syzAll[k] = t_syz;
                    278:   }
                    279:   t_syz = syzAll;
                    280:   Print("Done. betti="); Println(betti);
                    281:   /* Println(g);  g is in a format such as
                    282:     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
                    283:     [e_*x^2 , e_*x*y , 2*x*Dx*h , ...]
                    284:     [y-es*x , 3*es^4*y*Dy-es^5*x , 3*es^5*y*Dy-es^6*x , ...]
                    285:     [3*es^3*y*Dy-es^5*x ]
                    286:   */
                    287:   nexttower = Init(g);
                    288:   SturnOnHomogenization();
                    289:   return([[t_syz, nexttower],skel]);
                    290: }
                    291:
                    292:
                    293: def StotalDegree(f) {
                    294:   sm1(" [(grade) f] gbext (universalNumber) dc /FunctionValue set ");
                    295: }
                    296:
                    297: /* Sord_w(x^2*Dx*Dy,[x,-1,Dx,1]); */
                    298: def Sord_w(f,w) {
                    299:   local neww,i,n;
                    300:   n = Length(w);
                    301:   neww = NewArray(n);
                    302:   for (i=0; i<n; i=i+2) {
                    303:     neww[i] = ToString(w[i]);
                    304:   }
                    305:   for (i=1; i<n; i=i+2) {
                    306:     neww[i] = IntegerToSm1Integer(w[i]);
                    307:   }
                    308:   sm1(" f neww ord_w (universalNumber) dc /FunctionValue set ");
                    309: }
                    310:
                    311:
                    312: /* This is not satisfactory. */
                    313: def SinitOfArray(f) {
                    314:   local p,pos,top;
                    315:   if (IsArray(f)) {
                    316:      sm1(f," toes init /p set ");
                    317:      sm1(p," (es). degree (universalNumber) dc /pos set ");
                    318:      return([Init(f[pos]),pos]);
                    319:   } else {
                    320:      return(Init(f));
                    321:   }
                    322: }
                    323:
                    324: def test_SinitOfArray() {
                    325:   local f, frame,p,tower,i,j,k;
                    326:   Sweyl("x,y,z");
                    327:   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,
                    328:        -y^2*z^2 + x*z^3 + y*z^3, -z^4];
                    329:   p=SresolutionFrameWithTower(f);
                    330:   sm1_pmat(p);
                    331:   sm1_pmat(SgenerateTable(p[1]));
                    332:   return(p);
                    333:   frame = p[0];
                    334:   sm1_pmat(p[1]);
                    335:   sm1_pmat(frame);
                    336:   sm1_pmat(Map(frame[0],"SinitOfArray"));
                    337:   sm1_pmat(Map(frame[1],"SinitOfArray"));
                    338:   return(p);
                    339: }
                    340:
                    341: /* f is assumed to be a monomial with toes. */
                    342: def Sdegree(f,tower,level) {
1.6       takayama  343:   local i,ww, wd;
                    344:   /* extern WeightOfSweyl; */
                    345:   ww = WeightOfSweyl;
1.5       takayama  346:   f = Init(f);
1.1       takayama  347:   if (level <= 1) return(StotalDegree(f));
                    348:   i = Degree(f,es);
1.6       takayama  349: #ifdef TOTAL_STRATEGY
                    350:   return(StotalDegree(f)+Sdegree(tower[level-2,i],tower,level-1));
                    351: #endif
                    352:   /* Strategy must be compatible with ordering.  */
                    353:   /* Weight vector must be non-negative, too.  */
                    354:   /* See Sdegree, SgenerateTable, reductionTable. */
                    355:   wd = Sord_w(f,ww);
                    356:   return(wd+Sdegree(tower[level-2,i],tower,level-1));
                    357:
1.1       takayama  358: }
                    359:
                    360: def SgenerateTable(tower) {
                    361:   local height, n,i,j, ans, ans_at_each_floor;
                    362:   height = Length(tower);
                    363:   ans = NewArray(height);
                    364:   for (i=0; i<height; i++) {
                    365:     n = Length(tower[i]);
                    366:     ans_at_each_floor=NewArray(n);
                    367:     for (j=0; j<n; j++) {
1.6       takayama  368:       ans_at_each_floor[j] = Sdegree(tower[i,j],tower,i+1)-(i+1)
                    369:                             + OFFSET;
1.1       takayama  370:       /* Println([i,j,ans_at_each_floor[j]]); */
                    371:     }
                    372:     ans[i] = ans_at_each_floor;
                    373:   }
                    374:   return(ans);
                    375: }
                    376: Sweyl("x,y,z");
                    377: v=[[2*x*Dx + 3*y*Dy+6, 0],
                    378:    [3*x^2*Dy + 2*y*Dx, 0],
                    379:    [0,  x^2+y^2],
                    380:    [0,  x*y]];
                    381: /*  SresolutionFrameWithTower(v); */
                    382:
                    383: def SnewArrayOfFormat(p) {
                    384:   if (IsArray(p)) {
                    385:      return(Map(p,"SnewArrayOfFormat"));
                    386:   }else{
                    387:      return(null);
                    388:   }
                    389: }
1.4       takayama  390: def ScopyArray(a) {
                    391:   local n, i,ans;
                    392:   n = Length(a);
                    393:   ans = NewArray(n);
                    394:   for (i=0; i<n; i++) {
                    395:     ans[i] = a[i];
                    396:   }
                    397:   return(ans);
                    398: }
1.1       takayama  399: def SminOfStrategy(a) {
                    400:   local n,i,ans,tt;
                    401:   ans = 100000; /* very big number */
                    402:   if (IsArray(a)) {
                    403:     n = Length(a);
                    404:     for (i=0; i<n; i++) {
                    405:       if (IsArray(a[i])) {
                    406:         tt = SminOfStrategy(a[i]);
                    407:         if (tt < ans) ans = tt;
                    408:       }else{
                    409:         if (a[i] < ans) ans = a[i];
                    410:       }
                    411:     }
                    412:   }else{
                    413:      if (a < ans) ans = a;
                    414:   }
                    415:   return(ans);
                    416: }
                    417: def SmaxOfStrategy(a) {
                    418:   local n,i,ans,tt;
                    419:   ans = -100000; /* very small number */
                    420:   if (IsArray(a)) {
                    421:     n = Length(a);
                    422:     for (i=0; i<n; i++) {
                    423:       if (IsArray(a[i])) {
                    424:         tt = SmaxOfStrategy(a[i]);
                    425:         if (tt > ans) ans = tt;
                    426:       }else{
                    427:         if (a[i] > ans) ans = a[i];
                    428:       }
                    429:     }
                    430:   }else{
                    431:      if (a > ans) ans = a;
                    432:   }
                    433:   return(ans);
                    434: }
                    435:
                    436:
                    437: def SlaScala(g) {
                    438:   local rf, tower, reductionTable, skel, redundantTable, bases,
                    439:         strategy, maxOfStrategy, height, level, n, i,
                    440:         freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
1.4       takayama  441:         redundantTable_ordinary, redundant_seq_ordinary,
                    442:         reductionTable_tmp;
1.1       takayama  443:   /* extern WeightOfSweyl; */
                    444:   ww = WeightOfSweyl;
1.6       takayama  445:   Print("WeightOfSweyl="); Println(WeightOfSweyl);
1.1       takayama  446:   rf = SresolutionFrameWithTower(g);
                    447:   redundant_seq = 1;   redundant_seq_ordinary = 1;
                    448:   tower = rf[1];
                    449:   reductionTable = SgenerateTable(tower);
                    450:   skel = rf[2];
                    451:   redundantTable = SnewArrayOfFormat(rf[1]);
                    452:   redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
                    453:   reducer = SnewArrayOfFormat(rf[1]);
                    454:   freeRes = SnewArrayOfFormat(rf[1]);
                    455:   bettiTable = SsetBettiTable(rf[1],g);
                    456:
                    457:   strategy = SminOfStrategy( reductionTable );
                    458:   maxOfStrategy = SmaxOfStrategy( reductionTable );
                    459:   height = Length(reductionTable);
                    460:   while (strategy <= maxOfStrategy) {
                    461:     for (level = 0; level < height; level++) {
                    462:       n = Length(reductionTable[level]);
1.4       takayama  463:       reductionTable_tmp = ScopyArray(reductionTable[level]);
                    464:       while (SthereIs(reductionTable_tmp,strategy)) {
                    465:         i = SnextI(reductionTable_tmp,strategy,redundantTable,
                    466:                    skel,level,freeRes);
                    467:         Println([level,i]);
                    468:         reductionTable_tmp[i] = -200000;
1.1       takayama  469:         if (reductionTable[level,i] == strategy) {
                    470:            Print("Processing "); Print([level,i]);
                    471:            Print("   Strategy = "); Println(strategy);
                    472:            if (level == 0) {
                    473:              if (IsNull(redundantTable[level,i])) {
                    474:                bases = freeRes[level];
                    475:                /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
                    476:                pos = SwhereInGB(tower[0,i],rf[3,0]);
                    477:                bases[i] = rf[3,0,pos];
                    478:                redundantTable[level,i] = 0;
                    479:                redundantTable_ordinary[level,i] = 0;
                    480:                freeRes[level] = bases;
                    481:                /* Println(["GB=",i,bases,tower[0,i]]); */
                    482:              }
                    483:            }else{ /* level >= 1 */
                    484:              if (IsNull(redundantTable[level,i])) {
                    485:                bases = freeRes[level];
                    486:                f = SpairAndReduction(skel,level,i,freeRes,tower,ww);
                    487:                if (f[0] != Poly("0")) {
                    488:                   place = f[3];
                    489:                   /* (level-1, place) is the place for f[0],
                    490:                      which is a newly obtained  GB. */
                    491: #ifdef ORDINARY
                    492:                   redundantTable[level-1,place] = redundant_seq;
                    493:                   redundant_seq++;
                    494: #else
                    495:                   if (f[4] > f[5]) {
                    496:                     /* Zero in the gr-module */
                    497:                     Print("v-degree of [org,remainder] = ");
                    498:                     Println([f[4],f[5]]);
                    499:                     Print("[level,i] = "); Println([level,i]);
                    500:                     redundantTable[level-1,place] = 0;
                    501:                   }else{
                    502:                     redundantTable[level-1,place] = redundant_seq;
                    503:                     redundant_seq++;
                    504:                   }
                    505: #endif
                    506:                   redundantTable_ordinary[level-1,place]
                    507:                      =redundant_seq_ordinary;
                    508:                   redundant_seq_ordinary++;
                    509:                   bases[i] = SunitOfFormat(place,f[1])-f[1];  /* syzygy */
                    510:                   redundantTable[level,i] = 0;
                    511:                   redundantTable_ordinary[level,i] = 0;
                    512:                   /* i must be equal to f[2], I think. Double check. */
                    513:                   freeRes[level] = bases;
                    514:                   bases = freeRes[level-1];
                    515:                   bases[place] = f[0];
                    516:                   freeRes[level-1] = bases;
                    517:                   reducer[level-1,place] = f[1];
                    518:                }else{
                    519:                   redundantTable[level,i] = 0;
                    520:                   bases = freeRes[level];
                    521:                   bases[i] = f[1];  /* Put the syzygy. */
                    522:                   freeRes[level] = bases;
                    523:                }
                    524:              }
                    525:            } /* end of level >= 1 */
                    526:         }
                    527:       }
                    528:     }
                    529:     strategy++;
                    530:   }
                    531:   n = Length(freeRes);
                    532:   freeResV = SnewArrayOfFormat(freeRes);
                    533:   for (i=0; i<n; i++) {
                    534:     bases = freeRes[i];
                    535:     bases = Sbases_to_vec(bases,bettiTable[i]);
                    536:     freeResV[i] = bases;
                    537:   }
                    538:   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
                    539: }
1.4       takayama  540:
                    541: def SthereIs(reductionTable_tmp,strategy) {
                    542:   local n,i;
                    543:   n = Length(reductionTable_tmp);
                    544:   for (i=0; i<n; i++) {
                    545:     if (reductionTable_tmp[i] == strategy) {
                    546:       return(true);
                    547:     }
                    548:   }
                    549:   return(false);
                    550: }
                    551:
                    552: def SnextI(reductionTable_tmp,strategy,redundantTable,
                    553:                                   skel,level,freeRes)
                    554: {
                    555:    local ii,n,p,myindex,i,j,bases;
                    556:    n = Length(reductionTable_tmp);
                    557:    if (level == 0) {
                    558:      for (ii=0; ii<n; ii++) {
                    559:        if (reductionTable_tmp[ii] == strategy) {
                    560:           return(ii);
                    561:         }
                    562:       }
                    563:    }else{
                    564:      for (ii=0; ii<n; ii++) {
                    565:        if (reductionTable_tmp[ii] == strategy) {
                    566:          p = skel[level,ii];
                    567:          myindex = p[0];
                    568:          i = myindex[0]; j = myindex[1];
                    569:          bases = freeRes[level-1];
                    570:          if (IsNull(bases[i]) || IsNull(bases[j])) {
                    571:
                    572:          }else{
                    573:            return(ii);
                    574:          }
                    575:        }
                    576:      }
                    577:    }
1.5       takayama  578:    Print("reductionTable_tmp=");
1.4       takayama  579:    Println(reductionTable_tmp);
1.5       takayama  580:    Println("See also reductionTable, strategy, level,i");
1.4       takayama  581:    Error("SnextI: bases[i] or bases[j] is null for all combinations.");
                    582: }
                    583:
                    584:
1.1       takayama  585:
                    586: def SsetBettiTable(freeRes,g) {
                    587:   local level,i, n,bases,ans;
                    588:   ans = NewArray(Length(freeRes)+1);
                    589:   n = Length(freeRes);
                    590:   if (IsArray(g[0])) {
                    591:     ans[0] = Length(g[0]);
                    592:   }else{
                    593:     ans[0] = 1;
                    594:   }
                    595:   for (level=0; level<n; level++) {
                    596:     bases = freeRes[level];
                    597:     if (IsArray(bases)) {
                    598:       ans[level+1] = Length(bases);
                    599:     }else{
                    600:       ans[level+1] = 1;
                    601:     }
                    602:   }
                    603:   return(ans);
                    604: }
                    605:
                    606: def SwhereInGB(f,tower) {
                    607:   local i,n,p,q;
                    608:   n = Length(tower);
                    609:   for (i=0; i<n; i++) {
                    610:     p = MonomialPart(tower[i]);
                    611:     q = MonomialPart(f);
                    612:     if (p == q) return(i);
                    613:   }
                    614:   Println([f,tower]);
                    615:   Error("whereInGB : [f,myset]: f could not be found in the myset.");
                    616: }
                    617: def SunitOfFormat(pos,forms) {
                    618:   local ans,i,n;
                    619:   n = Length(forms);
                    620:   ans = NewArray(n);
                    621:   for (i=0; i<n; i++) {
                    622:     if (i != pos) {
                    623:       ans[i] = Poly("0");
                    624:     }else{
                    625:       ans[i] = Poly("1");
                    626:     }
                    627:   }
                    628:   return(ans);
                    629: }
                    630:
                    631: def Error(s) {
                    632:   sm1(" s error ");
                    633: }
                    634:
                    635: def IsNull(s) {
                    636:   if (Stag(s) == 0) return(true);
                    637:   else return(false);
                    638: }
                    639:
                    640: def StowerOf(tower,level) {
                    641:   local ans,i;
                    642:   ans = [ ];
                    643:   if (level == 0) return([[]]);
                    644:   for (i=0; i<level; i++) {
                    645:     ans = Append(ans,tower[i]);
                    646:   }
                    647:   return(Reverse(ans));
                    648: }
                    649:
                    650: def Sspolynomial(f,g) {
                    651:   if (IsArray(f)) {
                    652:     f = Stoes_vec(f);
                    653:   }
                    654:   if (IsArray(g)) {
                    655:     g = Stoes_vec(g);
                    656:   }
                    657:   sm1("f g spol /FunctionValue set");
                    658: }
                    659:
                    660: def MonomialPart(f) {
                    661:   sm1(" [(lmonom) f] gbext /FunctionValue set ");
                    662: }
                    663:
                    664: def SwhereInTower(f,tower) {
                    665:   local i,n,p,q;
                    666:   if (f == Poly("0")) return(-1);
                    667:   n = Length(tower);
                    668:   for (i=0; i<n; i++) {
                    669:     p = MonomialPart(tower[i]);
                    670:     q = MonomialPart(f);
                    671:     if (p == q) return(i);
                    672:   }
                    673:   Println([f,tower]);
                    674:   Error("[f,tower]: f could not be found in the tower.");
                    675: }
                    676:
                    677: def Stag(f) {
                    678:   sm1(f," tag (universalNumber) dc /FunctionValue set");
                    679: }
                    680:
                    681: def SpairAndReduction(skel,level,ii,freeRes,tower,ww) {
                    682:   local i, j, myindex, p, bases, tower2, gi, gj,
                    683:        si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
                    684:        vdeg,vdeg_reduced;
                    685:   Println("SpairAndReduction:");
                    686:
                    687:   if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
                    688:   p = skel[level,ii];
                    689:   myindex = p[0];
                    690:   i = myindex[0]; j = myindex[1];
                    691:   bases = freeRes[level-1];
                    692:   Println(["p and bases ",p,bases]);
                    693:   if (IsNull(bases[i]) || IsNull(bases[j])) {
                    694:     Println([level,i,j,bases[i],bases[j]]);
                    695:     Error("level, i, j : bases[i], bases[j]  must not be NULL.");
                    696:   }
                    697:
                    698:   tower2 = StowerOf(tower,level-1);
                    699:   SsetTower(tower2);
                    700:   /** sm1(" show_ring ");   */
                    701:
                    702:   gi = Stoes_vec(bases[i]);
                    703:   gj = Stoes_vec(bases[j]);
                    704:
                    705:   ssp = Sspolynomial(gi,gj);
                    706:   si = ssp[0,0];
                    707:   sj = ssp[0,1];
                    708:   syzHead = si*es^i;
                    709:   /* This will be the head term, I think. But, double check. */
                    710:   Println([si*es^i,sj*es^j]);
                    711:
                    712:   Print("[gi, gj] = "); Println([gi,gj]);
                    713:   sm1(" [(Homogenize)] system_variable message ");
                    714:   Print("Reduce the element "); Println(si*gi+sj*gj);
                    715:   Print("by  "); Println(bases);
                    716:
                    717:   tmp = Sreduction(si*gi+sj*gj, bases);
                    718:
                    719:   Print("result is "); Println(tmp);
                    720:
1.3       takayama  721:   /* This is essential part for V-minimal resolution. */
                    722:   /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
                    723:   vdeg = SvDegree(si*gi,tower,level-1,ww);
1.1       takayama  724:   vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
                    725:   Print("vdegree of the original = "); Println(vdeg);
                    726:   Print("vdegree of the remainder = "); Println(vdeg_reduced);
                    727:
                    728:   t_syz = tmp[2];
                    729:   si = si*tmp[1]+t_syz[i];
                    730:   sj = sj*tmp[1]+t_syz[j];
                    731:   t_syz[i] = si;
                    732:   t_syz[j] = sj;
                    733:   pos = SwhereInTower(syzHead,tower[level]);
                    734:   pos2 = SwhereInTower(tmp[0],tower[level-1]);
                    735:   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced];
                    736:   /* pos is the place to put syzygy at level. */
                    737:   /* pos2 is the place to put a new GB at level-1. */
                    738:   Println(ans);
                    739:   return(ans);
                    740: }
                    741:
                    742: def Sreduction(f,myset) {
                    743:   local n, indexTable, set2, i, j, tmp, t_syz;
                    744:   n = Length(myset);
                    745:   indexTable = NewArray(n);
                    746:   set2 = [ ];
                    747:   j = 0;
                    748:   for (i=0; i<n; i++) {
                    749:     if (IsNull(myset[i])) {
                    750:       indexTable[i] = -1;
                    751: /*    }else if (myset[i] == Poly("0")) {
                    752:       indexTable[i] = -1;  */
                    753:     }else{
                    754:       set2 = Append(set2,Stoes_vec(myset[i]));
                    755:       indexTable[i] = j;
                    756:       j++;
                    757:     }
                    758:   }
                    759:   sm1(" f toes set2 (gradedPolySet) dc reduction /tmp set ");
                    760:   t_syz = NewArray(n);
                    761:   for (i=0; i<n; i++) {
                    762:     if (indexTable[i] != -1) {
                    763:       t_syz[i] = tmp[2, indexTable[i]];
                    764:     }else{
                    765:       t_syz[i] = Poly("0");
                    766:     }
                    767:   }
                    768:   return([tmp[0],tmp[1],t_syz]);
                    769: }
                    770:
                    771: def Warning(s) {
                    772:   Print("Warning: ");
                    773:   Println(s);
                    774: }
                    775: def RingOf(f) {
                    776:   local r;
                    777:   if (IsPolynomial(f)) {
                    778:     if (f != Poly("0")) {
                    779:       sm1(f," (ring) dc /r set ");
                    780:     }else{
                    781:       sm1(" [(CurrentRingp)] system_variable /r set ");
                    782:     }
                    783:   }else{
                    784:     Warning("RingOf(f): the argument f must be a polynomial. Return the current ring.");
                    785:     sm1(" [(CurrentRingp)] system_variable /r set ");
                    786:   }
                    787:   return(r);
                    788: }
                    789:
                    790: def Sfrom_es(f,size) {
                    791:   local c,ans, i, d, myes, myee, j,n,r,ans2;
                    792:   if (Length(Arglist) < 2) size = -1;
                    793:   if (IsArray(f)) return(f);
                    794:   r = RingOf(f);
                    795:   myes = PolyR("es",r);
                    796:   myee = PolyR("e_",r);
                    797:   if (Degree(f,myee) > 0 && size == -1) {
                    798:     if (size == -1) {
                    799:        sm1(f," (array) dc /ans set");
                    800:        return(ans);
                    801:     }
                    802:   }
                    803:
                    804: /*
                    805:     Coefficients(x^2-1,x):
                    806:     [    [    2 , 0 ]  , [    1 , -1 ]  ]
                    807: */
                    808:   if (Degree(f,myee) > 0) {
                    809:     c = Coefficients(f,myee);
                    810:   }else{
                    811:     c = Coefficients(f,myes);
                    812:   }
                    813:   if (size < 0) {
                    814:     size = c[0,0]+1;
                    815:   }
                    816:   ans = NewArray(size);
                    817:   for (i=0; i<size; i++) {ans[i] = 0;}
                    818:   n = Length(c[0]);
                    819:   for (j=0; j<n; j++) {
                    820:     d = c[0,j];
                    821:     ans[d] = c[1,j];
                    822:   }
                    823:   return(ans);
                    824: }
                    825:
                    826: def Sbases_to_vec(bases,size) {
                    827:   local n, giveSize, newbases,i;
                    828:   /*  bases = [1+es*x, [1,2,3*x]] */
                    829:   if (Length(Arglist) > 1) {
                    830:     giveSize = true;
                    831:   }else{
                    832:     giveSize = false;
                    833:   }
                    834:   n = Length(bases);
                    835:   newbases = NewArray(n);
                    836:   for (i=0; i<n; i++) {
                    837:      if (giveSize) {
                    838:        newbases[i] = Sfrom_es(bases[i], size);
                    839:      }else{
                    840:        newbases[i] = Sfrom_es(bases[i]);
                    841:      }
                    842:   }
                    843:   return(newbases);
                    844: }
                    845:
                    846: def Sminimal(g) {
                    847:   local r, freeRes, redundantTable, reducer, maxLevel,
                    848:         minRes, seq, maxSeq, level, betti, q, bases, dr,
                    849:         betti_levelplus, newbases, i, j,qq;
                    850:   r = SlaScala(g);
                    851:   /* Should I turn off the tower?? */
                    852:   freeRes = r[0];
                    853:   redundantTable = r[1];
                    854:   reducer = r[2];
                    855:   minRes = SnewArrayOfFormat(freeRes);
                    856:   seq = 0;
                    857:   maxSeq = SgetMaxSeq(redundantTable);
                    858:   maxLevel = Length(freeRes);
                    859:   for (level = 0; level < maxLevel; level++) {
                    860:     minRes[level] = freeRes[level];
                    861:   }
                    862:   seq=maxSeq+1;
                    863:   while (seq > 1) {
                    864:     seq--;
                    865:     for (level = 0; level < maxLevel; level++) {
                    866:       betti = Length(freeRes[level]);
                    867:       for (q = 0; q<betti; q++) {
                    868:         if (redundantTable[level,q] == seq) {
                    869:           Print("[seq,level,q]="); Println([seq,level,q]);
                    870:           if (level < maxLevel-1) {
                    871:             bases = freeRes[level+1];
                    872:             dr = reducer[level,q];
                    873:             dr[q] = -1;
                    874:             newbases = SnewArrayOfFormat(bases);
                    875:             betti_levelplus = Length(bases);
                    876:             /*
                    877:                bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
                    878:             */
                    879:             for (i=0; i<betti_levelplus; i++) {
                    880:               newbases[i] = bases[i] + bases[i,q]*dr;
                    881:             }
                    882:             Println(["level, q =", level,q]);
                    883:             Println("bases="); sm1_pmat(bases);
                    884:             Println("dr="); sm1_pmat(dr);
                    885:             Println("newbases="); sm1_pmat(newbases);
                    886:             minRes[level+1] = newbases;
                    887:             freeRes = minRes;
                    888: #ifdef DEBUG
                    889:             for (qq=0; qq<betti; qq++) {
                    890:               if ((redundantTable[level,qq] >= seq) &&
                    891:                   (redundantTable[level,qq] <= maxSeq)) {
                    892:                 for (i=0; i<betti_levelplus; i++) {
                    893:                   if (!IsZero(newbases[i,qq])) {
                    894:                     Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
                    895:                     Print("redundantTable ="); sm1_pmat(redundantTable[level]);
                    896:                     Error("Stop in Sminimal for debugging.");
                    897:                   }
                    898:                 }
                    899:               }
                    900:             }
                    901: #endif
                    902:           }
                    903:         }
                    904:       }
                    905:     }
                    906:    }
                    907:    return([Stetris(minRes,redundantTable),
1.3       takayama  908:           [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
1.1       takayama  909:   /* r[4] is the redundantTable_ordinary */
1.3       takayama  910:   /* r[0] is the freeResolution */
1.1       takayama  911: }
                    912:
                    913:
                    914: def IsZero(f) {
                    915:   if (IsPolynomial(f)) {
                    916:     return( f == Poly("0"));
                    917:   }else if (IsInteger(f)) {
                    918:     return( f == 0);
                    919:   }else if (IsSm1Integer(f)) {
                    920:     return( f == true );
                    921:   }else if (IsDouble(f)) {
                    922:     return( f == 0.0 );
                    923:   }else if (IsRational(f)) {
                    924:     return(IsZero(Denominator(f)));
                    925:   }else{
                    926:     Error("IsZero: cannot deal with this data type.");
                    927:   }
                    928: }
                    929: def SgetMaxSeq(redundantTable) {
                    930:    local level,i,n,ans, levelMax,bases;
                    931:    levelMax = Length( redundantTable );
                    932:    ans = 0;
                    933:    for (level = 0; level < levelMax; level++) {
                    934:      bases = redundantTable[level];
                    935:      n = Length(bases);
                    936:      for (i=0; i<n; i++) {
                    937:        if (IsInteger( bases[i] )) {
                    938:           if (bases[i] > ans) {
                    939:              ans = bases[i];
                    940:           }
                    941:        }
                    942:      }
                    943:    }
                    944:    return(ans);
                    945: }
                    946:
                    947: def Stetris(freeRes,redundantTable) {
                    948:   local level, i, j, resLength, minRes,
                    949:         bases, newbases, newbases2;
                    950:   minRes = SnewArrayOfFormat(freeRes);
                    951:   resLength = Length( freeRes );
                    952:   for (level=0; level<resLength; level++) {
                    953:     bases = freeRes[level];
                    954:     newbases = SnewArrayOfFormat(bases);
                    955:     betti = Length(bases); j = 0;
                    956:     /* Delete rows */
                    957:     for (i=0; i<betti; i++) {
                    958:       if (redundantTable[level,i] < 1) {
                    959:          newbases[j] = bases[i];
                    960:          j++;
                    961:       }
                    962:     }
                    963:     bases = SfirstN(newbases,j);
                    964:     if (level > 0) {
                    965:       /* Delete columns */
                    966:       newbases = Transpose(bases);
                    967:       betti = Length(newbases); j = 0;
                    968:       newbases2 = SnewArrayOfFormat(newbases);
                    969:       for (i=0; i<betti; i++) {
                    970:         if (redundantTable[level-1,i] < 1) {
                    971:            newbases2[j] = newbases[i];
                    972:            j++;
                    973:         }
                    974:       }
                    975:       newbases = Transpose(SfirstN(newbases2,j));
                    976:     }else{
                    977:       newbases = bases;
                    978:     }
                    979:     Println(["level=", level]);
                    980:     sm1_pmat(bases);
                    981:     sm1_pmat(newbases);
                    982:
                    983:     minRes[level] = newbases;
                    984:   }
                    985:   return(minRes);
                    986: }
                    987:
                    988: def SfirstN(bases,k) {
                    989:    local ans,i;
                    990:    ans = NewArray(k);
                    991:    for (i=0; i<k; i++) {
                    992:      ans[i] = bases[i];
                    993:    }
                    994:    return(ans);
                    995: }
                    996:
                    997:
                    998: /* usage:  tt is tower. ww is weight.
                    999:     a = SresolutionFrameWithTower(v);
                   1000:     tt = a[1];
                   1001:     ww = [x,1,y,1,Dx,1,Dy,1];
                   1002:     SvDegree(x*es,tt,1,ww):
                   1003:
                   1004: In(17)=tt:
                   1005: [[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 ]  ,
                   1006:  [es*y , 3*es^3*y*Dy , 3*es^5*y*Dy , 3*x*Dy , es^2*y^2 , 9*y*Dy^2 ]  ,
                   1007:  [3*es^3*y*Dy ]  ]
                   1008: In(18)=SvDegree(x*es,tt,1,ww):
                   1009: 3
                   1010: In(19)=SvDegree(x*es^3,tt,1,ww):
                   1011: 4
                   1012: In(20)=SvDegree(x,tt,2,ww):
                   1013: 4
                   1014:
                   1015: */
                   1016: def SvDegree(f,tower,level,w) {
                   1017:   local i,ans;
                   1018:   if (IsZero(f)) return(null);
1.3       takayama 1019:   f = Init(f);
1.1       takayama 1020:   if (level <= 0) {
                   1021:     return(Sord_w(f,w));
                   1022:   }
                   1023:   i = Degree(f,es);
                   1024:   ans = Sord_w(f,w) +
                   1025:         SvDegree(tower[level-1,i],tower,level-1,w);
                   1026:   return(ans);
                   1027: }
                   1028:
1.2       takayama 1029: def Sannfs(f,v) {
                   1030:   local f2;
                   1031:   f2 = ToString(f);
                   1032:   if (IsArray(v)) {
                   1033:      v = Map(v,"ToString");
                   1034:   }
                   1035:   sm1(" [f2 v] annfs /FunctionValue set ");
                   1036: }
                   1037:
                   1038: /* Sannfs2("x^3-y^2"); */
                   1039: def Sannfs2(f) {
                   1040:   local p,pp;
                   1041:   p = Sannfs(f,"x,y");
1.6       takayama 1042:   sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
1.5       takayama 1043: /*
                   1044:   Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1],
                   1045:                ["x",-1,"y",-1,"Dx",1,"Dy",1]]); */
1.6       takayama 1046:   /* Sweyl("x,y",[["x",1,"y",1,"Dx",1,"Dy",1,"h",1]]); */
                   1047:   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
                   1048:   pp = Map(p,"Spoly");
                   1049:   return(Sminimal_v(pp));
                   1050:   /* return(Sminimal(pp)); */
                   1051: }
                   1052:
                   1053: /* Do not forget to turn on TOTAL_STRATEGY */
                   1054: def Sannfs2_laScala(f) {
                   1055:   local p,pp;
                   1056:   p = Sannfs(f,"x,y");
                   1057:   /*   Do not make laplace transform.
                   1058:     sm1(" p 0 get { [(x) (y) (Dx) (Dy)] laplace0 } map /p set ");
                   1059:     p = [p];
                   1060:   */
                   1061:   Sweyl("x,y",[["x",-1,"y",-1,"Dx",1,"Dy",1]]);
1.2       takayama 1062:   pp = Map(p[0],"Spoly");
                   1063:   return(Sminimal(pp));
                   1064: }
                   1065:
1.3       takayama 1066: def Sannfs3(f) {
                   1067:   local p,pp;
                   1068:   p = Sannfs(f,"x,y,z");
1.6       takayama 1069:   sm1(" p 0 get { [(x) (y) (z) (Dx) (Dy) (Dz)] laplace0 } map /p set ");
1.3       takayama 1070:   Sweyl("x,y,z",[["x",-1,"y",-1,"z",-1,"Dx",1,"Dy",1,"Dz",1]]);
1.6       takayama 1071:   pp = Map(p,"Spoly");
                   1072:   return(Sminimal_v(pp));
1.3       takayama 1073: }
                   1074:
1.2       takayama 1075: /*
                   1076:   The betti numbers of most examples are 2,1. (0-th and 1-th).
                   1077:   a=Sannfs2("x*y*(x+y-1)"); ==> The betti numbers are 3, 2.
                   1078:   a=Sannfs2("x^3-y^2-x");    : it causes an error. It should be fixed.
1.3       takayama 1079:   a=Sannfs2("x*y*(x-y)");    : it causes an error. It should be fixed.
1.2       takayama 1080:
                   1081: */
                   1082:
1.5       takayama 1083:
                   1084:
1.6       takayama 1085: /*  The below does not use LaScala-Stillman's algorithm. */
1.5       takayama 1086: def Sschreyer(g) {
                   1087:   local rf, tower, reductionTable, skel, redundantTable, bases,
                   1088:         strategy, maxOfStrategy, height, level, n, i,
                   1089:         freeRes,place, f, reducer,pos, redundant_seq,bettiTable,freeResV,ww,
                   1090:         redundantTable_ordinary, redundant_seq_ordinary,
1.6       takayama 1091:         reductionTable_tmp,c2,ii,nn, m,ii, jj, reducerBase;
1.5       takayama 1092:   /* extern WeightOfSweyl; */
                   1093:   ww = WeightOfSweyl;
                   1094:   Print("WeghtOfSweyl="); Println(WeightOfSweyl);
                   1095:   rf = SresolutionFrameWithTower(g);
                   1096:   redundant_seq = 1;   redundant_seq_ordinary = 1;
                   1097:   tower = rf[1];
                   1098:   reductionTable = SgenerateTable(tower);
                   1099:   skel = rf[2];
                   1100:   redundantTable = SnewArrayOfFormat(rf[1]);
                   1101:   redundantTable_ordinary = SnewArrayOfFormat(rf[1]);
                   1102:   reducer = SnewArrayOfFormat(rf[1]);
                   1103:   freeRes = SnewArrayOfFormat(rf[1]);
                   1104:   bettiTable = SsetBettiTable(rf[1],g);
                   1105:
                   1106:   height = Length(reductionTable);
                   1107:   for (level = 0; level < height; level++) {
                   1108:       n = Length(reductionTable[level]);
                   1109:       for (i=0; i<n; i++) {
                   1110:            Println([level,i]);
                   1111:            Print("Processing "); Print([level,i]);
                   1112:            if (level == 0) {
                   1113:              if (IsNull(redundantTable[level,i])) {
                   1114:                bases = freeRes[level];
                   1115:                /* Println(["At floor : GB=",i,bases,tower[0,i]]); */
                   1116:                pos = SwhereInGB(tower[0,i],rf[3,0]);
                   1117:                bases[i] = rf[3,0,pos];
                   1118:                /* redundantTable[level,i] = 0;
                   1119:                redundantTable_ordinary[level,i] = 0; */
                   1120:                freeRes[level] = bases;
                   1121:                /* Println(["GB=",i,bases,tower[0,i]]); */
                   1122:              }
                   1123:            }else{ /* level >= 1 */
                   1124:              if (IsNull(redundantTable[level,i])) {
                   1125:                bases = freeRes[level];
                   1126:                f = SpairAndReduction2(skel,level,i,freeRes,tower,
                   1127:                                       ww,redundantTable);
                   1128:                if (f[0] != Poly("0")) {
                   1129:                   place = f[3];
                   1130:                   /* (level-1, place) is the place for f[0],
                   1131:                      which is a newly obtained  GB. */
                   1132: #ifdef ORDINARY
                   1133:                   redundantTable[level-1,place] = redundant_seq;
                   1134:                   redundant_seq++;
                   1135: #else
                   1136:                   if (f[4] > f[5]) {
                   1137:                     /* Zero in the gr-module */
                   1138:                     Print("v-degree of [org,remainder] = ");
                   1139:                     Println([f[4],f[5]]);
                   1140:                     Print("[level,i] = "); Println([level,i]);
                   1141:                     redundantTable[level-1,place] = 0;
                   1142:                   }else{
                   1143:                     redundantTable[level-1,place] = redundant_seq;
                   1144:                     redundant_seq++;
                   1145:                   }
                   1146: #endif
                   1147:                   redundantTable_ordinary[level-1,place]
                   1148:                      =redundant_seq_ordinary;
                   1149:                   redundant_seq_ordinary++;
                   1150:                   bases[i] = SunitOfFormat(place,f[1])-f[1];  /* syzygy */
                   1151:                   /* redundantTable[level,i] = 0;
                   1152:                   redundantTable_ordinary[level,i] = 0; */
                   1153:                   /* i must be equal to f[2], I think. Double check. */
                   1154:
                   1155:                   /* Correction Of Constant */
1.7     ! takayama 1156:                   c2 = f[6];  /* or -f[6]?  Double check. */
        !          1157:                   Print("c2="); Println(c2);
1.5       takayama 1158:                   nn = Length(bases);
                   1159:                   for (ii=0; ii<nn;ii++) {
1.6       takayama 1160:                      if ((ii != place) && (! IsNull(bases[ii]))) {
1.7     ! takayama 1161:                        m = Length(bases[ii]);
        !          1162:                        for (jj=0; jj<m; jj++) {
        !          1163:                          if (jj != place) {
        !          1164:                            bases[ii,jj] = bases[ii,jj]*c2;
        !          1165:                          }
        !          1166:                        }
1.5       takayama 1167:                      }
                   1168:                   }
                   1169:
1.7     ! takayama 1170:                   Print("Old freeRes[level] = "); sm1_pmat(freeRes[level]);
1.5       takayama 1171:                   freeRes[level] = bases;
1.7     ! takayama 1172:                   Print("New freeRes[level] = "); sm1_pmat(freeRes[level]);
1.6       takayama 1173:
                   1174:                  /* Update the freeRes[level-1] */
1.7     ! takayama 1175:                   Print("Old freeRes[level-1] = "); sm1_pmat(freeRes[level-1]);
1.6       takayama 1176:                   bases = freeRes[level-1];
                   1177:                   bases[place] = f[0];
                   1178:                   freeRes[level-1] = bases;
1.7     ! takayama 1179:                   Print("New freeRes[level-1] = "); sm1_pmat(freeRes[level-1]);
1.6       takayama 1180:
1.5       takayama 1181:                   reducer[level-1,place] = f[1];
                   1182:                }else{
                   1183:                   /* redundantTable[level,i] = 0; */
                   1184:                   bases = freeRes[level];
                   1185:                   bases[i] = f[1];  /* Put the syzygy. */
                   1186:                   freeRes[level] = bases;
                   1187:                }
                   1188:              }  /* end of level >= 1 */
                   1189:           }
                   1190:     } /* i loop */
1.6       takayama 1191:
                   1192:     /* Triangulate reducer */
                   1193:     if (level >= 1) {
                   1194:       Println(" ");
                   1195:       Print("Triangulating reducer at level "); Println(level-1);
                   1196:       reducerBase = reducer[level-1];
                   1197:       Print("reducerBase=");  Println(reducerBase);
                   1198:       m = Length(reducerBase);
                   1199:       for (ii=m-1; ii>=0; ii--) {
                   1200:         if (!IsNull(reducerBase[ii])) {
                   1201:            for (jj=ii-1; jj>=0; jj--) {
                   1202:              if (!IsNull(reducerBase[jj])) {
                   1203:               if (!IsZero(reducerBase[jj,ii])) {
                   1204:                 reducerBase[jj] = reducerBase[jj]-reducerBase[jj,ii]*reducerBase[ii];
                   1205:               }
                   1206:              }
                   1207:            }
                   1208:          }
                   1209:        }
                   1210:        Println("New reducer");
                   1211:        sm1_pmat(reducerBase);
                   1212:        reducer[level-1] = reducerBase;
                   1213:     }
                   1214:
1.5       takayama 1215:   } /* level loop */
                   1216:   n = Length(freeRes);
                   1217:   freeResV = SnewArrayOfFormat(freeRes);
                   1218:   for (i=0; i<n; i++) {
                   1219:     bases = freeRes[i];
                   1220:     bases = Sbases_to_vec(bases,bettiTable[i]);
                   1221:     freeResV[i] = bases;
                   1222:   }
1.6       takayama 1223:
                   1224:   /* Mark the non-redundant elements. */
                   1225:   for (i=0; i<n; i++) {
                   1226:     m = Length(redundantTable[i]);
                   1227:     for (jj=0; jj<m; jj++) {
                   1228:       if (IsNull(redundantTable[i,jj])) {
                   1229:         redundantTable[i,jj] = 0;
                   1230:       }
                   1231:     }
                   1232:   }
                   1233:
                   1234:
1.5       takayama 1235:   return([freeResV, redundantTable,reducer,bettiTable,redundantTable_ordinary]);
                   1236: }
                   1237:
                   1238: def SpairAndReduction2(skel,level,ii,freeRes,tower,ww,redundantTable) {
                   1239:   local i, j, myindex, p, bases, tower2, gi, gj,
                   1240:        si, sj, tmp, t_syz, pos, ans, ssp, syzHead,pos2,
                   1241:        vdeg,vdeg_reduced,n,c2;
1.6       takayama 1242:   Println("SpairAndReduction2 : -------------------------");
1.5       takayama 1243:
                   1244:   if (level < 1) Error("level should be >= 1 in SpairAndReduction.");
                   1245:   p = skel[level,ii];
                   1246:   myindex = p[0];
                   1247:   i = myindex[0]; j = myindex[1];
                   1248:   bases = freeRes[level-1];
                   1249:   Println(["p and bases ",p,bases]);
                   1250:   if (IsNull(bases[i]) || IsNull(bases[j])) {
                   1251:     Println([level,i,j,bases[i],bases[j]]);
                   1252:     Error("level, i, j : bases[i], bases[j]  must not be NULL.");
                   1253:   }
                   1254:
                   1255:   tower2 = StowerOf(tower,level-1);
                   1256:   SsetTower(tower2);
                   1257:   /** sm1(" show_ring ");   */
                   1258:
                   1259:   gi = Stoes_vec(bases[i]);
                   1260:   gj = Stoes_vec(bases[j]);
                   1261:
                   1262:   ssp = Sspolynomial(gi,gj);
                   1263:   si = ssp[0,0];
                   1264:   sj = ssp[0,1];
                   1265:   syzHead = si*es^i;
                   1266:   /* This will be the head term, I think. But, double check. */
                   1267:   Println([si*es^i,sj*es^j]);
                   1268:
                   1269:   Print("[gi, gj] = "); Println([gi,gj]);
                   1270:   sm1(" [(Homogenize)] system_variable message ");
                   1271:   Print("Reduce the element "); Println(si*gi+sj*gj);
                   1272:   Print("by  "); Println(bases);
                   1273:
                   1274:   tmp = Sreduction(si*gi+sj*gj, bases);
                   1275:
                   1276:   Print("result is "); Println(tmp);
1.6       takayama 1277:   if (!IsZero(tmp[0])) {
                   1278:     Print("Error: base = ");
                   1279:     Println(Map(bases,"Stoes_vec"));
                   1280:     Error("SpairAndReduction2: the remainder should be zero. See tmp. tower2. show_ring.");
                   1281:   }
1.5       takayama 1282:   t_syz = tmp[2];
                   1283:   si = si*tmp[1]+t_syz[i];
                   1284:   sj = sj*tmp[1]+t_syz[j];
                   1285:   t_syz[i] = si;
                   1286:   t_syz[j] = sj;
                   1287:
                   1288:   c2 = null;
                   1289:   /* tmp[0] must be zero */
                   1290:   n = Length(t_syz);
                   1291:   for (i=0; i<n; i++) {
1.6       takayama 1292:      if (IsConstant(t_syz[i])){
                   1293:       if (!IsZero(t_syz[i])) {
1.5       takayama 1294:        if (IsNull(redundantTable[level-1,i])) {
                   1295:          /* i must equal to pos2 below. */
                   1296:          c2 = -t_syz[i];
1.6       takayama 1297:          tmp[0] = c2*Stoes_vec(freeRes[level-1,i]);
1.5       takayama 1298:          t_syz[i] = 0;
1.6       takayama 1299:          /* tmp[0] = t_syz . g */
1.5       takayama 1300:          /* break; does not work. Use */
                   1301:          i = n;
                   1302:        }
1.6       takayama 1303:       }
1.5       takayama 1304:      }
                   1305:   }
                   1306:
                   1307:   /* This is essential part for V-minimal resolution. */
                   1308:   /* vdeg = SvDegree(si*gi+sj*gj,tower,level-1,ww); */
                   1309:   vdeg = SvDegree(si*gi,tower,level-1,ww);
                   1310:   vdeg_reduced = SvDegree(tmp[0],tower,level-1,ww);
                   1311:   Print("vdegree of the original = "); Println(vdeg);
                   1312:   Print("vdegree of the remainder = "); Println(vdeg_reduced);
                   1313:
                   1314:   pos = SwhereInTower(syzHead,tower[level]);
                   1315:   pos2 = SwhereInTower(tmp[0],tower[level-1]);
                   1316:   ans = [tmp[0],t_syz,pos,pos2,vdeg,vdeg_reduced,c2];
                   1317:   /* pos is the place to put syzygy at level. */
                   1318:   /* pos2 is the place to put a new GB at level-1. */
                   1319:   Println(ans);
1.6       takayama 1320:   Println("  ");
1.5       takayama 1321:   return(ans);
                   1322: }
1.6       takayama 1323:
                   1324: def Sminimal_v(g) {
                   1325:   local r, freeRes, redundantTable, reducer, maxLevel,
                   1326:         minRes, seq, maxSeq, level, betti, q, bases, dr,
                   1327:         betti_levelplus, newbases, i, j,qq;
                   1328:   r = Sschreyer(g);
                   1329:   sm1_pmat(r);
                   1330:   Debug_Sminimal_v = r;
                   1331:   Println(" Return value of Schreyer(g) is set to Debug_Sminimal_v");
                   1332:   /* Should I turn off the tower?? */
                   1333:   freeRes = r[0];
                   1334:   redundantTable = r[1];
                   1335:   reducer = r[2];
                   1336:   minRes = SnewArrayOfFormat(freeRes);
                   1337:   seq = 0;
                   1338:   maxSeq = SgetMaxSeq(redundantTable);
                   1339:   maxLevel = Length(freeRes);
                   1340:   for (level = 0; level < maxLevel; level++) {
                   1341:     minRes[level] = freeRes[level];
                   1342:   }
                   1343:   for (level = 0; level < maxLevel; level++) {
                   1344:       betti = Length(freeRes[level]);
                   1345:       for (q = betti-1; q>=0; q--) {
                   1346:         if (redundantTable[level,q] > 0) {
                   1347:           Print("[seq,level,q]="); Println([seq,level,q]);
                   1348:           if (level < maxLevel-1) {
                   1349:             bases = freeRes[level+1];
                   1350:             dr = reducer[level,q];
                   1351:             dr[q] = -1;
                   1352:             newbases = SnewArrayOfFormat(bases);
                   1353:             betti_levelplus = Length(bases);
                   1354:             /*
                   1355:                bases[i,j] ---> bases[i,j]+bases[i,q]*dr[j]
                   1356:             */
                   1357:             for (i=0; i<betti_levelplus; i++) {
                   1358:               newbases[i] = bases[i] + bases[i,q]*dr;
                   1359:             }
                   1360:             Println(["level, q =", level,q]);
                   1361:             Println("bases="); sm1_pmat(bases);
                   1362:             Println("dr="); sm1_pmat(dr);
                   1363:             Println("newbases="); sm1_pmat(newbases);
                   1364:             minRes[level+1] = newbases;
                   1365:             freeRes = minRes;
                   1366: #ifdef DEBUG
                   1367: /*  Do it later.
                   1368:             for (qq=0; qq<betti; qq++) {
                   1369:                 for (i=0; i<betti_levelplus; i++) {
                   1370:                   if (!IsZero(newbases[i,qq])) {
                   1371:                     Println(["[i,qq]=",[i,qq]," is not zero in newbases."]);
                   1372:                     Print("redundantTable ="); sm1_pmat(redundantTable[level]);
                   1373:                     Error("Stop in Sminimal for debugging.");
                   1374:                   }
                   1375:                 }
                   1376:             }
                   1377: */
                   1378: #endif
                   1379:           }
                   1380:         }
                   1381:       }
                   1382:    }
                   1383:    return([Stetris(minRes,redundantTable),
                   1384:           [ minRes, redundantTable, reducer,r[3],r[4]],r[0]]);
                   1385:   /* r[4] is the redundantTable_ordinary */
                   1386:   /* r[0] is the freeResolution */
                   1387: }
                   1388:
                   1389: /* Sannfs2("x*y*(x-y)*(x+y)"); is a test problem */

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