[BACK]Return to module_syz.rr CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / testing / noro

Annotation of OpenXM/src/asir-contrib/testing/noro/module_syz.rr, Revision 1.8

1.1       noro        1: module newsyz;
                      2:
1.4       noro        3: localf module_syz, module_syz_old;
1.3       noro        4: localf simplify_syz, icont, mod, remove_cont,ordcheck;
                      5: localf complsb, complsb_sd, sortlsb, find_pos, find_pos, reduce, lres_setup, dpm_sort1, comp_pos;
                      6: localf fres,minres,sres,minsres,lres, create_base_ord, simplify_k, simplify_by_k, remove_k, remove_k1, extract_nonzero;
                      7: localf nonzero, phi, syz_check, renumber_pos, compress, compress_h;
1.8     ! noro        8: localf syz_check0,phi0,todpmlist,dpmlisttollist,comp_lex;
1.1       noro        9:
                     10: /* F : a list of (lists or polynomials),
                     11:    V : a variable list, H >1=> over GF(H), H=0,1=> over Q
                     12:    O : term order
                     13:    return: [GS,G]
                     14:    GS : a GB of syz(F) wrt [1,O] (POT), G: a GB of F wrt [1,O]
                     15: */
                     16:
1.3       noro       17: // return [dpmlist F,rank N]
                     18: def todpmlist(F,V)
1.1       noro       19: {
1.3       noro       20:   K = length(F);
                     21:   for ( I = 0; I < K; I++ ) if ( F[I] ) break;
                     22:   if ( I == K ) return [];
                     23:   if ( type(F[I]) <= 2 ) {
                     24:     // F is a polynimial list
                     25:     F = map(dp_ptod,F,V);
                     26:     F = map(dpm_dptodpm,F,1);
                     27:     N = 1;
                     28:   } else if ( type(F[I]) == 9 ) {
                     29:     // F is a dpoly list
                     30:     F = map(dpm_dptodpm,F,1);
                     31:     N = 1;
                     32:   } else if ( type(F[I]) == 4 ) {
                     33:     // F is a list of poly lists
                     34:     N = length(F[0]);
                     35:     F = map(dpm_ltod,F,V);
                     36:   } else if ( type(F[I]) == 26 ) {
                     37:     // F is a DPM list
                     38:     for ( N = 0, T = F; T != []; T = cdr(T) ) {
                     39:       for ( A = car(T); A; A = dpm_rest(A) ) {
                     40:         N1 = dpm_hp(A);
                     41:         if ( N1 > N ) N = N1;
                     42:       }
                     43:     }
                     44:   } else {
                     45:     error("todpmlist: arugument type is invalid.");
                     46:   }
                     47:   return [F,N];
                     48: }
                     49:
                     50: def module_syz(F,V,H,Ord)
                     51: {
                     52:   if ( type(Weyl=getopt(weyl)) == -1 ) Weyl = 0;
                     53:   if ( type(DP=getopt(dp)) == -1 ) DP = 0;
                     54:   if ( type(F4=getopt(f4)) == -1 ) F4 = 0;
                     55:   dp_ord(Ord);
                     56:   K = length(F);
                     57:   for ( I = 0; I < K; I++ ) if ( F[I] ) break;
                     58:   if ( I == K ) return [[],[],[]];
                     59:   L = todpmlist(F,V);
                     60:   F = L[0]; N = L[1];
                     61:   dp_ord([1,Ord]);
                     62:   B = [];
                     63:   for ( I = 0; I < K; I++ ) {
                     64:     B = cons(F[I]+dpm_dptodpm(dp_ptod(1,V),N+I+1),B);
                     65:   }
                     66:   B = reverse(B);
1.4       noro       67:   if ( H >= 2 ) {
                     68:     // finite field
                     69:     if ( Weyl )
                     70:       G = nd_weyl_gr(B,V,H,[1,Ord]|dp=1);
                     71:     else if ( F4 )
                     72:       G = nd_f4(B,V,H,[1,Ord]|dp=1);
                     73:     else
                     74:       G = nd_gr(B,V,H,[1,Ord]|dp=1);
                     75:   } else {
                     76:     if ( Weyl )
                     77:       G = nd_weyl_gr(B,V,0,[1,Ord]|dp=1,homo=H);
1.7       noro       78:     else {
1.4       noro       79:       Ind = 0;
                     80:       while ( 1 ) {
1.7       noro       81:         if ( F4 )
                     82:           G = nd_f4_trace(B,V,H,-lprime(Ind),[1,Ord]|dp=1);
                     83:         else
                     84:           G = nd_gr_trace(B,V,H,-lprime(Ind),[1,Ord]|dp=1);
1.4       noro       85:         if ( G ) break;
                     86:         else Ind++;
                     87:       }
1.7       noro       88:     }
1.4       noro       89:   }
1.3       noro       90:   G0 = []; S0 = []; Gen0 = [];
                     91:   for ( T = G; T != []; T = cdr(T) ) {
                     92:     H = car(T);
                     93:     if ( dpm_hp(H) > N ) {
                     94:       S0 = cons(dpm_shift(H,N),S0);
                     95:     } else {
                     96:       L = dpm_split(H,N);
                     97:       G0 = cons(L[0],G0);
                     98:       Gen0 = cons(dpm_shift(L[1],N),Gen0);
1.1       noro       99:        }
1.3       noro      100:   }
                    101: #if 0
                    102:   S0 = reverse(S0); G0 = reverse(G0); Gen0 = reverse(Gen0);
                    103: #endif
                    104:   if ( !DP ) {
                    105:     S0 = map(dpm_dtol,S0,V); G0 = map(dpm_dtol,G0,V); Gen0 = map(dpm_dtol,Gen0,V);
                    106:   }
                    107:   return [S0,G0,Gen0];
                    108: }
                    109:
1.4       noro      110: def module_syz_old(F,V,H,O)
                    111: {
                    112:        Weyl = type(getopt(weyl)) != -1 ? 1 : 0;
                    113:        K = length(F);
                    114:        if ( type(F[0]) <= 2 ) {
                    115:                for ( T = [], S = F; S != []; S = cdr(S) )
                    116:                        T = cons([car(S)],T);
                    117:                F = reverse(T);
                    118:        }
                    119:        N = length(F[0]);
                    120:        B = [];
                    121:        for ( I = 0; I < K; I++ ) {
                    122:                E = vector(N+K);
                    123:                for ( J = 0; J < N; J++ ) E[J] = F[I][J];
                    124:                E[N+I] = 1;
                    125:                B = cons(vtol(E),B);
                    126:        }
                    127:        B = reverse(B);
                    128:        if ( H >= 2 ) {
                    129:                if ( Weyl )
                    130:                        G = nd_weyl_gr(B,V,H,[1,O]);
                    131:                else
                    132:                        G = nd_gr(B,V,H,[1,O]);
                    133:        } else {
                    134:                if ( Weyl )
                    135:                        G = nd_weyl_gr_trace(B,V,H,-1,[1,O]);
                    136:                else
                    137:                        G = nd_gr_trace(B,V,H,-1,[1,O]);
                    138:        }
                    139:        G0 = []; S0 = []; Gen0 = [];
                    140:        for ( T = G; T != []; T = cdr(T) ) {
                    141:                H = car(T);
                    142:                for ( J = 0; J < N; J++ ) if ( H[J] ) break;
                    143:                if ( J == N ) {
                    144:                        H1 = vector(K);
                    145:                        for ( J = 0; J < K; J++ ) H1[J] = H[N+J];
                    146:                        S0 = cons(vtol(H1),S0);
                    147:                } else {
                    148:                        H1 = vector(N);
                    149:                        for ( J = 0; J < N; J++ ) H1[J] = H[J];
                    150:                        G0 = cons(vtol(H1),G0);
                    151:                        H1 = vector(K);
                    152:                        for ( J = 0; J < K; J++ ) H1[J] = H[N+J];
                    153:                        Gen0 = cons(vtol(H1),Gen0);
                    154:                }
                    155:        }
                    156:        return [S0,G0,Gen0];
                    157: }
                    158:
1.3       noro      159: def fres(F,V,H,O)
                    160: {
                    161:   if ( type(Weyl=getopt(weyl)) == -1 ) Weyl = 0;
                    162:   if ( type(F4=getopt(f4)) == -1 ) F4 = 0;
                    163:   if ( type(DP=getopt(dp)) == -1 ) DP = 0;
                    164:   L = todpmlist(F,V); F = L[0];
                    165:   R = [F];
                    166:   while ( 1 ) {
                    167:     if ( Weyl )
                    168:       L = module_syz(car(R),V,H,O|weyl=1,dp=1,f4=F4);
                    169:     else
                    170:       L = module_syz(car(R),V,H,O|dp=1,f4=F4);
                    171:     L = L[0];
                    172:     L = ordcheck(L,V);
                    173:     if ( L == [] ) {
1.5       noro      174:       R = reverse(R);
1.3       noro      175:       if ( DP ) return R;
                    176:       else return map(dpmlisttollist,R,V);
                    177:     }
                    178:     R = cons(L,R);
                    179:   }
                    180: }
                    181:
                    182: def minres(F,V,H,O)
                    183: {
                    184:   if ( type(Weyl=getopt(weyl)) == -1 ) Weyl = 0;
                    185:   if ( type(F4=getopt(f4)) == -1 ) F4 = 0;
                    186:   if ( type(DP=getopt(dp)) == -1 ) DP = 0;
                    187:   L = todpmlist(F,V); F = L[0];
                    188:   R = [F];
                    189:   while ( 1 ) {
                    190:     if ( Weyl )
                    191:       L = module_syz(car(R),V,H,O|weyl=1,dp=1,f4=F4);
                    192:     else
                    193:       L = module_syz(car(R),V,H,O|dp=1,f4=F4);
                    194:     L = L[0];
                    195:     L = ordcheck(L,V);
                    196:     if ( L == [] ) break;
                    197:     S = dpm_simplify_syz(L,R[0]);
                    198:     if ( S[0] == [] && S[1] == [] ) {
                    199:       R = cdr(R); break;
                    200:     }
                    201:     R = append([S[0],S[1]],cdr(R));
                    202:     if ( R[0] == [] ) {
                    203:       R = cdr(R); break;
                    204:     }
                    205:   }
1.5       noro      206:   R = reverse(R);
1.3       noro      207:   if ( DP ) return R;
                    208:   else return map(dpmlisttollist,R,V);
                    209: }
                    210:
                    211: def dpmlisttollist(F,V)
                    212: {
                    213:   return map(dpm_dtol,F,V);
                    214: }
                    215:
                    216: def sres(F,V,H,Ord)
                    217: {
                    218:   if ( type(DP=getopt(dp)) == -1 ) DP = 0;
                    219:   dpm_set_schreyer(0);
                    220:   dp_ord(Ord);
                    221:   K = length(F);
                    222:   K = length(F);
                    223:   for ( I = 0; I < K; I++ ) if ( F[I] ) break;
                    224:   if ( I == K ) return [[],[],[]];
                    225:   L = todpmlist(F,V);
                    226:   F = L[0]; N = L[1];
1.6       noro      227: #if 0
1.3       noro      228:   G = nd_gr(F,V,H,[0,Ord]|dp=1);
1.6       noro      229: #else
                    230:   G = nd_gr_trace(F,V,H,1,[0,Ord]|dp=1);
                    231: #endif
1.3       noro      232:   G = reverse(G);
                    233:   R = [G];
                    234:   dp_ord([0,Ord]);
                    235:   while ( 1 ) {
                    236:     S = dpm_schreyer_base(R[0]);
1.7       noro      237:     if ( dp_gr_print() ) print(["length",length(S)]);
1.3       noro      238:     if ( S == [] ) break;
                    239:     else R = cons(S,R);
                    240:   }
                    241:   dp_ord([0,0]);
1.5       noro      242:   R = reverse(R);
1.3       noro      243:   if ( DP ) return R;
                    244:   else return map(dpmlisttollist,R,V);
                    245: }
                    246:
                    247: def minsres(F,V,H,Ord)
                    248: {
                    249:   if ( type(DP=getopt(dp)) == -1 ) DP = 0;
                    250:   R = sres(F,V,H,Ord|dp=1);
1.5       noro      251:   R = ltov(reverse(R));
1.3       noro      252:   M = length(R);
                    253:   for ( I = 0; I < M; I++ ) R[I] = map(dpm_sort,R[I]);
                    254:   R = vtol(R);
                    255:   for ( T = R, U = []; length(T) >= 2; ) {
                    256:     Z = dpm_simplify_syz(T[0],T[1]);
                    257:     U = cons(Z[0],U);
                    258:     T = cons(Z[1],cdr(cdr(T)));
                    259:   }
                    260:   U = cons(T[0],U);
                    261:   if ( DP ) return U;
                    262:   else return map(dpmlisttollist,U,V);
                    263: }
                    264:
                    265: def ordcheck(D,V)
                    266: {
                    267:   B=map(dpm_dtol,D,V);
                    268:   dp_ord([1,0]);
                    269:   D = map(dpm_sort,D);
                    270:   D1 = map(dpm_ltod,B,V);
                    271:   if ( D != D1 )
                    272:     print("afo");
                    273:   return D1;
                    274: }
                    275:
                    276: def complsb(A,B)
                    277: {
                    278:   AL = A[3]; AD = A[4];
                    279:   BL = B[3]; BD = B[4];
                    280:   if ( AD > BD ) return 1;
                    281:   else if ( AD < BD ) return -1;
                    282:   else if ( AL < BL ) return 1;
                    283:   else if ( AL > BL ) return -1;
                    284:   else return 0;
                    285: }
                    286:
                    287: def complsb_sd(A,B)
                    288: {
                    289:   AL = A[3]; AD = A[4];
                    290:   BL = B[3]; BD = B[4];
                    291:   if ( AD-AL > BD-BL ) return 1;
                    292:   else if ( AD-AL < BD-BL ) return -1;
                    293:   else if ( AL > BL ) return 1;
                    294:   else if ( AL < BL ) return -1;
                    295:   else return 0;
                    296: }
                    297:
                    298: def sortlsb(A)
                    299: {
                    300:   S = [];
                    301:   N = length(A);
                    302:   for ( I = 0; I < N; I++ ) S = append(cdr(vtol(A[I])),S);
                    303: //  S = qsort(S,newsyz.complsb_sd);
                    304:   S = qsort(S,newsyz.complsb);
                    305:   return S;
                    306: }
                    307: /* D=[in,sum,deg,lev,ind] */
                    308: /* B=[0 B1 B2 ...] */
                    309: /* C=[0 C1 C2 ...] */
                    310: /* H=[0 H1 H2 ...] */
                    311: /* Z=[0 Z1 Z2 ...] */
                    312: /* Zi=[0 L1 L2 ...] */
                    313: /* Lj=[I1,I2,...] */
                    314: /* One=<<0,...,0>> */
                    315:
                    316: extern Rtime,Stime,Ptime$
                    317:
                    318: // B is sorted wrt positon
                    319: def find_pos(HF,B)
                    320: {
                    321:   Len = length(B);
                    322:   Pos = dpm_hp(HF);
                    323:   First = 0; Last = Len-1;
                    324:   if ( B[First][0] == HF ) return B[First][2];
                    325:   if ( B[Last][0] == HF ) return B[Last][2];
                    326:   while ( 1 ) {
                    327:     Mid = idiv(First+Last,2);
                    328:     PosM = dpm_hp(B[Mid][0]);
                    329:     if ( PosM == Pos ) break;
                    330:     else if ( PosM < Pos )
                    331:       First = Mid;
                    332:     else
                    333:       Last = Mid;
                    334:   }
                    335:   for ( I = Mid; I < Len && dpm_hp(B[I][0]) == Pos; I++ )
                    336:     if ( HF == B[I][0] ) return B[I][2];
                    337:   for ( I = Mid-1; I >= 1 && dpm_hp(B[I][0]) == Pos; I-- )
                    338:     if ( HF == B[I][0] ) return B[I][2];
                    339:   error("find_pos : cannot happen");
                    340: }
                    341:
                    342: def reduce(D,B,Bpos,C,H,Z,K,Kind,G,One,Top)
                    343: {
                    344:   M = D[0]; Ind = D[2]; I = D[3];
                    345:   if ( I == 1 ) {
                    346:     C[I][Ind] = G[Ind];
                    347:     dpm_insert_to_zlist(Z[1],dpm_hp(G[Ind]),Ind);
                    348:     H[I][Ind] = G[Ind];
                    349:   } else {
                    350:     /* M is in F(I-1) => phi(M) is in F(I-2) */
                    351:     /* reduction in F(I-2) */
                    352:     T0 = time()[0];
                    353:     dpm_set_schreyer_level(I-2);
                    354:     CI = C[I-1]; ZI = Z[I-1]; BI = B[I-1]; BposI = Bpos[I-1];
                    355:     Len = size(CI);
                    356:     T=dpm_hc(M); EM=dpm_hp(M);
                    357:     XiM = T*dpm_ht(CI[EM]);
                    358:     EXiM = dpm_hp(XiM);
                    359:     ZIE = ZI[EXiM];
                    360:     for ( ZIE = ZI[EXiM]; ZIE != []; ZIE = cdr(ZIE) ) {
                    361:       J = car(ZIE);
                    362:       if ( J > EM && dpm_redble(XiM,dpm_ht(CI[J])) ) break;
                    363:     }
                    364:     Ptime += time()[0]-T0;
                    365:     T0 = time()[0];
                    366:     QR = dpm_sp_nf(CI,ZI,EM,J|top=Top);
                    367:     Rtime += time()[0]-T0;
                    368:     G = QR[0]; F = QR[1];
                    369:     if ( F ) {
                    370:       HF = dpm_ht(F); EF = dpm_hp(HF);
                    371:       /* find HF in B[I-1] */
                    372:     T0 = time()[0];
                    373:     J = find_pos(HF,BposI);
                    374:     Stime += time()[0]-T0;
                    375:       /* F=Ret[0]*Ret[1] */
                    376:       Ret = dpm_remove_cont(F);
                    377:       CI[J] = Ret[1];
                    378:       Kind[I-1][J] = 2;
                    379:       dpm_insert_to_zlist(ZI,EF,J);
                    380:       dpm_set_schreyer_level(I-1);
                    381:       Tail = -Ret[0]*dpm_dptodpm(One,J);
                    382:       G += Tail;
                    383:       dpm_set_schreyer_level(0);
                    384:       Ret = dpm_remove_cont(G); G = Ret[1];
                    385:       C[I][Ind] = G;
                    386:       Kind[I][Ind] = 1;
                    387:       K[I] = cons([G,Tail],K[I]);
                    388:       dpm_insert_to_zlist(Z[I],EM,Ind);
                    389:       /* level <- 0 */
                    390:       BI[J][3] = 0;
                    391:     } else {
                    392:       Kind[I][Ind] = 0;
                    393:       Ret = dpm_remove_cont(G); G = Ret[1];
                    394:       H[I][Ind] = G;
                    395:       C[I][Ind] = G;
                    396:       dpm_insert_to_zlist(Z[I],EM,Ind);
                    397:     }
                    398:   }
                    399: }
                    400:
1.8     ! noro      401: def comp_lex(A,B)
        !           402: {
        !           403:   HA = dpm_hc(A); HB = dpm_hc(B);
        !           404:   if ( HA > HB ) return 1;
        !           405:   else if ( HA < HB ) return -1;
        !           406:   else return 0;
        !           407: }
        !           408:
1.3       noro      409: def lres_setup(F,V,H,Ord)
                    410: {
1.8     ! noro      411:   if ( type(Lex=getopt(lex)) == -1 ) Lex = 0;
1.3       noro      412:   dpm_set_schreyer(0);
                    413:   dp_ord(Ord);
                    414:   K = length(F);
                    415:   if ( type(F[0]) <= 2 ) {
                    416:     // F is a polynimial list
                    417:     F = map(dp_ptod,F,V);
                    418:     F = map(dpm_dptodpm,F,1);
                    419:     N = 1;
                    420:   } else if ( type(F[0]) == 9 ) {
                    421:     // F is a dpoly list
                    422:     F = map(dpm_dptodpm,F,1);
                    423:     N = 1;
                    424:   } else if ( type(F[0]) == 4 ) {
                    425:     // F is a list of poly lists
                    426:     N = length(F[0]);
                    427:     F = map(dpm_ltod,F,V);
                    428:   } else if ( type(F[0]) == 26 ) {
                    429:     // F is a DPM list
                    430:     for ( N = 0, T = F; T != []; T = cdr(T) ) {
                    431:       for ( A = car(T); A; A = dpm_rest(A) ) {
                    432:         N1 = dpm_hp(A);
                    433:         if ( N1 > N ) N = N1;
                    434:       }
                    435:     }
                    436:   } else {
                    437:     error("lres_setup: arugument type is invalid.");
                    438:   }
1.6       noro      439:   dp_ord([0,Ord]);
                    440:   F = map(dpm_sort,F);
                    441: #if 0
1.3       noro      442:   G = nd_gr(F,V,H,[0,Ord]|dp=1);
1.6       noro      443: #else
                    444:   G = nd_gr_trace(F,V,H,1,[0,Ord]|dp=1);
                    445: #endif
1.8     ! noro      446:   if ( Lex ) {
        !           447:     dp_ord(2);
        !           448:     G = qsort(G,newsyz.comp_lex);
        !           449:   }
1.3       noro      450:   G = reverse(G);
                    451:   dp_ord([0,Ord]);
                    452:   One = dp_ptod(1,V);
                    453:   return [G,One];
                    454: }
                    455:
                    456: def dpm_sort1(L)
                    457: {
                    458:   return [dpm_sort(L[0]),L[1]];
                    459: }
                    460:
                    461: def comp_pos(A,B)
                    462: {
                    463:   PA = dpm_hp(A[0]); PB = dpm_hp(B[0]);
                    464:   if ( PA > PB ) return 1;
                    465:   else if ( PA < PB ) return -1;
                    466:   else return 0;
                    467: }
                    468:
                    469: def lres(F,V,H,Ord)
                    470: {
                    471:   T0 = time();
                    472:   if ( type(Top=getopt(top)) == -1 ) Top = 0;
1.5       noro      473:   if ( type(DP=getopt(dp)) == -1 ) DP = 0;
1.3       noro      474:   if ( type(NoSimpK=getopt(nosimpk)) == -1 ) NoSimpK = 0;
                    475:   if ( type(NoPreProj=getopt(nopreproj)) == -1 ) NoPreProj = 0;
1.8     ! noro      476:   if ( type(Lex=getopt(lex)) == -1 ) Lex = 0;
1.3       noro      477:   Rtime = Stime = Ptime = 0;
                    478:   L = lres_setup(F,V,H,Ord);
                    479:   G = L[0];
                    480:   One = L[1];
1.8     ! noro      481:   F = dpm_schreyer_frame(G|lex=Lex);
1.3       noro      482:   G = ltov(cons(0,L[0]));
                    483:   F = reverse(F);
                    484:   F = ltov(F);
                    485:   N = length(F);
                    486:   for ( I = 0; I < N; I++ ) {
                    487:     FI = F[I]; FI[0] = [];
                    488:     FI = map(ltov,FI);
                    489:     F[I] = FI;
                    490:   }
                    491:   R = sortlsb(F);
                    492:   B = vector(N+1);
                    493:   Bpos = vector(N+1);
                    494:   C = vector(N+1);
                    495:   H = vector(N+1);
                    496:   Z = vector(N+1);
                    497:   K = vector(N+1);
                    498:   L = vector(N+1);
                    499:   K = vector(N+1);
                    500:   D = vector(N+1);
                    501:   Kind = vector(N+1);
                    502:
                    503:   for ( I = 1; I <= N; I++ ) {
                    504:     FI = F[I-1]; Len = length(FI);
                    505:     B[I] = FI;
                    506:     T = vector(Len-1);
                    507:     for ( J = 1; J < Len; J++ ) T[J-1] = FI[J];
                    508:     Bpos[I] = qsort(T,newsyz.comp_pos);
                    509:     C[I] = vector(Len);
                    510:     H[I] = vector(Len);
                    511:     Kind[I] = vector(Len);
                    512:     K[I] = [];
                    513:     Max = 0;
                    514:     for ( J = 1; J < Len; J++ )
                    515:       if ( (Pos = dpm_hp(FI[J][0])) > Max ) Max = Pos;
                    516:     Z[I] = ZI = vector(Max+1);
                    517:     for ( J = 1; J <= Max; J++ ) ZI[J] = [];
                    518:   }
                    519:   T1 = time(); Ftime = T1[0]-T0[0];
                    520:   R = ltov(R); Len = length(R);
1.7       noro      521:   if ( dp_gr_print() ) print(["Len",Len]);
1.3       noro      522:   for ( I = 0, NF = 0; I < Len; I++ ) {
1.7       noro      523:     if ( dp_gr_print() ) {
                    524:       if ( !((I+1)%100) ) print(".",2);
                    525:       if ( !((I+1)%10000) ) print(I+1);
                    526:     }
1.3       noro      527:     if ( !R[I][3] ) continue;
                    528:     NF++;
                    529:     reduce(R[I],B,Bpos,C,H,Z,K,Kind,G,One,Top);
                    530:   }
1.7       noro      531:   if ( dp_gr_print() ) {
                    532:     print("");
                    533:     print(["NF",NF]);
                    534:   }
1.3       noro      535:   T0 = time();
                    536:   dpm_set_schreyer_level(0);
                    537:   D[1] = map(dpm_sort,H[1]);
                    538:   for ( I = 2; I <= N; I++ ) {
                    539:     // HTT = [Head,Tab,TailTop]
                    540:     HTT = create_base_ord(K[I],length(Kind[I-1]));
                    541:     Head = HTT[0]; Tab = HTT[1]; TailTopPos = HTT[2];
                    542:     if ( !NoPreProj )
                    543:       Tab = map(remove_k,map(dpm_sort,Tab),Kind[I-1]);
                    544:     else
                    545:       Tab = map(dpm_sort,Tab);
                    546:     TailTop = dpm_dptodpm(One,TailTopPos);
                    547:     if ( !NoSimpK ) {
1.7       noro      548:       if ( dp_gr_print() ) print("simplify_k "+rtostr(I)+"...",2);
1.3       noro      549:       simplify_k(Head,Tab,TailTop,One);
1.7       noro      550:       if ( dp_gr_print() ) print("done");
1.3       noro      551:     }
                    552:     HI = map(remove_k,map(dpm_sort,H[I]),Kind[I-1]);
                    553:     Len = length(HI);
1.7       noro      554:     if ( dp_gr_print() ) print("simplify_by_k "+rtostr(I)+"...",2);
1.3       noro      555:     D[I] = vector(Len);
                    556:     for ( J = 0; J < Len; J++ ) {
                    557:       D[I][J] = simplify_by_k(HI[J],Tab,TailTop,One);
                    558:       if ( NoPreProj )
                    559:         D[I][J] = remove_k(D[I][J],Kind[I-1]);
                    560:     }
1.7       noro      561:     if ( dp_gr_print() ) print("done");
1.3       noro      562:   }
                    563:   dp_ord([1,0]);
                    564:   T1 = time();
1.7       noro      565:   if ( dp_gr_print() ) print(["Frame",Ftime,"Prep",Ptime,"Reduce",Rtime,"Search",Stime,"Minimalize",T1[0]-T0[0]]);
1.3       noro      566: //  return [C,H,K,Kind,D];
                    567:   D = compress_h(D);
1.5       noro      568:   if ( DP ) return D;
                    569:   else return vtol(map(dpmlisttollist,D,V));
1.3       noro      570: }
                    571:
                    572: def create_base_ord(K,N)
                    573: {
                    574:   Tab = vector(N);
                    575:   Ks = [];
                    576:   for ( T = K; T != []; T = cdr(T) ) {
                    577:     Ks = cons(I=dpm_hp(car(T)[1]),Ks);
                    578:     Tab[I] = car(T)[0];
                    579:   }
                    580:   Others = [];
                    581:   for ( I = N-1; I >= 1; I-- )
                    582:     if ( !Tab[I] ) Others = cons(I,Others);
                    583:   Ks = reverse(Ks);
                    584:   dp_ord([4,append(Ks,Others),0]);
                    585:   return [ltov(Ks),Tab,car(Others)];
                    586: }
                    587:
                    588: /* Head = [i1,i2,...], Tab[i1] = K[0], Tab[i2] = K[1],... */
                    589:
                    590: def simplify_k(Head,Tab,TailTop,One)
                    591: {
                    592:   N = length(Tab);
                    593:   Len = length(Head);
                    594:   for ( I = Len-1; I >= 0; I-- ) {
                    595:     M = 1;
                    596:     R = Tab[Head[I]];
                    597:     H = dpm_hm(R); R = dpm_rest(R);
                    598:     while ( R && dpm_dptodpm(One,dpm_hp(R)) > TailTop ) {
                    599:       Pos = dpm_hp(R); Red = Tab[Pos];
                    600:       CRed = dp_hc(dpm_hc(Red));
                    601:       CR = dpm_extract(R,Pos);
                    602:       L = dpm_remove_cont(CRed*R-CR*Red);
                    603:       R = L[1]; M *= CRed/L[0];
                    604:     }
                    605:     Tab[Head[I]] = M*H+R;
                    606:   }
                    607: }
                    608:
                    609: def simplify_by_k(F,Tab,TailTop,One)
                    610: {
                    611:   R = F;
                    612:   M = 1;
                    613:   while ( R && dpm_dptodpm(One,dpm_hp(R)) > TailTop ) {
                    614:     Pos = dpm_hp(R); Red = Tab[Pos];
                    615:     CRed = dp_hc(dpm_hc(Red));
                    616:     CR = dpm_extract(R,Pos);
                    617:     L = dpm_remove_cont(CRed*R-CR*Red);
                    618:     R = L[1]; M *= CRed/L[0];
                    619:   }
                    620:   return (1/M)*R;
                    621: }
                    622:
                    623: /* Kind[I]=0 => phi(ei)\in H, Kind[I]=1 => phi(ei)\in K */
                    624: def remove_k(F,Kind)
                    625: {
                    626:   R = [];
                    627:   for ( T = F; T; T = dpm_rest(T) )
                    628:     if ( Kind[dpm_hp(T)] != 1 ) R = cons(dpm_hm(T),R);
                    629:   for ( S = 0; R != []; R = cdr(R) )
                    630:     S += car(R);
                    631:   return S;
                    632: }
                    633:
                    634: def remove_k1(F,Kind)
                    635: {
                    636:   return [remove_k(F[0],Kind),F[1]];
                    637: }
                    638:
                    639: def extract_nonzero(A)
                    640: {
                    641:   N = length(A);
                    642:   R = [];
                    643:   for ( C = I = 0; I < N; I++ )
                    644:     if ( A[I] ) R=cons(A[I],R);
                    645:   return reverse(R);;
                    646: }
                    647:
                    648: def nonzero(A)
                    649: {
                    650:   N = length(A);
                    651:   for ( C = I = 0; I < N; I++ )
                    652:     if ( A[I] ) C++;
                    653:   return C;
                    654: }
                    655:
                    656: def phi(C,F)
                    657: {
                    658:   R = 0;
                    659:   for ( T = F; T; T = dpm_rest(T) ) {
                    660:     Coef = dpm_hc(T); Pos = dpm_hp(T);
1.6       noro      661:     R += Coef*C[Pos-1];
1.3       noro      662:   }
                    663:   return R;
1.1       noro      664: }
                    665:
1.3       noro      666: def syz_check(H,I)
1.1       noro      667: {
1.3       noro      668:   HI = H[I];
                    669: //  dpm_set_schreyer_level(I-1);
                    670:   HI1 = H[I+1];
                    671:   Len = size(HI1)[0];
                    672:   for ( J = 1; J < Len; J++ ) {
                    673:     F = HI1[J];
                    674:     R = phi(HI,F);
                    675:     if ( R ) print([J,"NG"]);
                    676:   }
1.1       noro      677: }
                    678:
1.3       noro      679: // for compressed C
                    680: def phi0(C,F)
1.1       noro      681: {
1.3       noro      682:   R = 0;
                    683:   for ( T = F; T; T = dpm_rest(T) ) {
                    684:     Coef = dpm_hc(T); Pos = dpm_hp(T);
                    685:     R += Coef*C[Pos-1];
                    686:   }
                    687:   return R;
1.1       noro      688: }
                    689:
1.3       noro      690: // for compressed H
                    691: def syz_check0(H,I)
1.1       noro      692: {
1.3       noro      693:   HI = H[I];
                    694: //  dpm_set_schreyer_level(I-1);
                    695:   HI1 = H[I+1];
                    696:   Len = length(HI1);
                    697:   for ( J = 0; J < Len; J++ ) {
                    698:     F = HI1[J];
                    699:     R = phi0(HI,F);
                    700:     if ( R ) print([J,"NG"]);
                    701:   }
1.1       noro      702: }
                    703:
1.3       noro      704: // renumber position
                    705: def renumber_pos(F,Tab)
1.1       noro      706: {
1.3       noro      707:   L = [];
                    708:   for ( T = F; T; T = dpm_rest(T) )
                    709:     L = cons(dpm_hm(T),L);
                    710:   R = 0;
                    711:   for ( T = L; T != []; T = cdr(T) )
                    712:     R += dpm_dptodpm(dpm_hc(car(T)),Tab[dpm_hp(car(T))]);
                    713:   return R;
1.1       noro      714: }
                    715:
1.3       noro      716: // compress H1 and renumber H
                    717: def compress(H,H1)
1.1       noro      718: {
1.3       noro      719:   // create index table for H1
                    720:   L1 = length(H1);
                    721:   Tmp = vector(L1);
                    722:   Tab = vector(L1);
                    723:   for ( I = 1, J = 1; I < L1; I++ )
                    724:     if ( H1[I] ) {
                    725:       Tab[I] = J;
                    726:       Tmp[J++] = H1[I];
                    727:     }
                    728:   NH1 = vector(J);
                    729:   for ( I = 1; I < J; I++ )
                    730:     NH1[I] = Tmp[I];
                    731:   if ( H )
                    732:     H = map(renumber_pos,H,Tab);
                    733:   return [H,NH1];
1.1       noro      734: }
                    735:
1.3       noro      736: def compress_h(H)
1.1       noro      737: {
1.3       noro      738:   // H = [0,H[1],...,H[L-1]]
                    739:   L = length(H);
                    740:   NH = vector(L-1);
                    741:   H1 = H[1];
                    742:   for ( I = 2; I < L; I++ ) {
                    743:     R = compress(H[I],H1);
                    744:     H1 = R[0];
                    745:     NH[I-2] = cdr(vtol(R[1]));
                    746:   }
                    747:   R = compress(0,H1);
                    748:   NH[L-2] = cdr(vtol(R[1]));
                    749:   return NH;
1.1       noro      750: }
1.3       noro      751: endmodule$
1.1       noro      752: end$

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