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

Annotation of OpenXM/src/asir-contrib/testing/noro/pd.rr, Revision 1.2

1.1       noro        1: import("gr")$
                      2: module noro_pd$
                      3:
                      4: static GBCheck,F4,Procs,SatHomo$
                      5:
                      6: localf init_procs, kill_procs, syca_dec, syc_dec, find_separating_ideal0$
                      7: localf find_separating_ideal1, find_separating_ideal2$
                      8: localf sy_dec, pseudo_dec, iso_comp, prima_dec$
                      9: localf prime_dec, prime_dec_main, lex_predec1, zprimedec, zprimadec$
                     10: localf complete_qdecomp, partial_qdecomp, partial_qdecomp0, complete_decomp$
                     11: localf partial_decomp, partial_decomp0, zprimacomp, zprimecomp$
                     12: localf fast_gb, elim_gb, ldim, make_mod_subst$
                     13: localf rsgn, find_npos, gen_minipoly, indepset$
                     14: localf maxindep, contraction, ideal_list_intersection, ideal_intersection$
                     15: localf radical_membership, quick_radical_membership, modular_radical_membership$
                     16: localf radical_membership_rep, ideal_product, saturation$
                     17: localf sat, satind, sat_ind, colon$
1.2     ! noro       18: localf ideal_colon, ideal_sat, ideal_inclusion, qd_simp_comp, qd_remove_redundant_comp$
1.1       noro       19: localf remove_redundant_comp, remove_redundant_comp_first, ppart, sq$
                     20: localf lcfactor, compute_deg0, compute_deg, member$
                     21: localf elimination, setintersection, setminus, sep_list$
                     22: localf first_element, comp_tdeg, tdeg, comp_by_ord, comp_by_second$
                     23: localf gbcheck,f4,sathomo$
                     24:
                     25: SatHomo=0$
                     26: GBCheck=1$
                     27:
                     28: #define MAX(a,b) ((a)>(b)?(a):(b))
                     29:
                     30: def gbcheck(A)
                     31: {
                     32:        if ( A ) GBCheck = 1;
                     33:        else GBcheck = -1;
                     34: }
                     35:
                     36: def f4(A)
                     37: {
                     38:        if ( A ) F4 = 1;
                     39:        else F4 = 0;
                     40: }
                     41:
                     42: def sathomo(A)
                     43: {
                     44:        if ( A ) SatHomo = 1;
                     45:        else SatHomo = 0;
                     46: }
                     47:
                     48: def init_procs()
                     49: {
                     50:        if ( type(NoX=getopt(nox)) == -1 ) NoX = 0;
                     51:        if ( !Procs ) {
                     52:                if ( NoX ) {
                     53:                        P0 = ox_launch_nox();
                     54:                        P1 = ox_launch_nox();
                     55:                } else {
                     56:                        P0 = ox_launch();
                     57:                        P1 = ox_launch();
                     58:                }
                     59:                Procs = [P0,P1];
                     60:        }
                     61: }
                     62:
                     63: def kill_procs()
                     64: {
                     65:        if ( Procs ) {
                     66:                ox_shutdown(Procs[0]);
                     67:                ox_shutdown(Procs[1]);
                     68:                Procs = 0;
                     69:        }
                     70: }
                     71:
                     72: /* SYC primary decomositions */
                     73:
                     74: def syca_dec(B,V)
                     75: {
                     76: T00 = time();
                     77:        if ( type(Nolexdec=getopt(nolexdec)) == -1 ) Nolexdec = 0;
                     78:        if ( type(SepIdeal=getopt(sepideal)) == -1 ) SepIdeal = 1;
                     79:        if ( type(NoSimp=getopt(nosimp)) == -1 ) NoSimp = 0;
                     80:        if ( type(Time=getopt(time)) == -1 ) Time = 0;
                     81:        Ord = 0;
                     82:        Gt = G0 = G = fast_gb(B,V,0,Ord);
                     83:        Q0 = Q = []; IntQ0 = IntQ = [1]; First = 1;
                     84:        C = 0;
                     85:
                     86:        Tass = Tiso = Tcolon = Tsep = Tirred = 0;
                     87:        Rass = Riso = Rcolon = Rsep = Rirred = 0;
                     88:        while ( 1 ) {
                     89:                if ( type(Gt[0])==1 ) break;
                     90:                T0 = time();
                     91:                Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec);
                     92:                T1 = time(); Tass += T1[0]-T0[0]+T1[1]-T0[1]; Rass += T1[3]-T0[3];
                     93:                T0 = time();
                     94:                Qt = iso_comp(Gt,Pt,V,Ord);
                     95:                T1 = time(); Tiso += T1[0]-T0[0]+T1[1]-T0[1]; Riso += T1[3]-T0[3];
                     96:                IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord);
                     97:                IntPt = ideal_list_intersection(map(first_element,Pt),V,Ord);
                     98:                if ( First ) {
                     99:                        IntQ0 = IntQ = IntQt; IntP = IntPt; Qi = Qt; First = 0;
                    100:                } else {
                    101:                        IntQ1 = ideal_intersection(IntQ,IntQt,V,Ord);
                    102:                        if ( gb_comp(IntQ,IntQ1) ) {
                    103:                                G = Gt; IntP = IntPt; Q = []; IntQ = [1]; C = 0;
                    104:                                continue;
                    105:                        } else {
                    106:                                IntQ = IntQ1;
                    107:                                IntQ1 = ideal_intersection(IntQ0,IntQt,V,Ord);
                    108:                                if ( !gb_comp(IntQ0,IntQ1) ) {
                    109:                                        IntQ0 = IntQ1;
                    110:                                        Q = append(Qt,Q); Q0 = append(Qt,Q0);
                    111:                                }
                    112:                        }
                    113:                }
                    114:                if ( gb_comp(IntQt,Gt) || gb_comp(IntQ,G) || gb_comp(IntQ0,G0) ) break;
                    115:                T0 = time();
                    116:                C1 = ideal_colon(G,IntQ,V);
                    117:                T1 = time(); Tcolon += T1[0]-T0[0]+T1[1]-T0[1]; Rcolon += T1[3]-T0[3];
                    118:                if ( C && gb_comp(C,C1) ) {
                    119:                        G = Gt; IntP = IntPt; Q = []; IntQ = [1]; C = 0;
                    120:                        continue;
                    121:                } else C = C1;
                    122:                T0 = time();
                    123:                if ( SepIdeal == 0 )
                    124:                        Ok = find_separating_ideal0(C,G,IntQ,IntP,V,Ord);
                    125:                else if ( SepIdeal == 1 )
                    126:                        Ok = find_separating_ideal1(C,G,IntQ,IntP,V,Ord);
                    127:                else if ( SepIdeal == 2 )
                    128:                        Ok = find_separating_ideal2(C,G,IntQ,IntP,V,Ord);
                    129:                G1 = append(Ok,G);
                    130:                Gt1 = fast_gb(G1,V,0,Ord);
                    131:                T1 = time(); Tsep += T1[0]-T0[0]+T1[1]-T0[1]; Rsep += T1[3]-T0[3];
                    132: #if 0
                    133:                if ( ideal_inclusion(Gt1,Gt,V,Ord) ) {
                    134:                        G = Gt; IntP = IntPt; Q = []; IntQ = [1]; C = 0;
                    135:                } else
                    136: #endif
                    137:                        Gt = Gt1;
                    138:        }
                    139:        T0 = time();
                    140:        if ( !NoSimp ) Q1 = qd_remove_redundant_comp(G0,Qi,Q0,V,Ord);
                    141:        else Q1 = Q0;
                    142:        if ( Time ) {
                    143:                T1 = time(); Tirred += T1[0]-T0[0]+T1[1]-T0[1]; Rirred += T1[3]-T0[3];
                    144:                Tall = T1[0]-T00[0]+T1[1]-T00[1]; Rall += T1[3]-T00[3];
                    145:                print(["total",Tall,"ass",Tass,"iso",Tiso, "colon",Tcolon,"sep",Tsep,"irred",Tirred]);
                    146:                print(["Rtotal",Rall,"Rass",Rass,"Riso",Riso, "Rcolon",Rcolon,"Rsep",Rsep,"Rirred",Rirred]);
                    147:                print(["iso",length(Qi),"emb",length(Q0),"->",length(Q1)]);
                    148:        }
                    149:        return [Qi,Q1];
                    150: }
                    151:
                    152: def syc_dec(B,V)
                    153: {
                    154: T00 = time();
                    155:        if ( type(Nolexdec=getopt(nolexdec)) == -1 ) Nolexdec = 0;
                    156:        if ( type(SepIdeal=getopt(sepideal)) == -1 ) SepIdeal = 1;
                    157:        if ( type(NoSimp=getopt(nosimp)) == -1 ) NoSimp = 0;
                    158:        if ( type(Time=getopt(time)) == -1 ) Time = 0;
                    159:        Ord = 0;
                    160:        G = fast_gb(B,V,0,Ord);
                    161:        Q = []; IntQ = [1]; Gt = G; First = 1;
                    162:        Tass = Tiso = Tcolon = Tsep = Tirred = 0;
                    163:        Rass = Riso = Rcolon = Rsep = Rirred = 0;
                    164:        while ( 1 ) {
                    165:                if ( type(Gt[0])==1 ) break;
                    166:                T0 = time();
                    167:                Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec);
                    168:                T1 = time(); Tass += T1[0]-T0[0]+T1[1]-T0[1]; Rass += T1[3]-T0[3];
                    169:                T0 = time();
                    170:                Qt = iso_comp(Gt,Pt,V,Ord);
                    171:                T1 = time(); Tiso += T1[0]-T0[0]+T1[1]-T0[1]; Riso += T1[3]-T0[3];
                    172:                IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord);
                    173:                IntPt = ideal_list_intersection(map(first_element,Pt),V,Ord);
                    174:                if ( First ) {
                    175:                        IntQ = IntQt; Qi = Qt; First = 0;
                    176:                } else {
                    177:                        IntQ1 = ideal_intersection(IntQ,IntQt,V,Ord);
                    178:                        if ( !gb_comp(IntQ1,IntQ) )
                    179:                                Q = append(Qt,Q);
                    180:                }
                    181:                if ( gb_comp(IntQ,G) || gb_comp(IntQt,Gt) )
                    182:                        break;
                    183:                T0 = time();
                    184:                C = ideal_colon(Gt,IntQt,V);
                    185:                T1 = time(); Tcolon += T1[0]-T0[0]+T1[1]-T0[1]; Rcolon += T1[3]-T0[3];
1.2     ! noro      186:                T0 = time();
1.1       noro      187:                if ( SepIdeal == 0 )
                    188:                        Ok = find_separating_ideal0(C,Gt,IntQt,IntPt,V,Ord);
                    189:                else if ( SepIdeal == 1 )
                    190:                        Ok = find_separating_ideal1(C,Gt,IntQt,IntPt,V,Ord);
                    191:                else if ( SepIdeal == 2 )
                    192:                        Ok = find_separating_ideal2(C,Gt,IntQt,IntPt,V,Ord);
                    193:                G1 = append(Ok,Gt);
                    194:                Gt = fast_gb(G1,V,0,Ord);
                    195:                T1 = time(); Tsep += T1[0]-T0[0]+T1[1]-T0[1]; Rsep += T1[3]-T0[3];
                    196:        }
                    197:        T0 = time();
                    198:        if ( !NoSimp ) Q1 = qd_remove_redundant_comp(G,Qi,Q,V,Ord);
                    199:        else Q1 = Q;
                    200:        T1 = time(); Tirred += T1[0]-T0[0]+T1[1]-T0[1]; Rirred += T1[3]-T0[3];
                    201:        Tall = T1[0]-T00[0]+T1[1]-T00[1]; Rall += T1[3]-T00[3];
                    202:        if ( Time ) {
                    203:                print(["total",Tall,"ass",Tass,"iso",Tiso, "colon",Tcolon,"sep",Tsep,"irred",Tirred]);
                    204:                print(["Rtotal",Rall,"Rass",Rass,"Riso",Riso, "Rcolon",Rcolon,"Rsep",Rsep,"Rirred",Rirred]);
                    205:                print(["iso",length(Qi),"emb",length(Q),"->",length(Q1)]);
                    206:        }
                    207:        return [Qi,Q1];
                    208: }
                    209:
                    210: /* C=G:Q, Rad=rad(Q), return J s.t. Q cap (G+J) = G */
                    211:
                    212: def find_separating_ideal0(C,G,Q,Rad,V,Ord) {
                    213:        for ( CI = C, I = 1; ; I++ ) {
                    214:                for ( T = CI, S = []; T != []; T = cdr(T) )
                    215:                        if ( nd_nf(car(T),Q,V,Ord,0) ) S = cons(car(T),S);
                    216:                if ( S == [] )
                    217:                        error("find_separating_ideal0 : cannot happen");
                    218:                G1 = append(S,G);
                    219:                Int = ideal_intersection(G1,Q,V,Ord);
                    220:                /* check whether (Q cap (G+S)) = G */
                    221:                if ( gb_comp(Int,G) ) return reverse(S);
                    222:                CI = ideal_product(CI,C,V);
                    223:        }
                    224: }
                    225:
                    226: def find_separating_ideal1(C,G,Q,Rad,V,Ord) {
                    227:        for ( T = C, S = []; T != []; T = cdr(T) )
                    228:                if ( nd_nf(car(T),Q,V,Ord,0) ) S = cons(car(T),S);
                    229:        if ( S == [] )
                    230:                error("find_separating_ideal1 : cannot happen");
                    231:        G1 = append(S,G);
                    232:        Int = ideal_intersection(G1,Q,V,Ord);
                    233:        /* check whether (Q cap (G+S)) = G */
                    234:        if ( gb_comp(Int,G) ) return reverse(S);
                    235:
                    236:        C = qsort(S,comp_tdeg);
                    237:        for ( T = C, S = []; T != []; T = cdr(T) ) {
                    238:                if ( !nd_nf(car(T),Rad,V,Ord,0) ) continue;
                    239:                Ui = U = car(T);
                    240:                for ( I = 1; ; I++ ) {
                    241:                        G1 = cons(Ui,G);
                    242:                        Int = ideal_intersection(G1,Q,V,Ord);
                    243:                        if ( gb_comp(Int,G) ) break;
                    244:                        else
                    245:                                Ui = nd_nf(Ui*U,G,V,Ord,0);
                    246:                }
                    247:                if ( length(S) ) {
                    248:                        G1 = append(cons(Ui,S),G);
                    249:                        Int = ideal_intersection(G1,Q,V,Ord);
                    250:                        if ( !gb_comp(Int,G) )
                    251:                                break;
                    252:                }
                    253:                S = cons(Ui,S);
                    254:        }
                    255:        return reverse(S);
                    256: }
                    257:
                    258: def find_separating_ideal2(C,G,Q,Rad,V,Ord) {
                    259:        for ( T = C, S = []; T != []; T = cdr(T) )
                    260:                if ( nd_nf(car(T),Q,V,Ord,0) ) S = cons(car(T),S);
                    261:        if ( S == [] )
                    262:                error("find_separating_ideal2 : cannot happen");
                    263:        G1 = append(S,G);
                    264:        Int = ideal_intersection(G1,Q,V,Ord);
                    265:        /* check whether (Q cap (G+S)) = G */
                    266:        if ( gb_comp(Int,G) ) return reverse(S);
                    267:
                    268:        C = qsort(C,comp_tdeg);
                    269:        for ( T = C, S = []; T != []; T = cdr(T) ) {
                    270:                if ( !nd_nf(car(T),Rad,V,Ord,0) ) continue;
                    271:                Ui = U = car(T);
                    272:                for ( I = 1; ; I++ ) {
                    273:                        G1 = cons(Ui,G);
                    274:                        Int = ideal_intersection(G1,Q,V,Ord);
                    275:                        if ( gb_comp(Int,G) ) break;
                    276:                        else
                    277:                                Ui = nd_nf(Ui*U,G,V,Ord,0);
                    278:                }
                    279:                S = cons(Ui,S);
                    280:        }
                    281:        S = reverse(S);
                    282:        Len = length(S);
                    283:        Ok = [S[0]];
                    284:        if ( Len > 1 ) {
                    285:                K = 2;
                    286:                while ( 1 ) {
                    287:                        for ( St = [], I = 0; I < K; I++ ) St = cons(S[I],St);
                    288:                        G1 = append(St,G);
                    289:                        Int = ideal_intersection(G1,Q,V,Ord);
                    290:                        if ( !gb_comp(Int,G) ) break;
                    291:                        Ok = St;
                    292:                        if ( K == Len ) break;
                    293:                        else {
                    294:                                K = 2*K;
                    295:                                if ( K > Len ) K = Len;
                    296:                        }
                    297:                }
                    298:        }
                    299:        return Ok;
                    300: }
                    301:
                    302: /* SY primary decompsition */
                    303:
                    304: def sy_dec(B,V)
                    305: {
                    306:        if ( type(Nolexdec=getopt(nolexdec)) == -1 ) Nolexdec = 0;
                    307:        Ord = 0;
                    308:        G = fast_gb(B,V,0,Ord);
                    309:        Q = [];
                    310:        IntQ = [1];
                    311:        Gt = G;
                    312:        First = 1;
                    313:        while ( 1 ) {
                    314:                if ( type(Gt[0]) == 1 ) break;
                    315:                Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec);
                    316:                L = pseudo_dec(Gt,Pt,V,Ord);
                    317:                Qt = L[0]; Rt = L[1]; St = L[2];
                    318:                IntQt = ideal_list_intersection(Qt,V,Ord);
                    319:                if ( First ) {
                    320:                        IntQ = IntQt;
                    321:                        Qi = Qt;
                    322:                        First = 0;
                    323:                } else {
                    324:                        IntQ = ideal_intersection(IntQ,IntQt,V,Ord);
                    325:                        Q = append(Qt,Q);
                    326:                }
                    327:                if ( gb_comp(IntQ,G) ) break;
                    328:                for ( T = Rt; T != []; T = cdr(T) ) {
                    329:                        if ( type(car(T)[0]) == 1 ) continue;
                    330:                        U = sy_dec(car(T),V|nolexdec=Nolexdec);
                    331:                        IntQ = ideal_list_intersection(cons(IntQ,U),V,Ord);
                    332:                        Q = append(U,Q);
                    333:                        if ( gb_comp(IntQ,G) ) break;
                    334:                }
                    335:                Gt = fast_gb(append(Gt,St),V,0,Ord);
                    336:        }
                    337:        Q = remove_redundant_comp(G,Qi,Q,V,Ord);
                    338:        return append(Qi,Q);
                    339: }
                    340:
                    341: def pseudo_dec(G,L,V,Ord)
                    342: {
                    343:        N = length(L);
                    344:        S = vector(N);
                    345:        Q = vector(N);
                    346:        R = vector(N);
                    347:        L0 = map(first_element,L);
                    348:        for ( I = 0; I < N; I++ ) {
                    349:                LI = setminus(L0,[L0[I]]);
                    350:                PI = ideal_list_intersection(LI,V,Ord);
                    351:                PI = qsort(PI,comp_tdeg);
                    352:                for ( T = PI; T != []; T = cdr(T) )
                    353:                        if ( p_nf(car(T),L0[I],V,Ord) ) break;
                    354:                if ( T == [] ) error("separator : cannot happen");
                    355:                SI = sat_ind(G,car(T),V);
                    356:                QI = SI[0];
                    357:                S[I] = car(T)^SI[1];
                    358:                PV = L[I][1];
                    359:                V0 = setminus(V,PV);
                    360: #if 0
                    361:                GI = fast_gb(QI,append(V0,PV),0,
                    362:                        [[Ord,length(V0)],[Ord,length(PV)]]);
                    363: #else
                    364:                GI = fast_gb(QI,append(V0,PV),0,
                    365:                        [[2,length(V0)],[Ord,length(PV)]]);
                    366: #endif
                    367:                LCFI = lcfactor(GI,V0,Ord);
                    368:                for ( F = 1, T = LCFI, Gt = QI; T != []; T = cdr(T) ) {
                    369:                        St = sat_ind(Gt,T[0],V);
                    370:                        Gt = St[0]; F *= T[0]^St[1];
                    371:                }
                    372:                Q[I] = Gt;
                    373:                R[I] = fast_gb(cons(F,QI),V,0,Ord);
                    374:        }
                    375:        return [vtol(Q),vtol(R),vtol(S)];
                    376: }
                    377:
                    378: def iso_comp(G,L,V,Ord)
                    379: {
                    380:        N = length(L);
                    381:        S = vector(N);
                    382:        Ind = vector(N);
                    383:        Q = vector(N);
                    384:        L0 = map(first_element,L);
                    385:        for ( I = 0; I < N; I++ ) {
                    386:                LI = setminus(L0,[L0[I]]);
                    387:                PI = ideal_list_intersection(LI,V,Ord);
                    388:                for ( T = PI; T != []; T = cdr(T) )
                    389:                        if ( p_nf(car(T),L0[I],V,Ord) ) break;
                    390:                if ( T == [] ) error("separator : cannot happen");
                    391:                S[I] = car(T);
                    392:                QI = sat(G,S[I],V);
                    393:                PV = L[I][1];
                    394:                V0 = setminus(V,PV);
                    395:                GI = elim_gb(QI,V0,PV,0,[[0,length(V0)],[0,length(PV)]]);
                    396:                Q[I] = [contraction(GI,V0),L0[I]];
                    397:        }
                    398:        return vtol(Q);
                    399: }
                    400:
                    401: /* GTZ primary decompsition */
                    402:
                    403: def prima_dec(B,V)
                    404: {
                    405:        G = nd_gr_trace(B,V,1,GBCheck,0);
                    406:        G0 = G;
                    407:        IntP = [1];
                    408:        QD = [];
                    409:        while ( 1 ) {
                    410:                if ( ideal_inclusion(IntP,G0,V,0) )
                    411:                        return QD;
                    412:                W = maxindep(G,V,0); NP = length(W);
                    413:                V0 = setminus(V,W); N = length(V0);
                    414:                V1 = append(V0,W);
                    415:                G1 = fast_gb(G,V1,0,[[0,N],[0,NP]]);
                    416:                LCF = lcfactor(G1,V0,0);
                    417:                L = zprimacomp(G,V0);
                    418:                F = 1;
                    419:                for ( T = LCF, G2 = G1; T != []; T = cdr(T) ) {
                    420:                        S = sat_ind(G2,T[0],V1);
                    421:                        G2 = S[0]; F *= T[0]^S[1];
                    422:                }
                    423:                for ( T = L, QL = []; T != []; T = cdr(T) )
                    424:                        QL = cons(car(T)[0],QL);
                    425:                Int = ideal_list_intersection(QL,V,0);
                    426:                IntP = ideal_intersection(IntP,Int,V,0);
                    427:                QD = append(QD,L);
                    428:                F = p_nf(F,G,V,0);
                    429:                G = cons(F,G);
                    430:        }
                    431: }
                    432:
                    433: /* SL prime decomposition */
                    434:
                    435: def prime_dec(B,V)
                    436: {
                    437:        if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
                    438:        if ( type(NoLexDec=getopt(nolexdec)) == -1 ) NoLexDec = 0;
                    439:        B = map(sq,B);
                    440:        if ( !NoLexDec )
                    441:                PD = lex_predec1(B,V);
                    442:        else
                    443:                PD = [B];
                    444:        G = ideal_list_intersection(PD,V,0);
                    445:        PD = remove_redundant_comp(G,[],PD,V,0);
                    446:        R = [];
                    447:        for ( T = PD; T != []; T = cdr(T) )
                    448:                R = append(prime_dec_main(car(T),V|indep=Indep),R);
                    449:        if ( Indep ) {
                    450:                G = ideal_list_intersection(map(first_element,R),V,0);
                    451:                R = remove_redundant_comp_first(G,R,V,0);
                    452:        } else {
                    453:                G = ideal_list_intersection(R,V,0);
                    454:                R = remove_redundant_comp(G,[],R,V,0);
                    455:        }
                    456:        return R;
                    457: }
                    458:
                    459: def prime_dec_main(B,V)
                    460: {
                    461:        if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
                    462:        G = nd_gr_trace(B,V,1,GBCheck,0);
                    463:        IntP = [1];
                    464:        PD = [];
                    465:        while ( 1 ) {
                    466:                /* rad(G) subset IntP */
                    467:                /* check if IntP subset rad(G) */
                    468:                for ( T = IntP; T != []; T = cdr(T) ) {
                    469:                        if ( (GNV = modular_radical_membership(car(T),G,V)) ) {
                    470:                                F = car(T);
                    471:                                break;
                    472:                        }
                    473:                }
                    474:                if ( T == [] ) return PD;
                    475:
                    476:                /* GNV = [GB(<NV*F-1,G>),NV] */
                    477:                G1 = nd_gr_trace(GNV[0],cons(GNV[1],V),1,GBCheck,[[0,1],[0,length(V)]]);
                    478:                G0 = elimination(G1,V);
                    479:                PD0 = zprimecomp(G0,V,Indep);
                    480:                if ( Indep ) {
                    481:                        Int = ideal_list_intersection(PD0[0],V,0);
                    482:                        IndepSet = PD0[1];
                    483:                        for ( PD1 = [], T = PD0[0]; T != []; T = cdr(T) )
                    484:                                PD1 = cons([car(T),IndepSet],PD1);
                    485:                        PD = append(PD,reverse(PD1));
                    486:                } else {
                    487:                        Int = ideal_list_intersection(PD0,V,0);
                    488:                        PD = append(PD,PD0);
                    489:                }
                    490:                IntP = ideal_intersection(IntP,Int,V,0);
                    491:        }
                    492: }
                    493:
                    494: /* pre-decomposition */
                    495:
                    496: def lex_predec1(B,V)
                    497: {
                    498:        G = nd_gr_trace(B,V,1,GBCheck,2);
                    499:        for ( T = G; T != []; T = cdr(T) ) {
                    500:                F = fctr(car(T));
                    501:                if ( length(F) > 2 || length(F) == 2 && F[1][1] > 1 ) {
                    502:                        for ( R = [], S = cdr(F); S != []; S = cdr(S) ) {
                    503:                                Ft = car(S)[0];
                    504:                                Gt = map(ptozp,map(p_nf,G,[Ft],V,0));
                    505:                                R1 = nd_gr_trace(cons(Ft,Gt),V,1,GBCheck,0);
                    506:                                R = cons(R1,R);
                    507:                        }
                    508:                        return R;
                    509:                }
                    510:        }
                    511:        return [G];
                    512: }
                    513:
                    514: /* zero-dimensional prime/primary decomosition */
                    515:
                    516: def zprimedec(B,V,Mod)
                    517: {
                    518:        L = partial_decomp(B,V,Mod);
                    519:        P = L[0]; NP = L[1];
                    520:        R = [];
                    521:        for ( ; P != []; P = cdr(P) ) R = cons(car(car(P)),R);
                    522:        for ( T = NP; T != []; T = cdr(T) ) {
                    523:                R1 = complete_decomp(car(T),V,Mod);
                    524:                R = append(R1,R);
                    525:        }
                    526:        return R;
                    527: }
                    528:
                    529: def zprimadec(B,V,Mod)
                    530: {
                    531:        L = partial_qdecomp(B,V,Mod);
                    532:        Q = L[0]; NQ = L[1];
                    533:        R = [];
                    534:        for ( ; Q != []; Q = cdr(Q) ) {
                    535:                T = car(Q); R = cons([T[0],T[1]],R);
                    536:        }
                    537:        for ( T = NQ; T != []; T = cdr(T) ) {
                    538:                R1 = complete_qdecomp(car(T),V,Mod);
                    539:                R = append(R1,R);
                    540:        }
                    541:        return R;
                    542: }
                    543:
                    544: def complete_qdecomp(GD,V,Mod)
                    545: {
                    546:        GQ = GD[0]; GP = GD[1]; D = GD[2];
                    547:        W = vars(GP);
                    548:        PV = setminus(W,V);
                    549:        N = length(V); PN = length(PV);
                    550:        U = find_npos([GP,D],V,PV,Mod);
                    551:        NV = ttttt;
                    552:        M = gen_minipoly(cons(NV-U,GQ),cons(NV,V),PV,0,NV,Mod);
                    553:        M = ppart(M,NV,Mod);
                    554:        MF = Mod ? modfctr(M) : fctr(M);
                    555:        R = [];
                    556:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                    557:                S = car(T);
                    558:                Mt = subst(S[0],NV,U);
                    559:                GP1 = fast_gb(cons(Mt,GP),W,Mod,0);
                    560:                GQ1 = fast_gb(cons(Mt^S[1],GQ),W,Mod,0);
                    561:                if ( PV != [] ) {
                    562:                        GP1 = elim_gb(GP1,V,PV,Mod,[[0,N],[0,PN]]);
                    563:                        GQ1 = elim_gb(GQ1,V,PV,Mod,[[0,N],[0,PN]]);
                    564:                }
                    565:                R = cons([GQ1,GP1],R);
                    566:        }
                    567:        return R;
                    568: }
                    569:
                    570: def partial_qdecomp(B,V,Mod)
                    571: {
                    572:        Elim = (Elim=getopt(elim))&&type(Elim)!=-1 ? 1 : 0;
                    573:        N = length(V);
                    574:        W = vars(B);
                    575:        PV = setminus(W,V);
                    576:        NP = length(PV);
                    577:        W = append(V,PV);
                    578:        if ( Elim && PV != [] ) Ord = [[0,N],[0,NP]];
                    579:        else Ord = 0;
                    580:        if ( Mod )
                    581:                B = nd_f4(B,W,Mod,Ord);
                    582:        else
                    583:                B = nd_gr_trace(B,W,1,GBCheck,Ord);
                    584:        Q = []; NQ = [[B,B,vector(N+1)]];
                    585:        for ( I = length(V)-1; I >= 0; I-- ) {
                    586:                NQ1 = [];
                    587:                for ( T = NQ; T != []; T = cdr(T) ) {
                    588:                        L = partial_qdecomp0(car(T),V,PV,Ord,I,Mod);
                    589:                        Q = append(L[0],Q);
                    590:                        NQ1 = append(L[1],NQ1);
                    591:                }
                    592:                NQ = NQ1;
                    593:        }
                    594:        return [Q,NQ];
                    595: }
                    596:
                    597: def partial_qdecomp0(GD,V,PV,Ord,I,Mod)
                    598: {
                    599:        GQ = GD[0]; GP = GD[1]; D = GD[2];
                    600:        N = length(V); PN = length(PV);
                    601:        W = append(V,PV);
                    602:        VI = V[I];
                    603:        M = gen_minipoly(GQ,V,PV,Ord,VI,Mod);
                    604:        M = ppart(M,VI,Mod);
                    605:        if ( Mod )
                    606:                MF = modfctr(M,Mod);
                    607:        else
                    608:                MF = fctr(M);
                    609:        Q = []; NQ = [];
                    610:        if ( length(MF) == 2 && MF[1][1] == 1 ) {
                    611:                D1 = D*1; D1[I] = M;
                    612:                GQelim = elim_gb(GQ,V,PV,Mod,Ord);
                    613:                GPelim = elim_gb(GP,V,PV,Mod,Ord);
                    614:                LD = ldim(GQelim,V);
                    615:                if ( deg(M,VI) == LD )
                    616:                        Q = cons([GQelim,GPelim,D1],Q);
                    617:                else
                    618:                        NQ = cons([GQelim,GPelim,D1],NQ);
                    619:                return [Q,NQ];
                    620:        }
                    621:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                    622:                S = car(T); Mt = S[0]; D1 = D*1; D1[I] = Mt;
                    623:
                    624:                GQ1 = fast_gb(cons(Mt^S[1],GQ),W,Mod,Ord);
                    625:                GQelim = elim_gb(GQ1,V,PV,Mod,Ord);
                    626:                GP1 = fast_gb(cons(Mt,GP),W,Mod,Ord);
                    627:                GPelim = elim_gb(GP1,V,PV,Mod,Ord);
                    628:
                    629:                D1[N] = LD = ldim(GPelim,V);
                    630:
                    631:                for ( J = 0; J < N; J++ )
                    632:                        if ( D1[J] && deg(D1[J],V[J]) == LD ) break;
                    633:                if ( J < N )
                    634:                        Q = cons([GQelim,GPelim,D1],Q);
                    635:                else
                    636:                        NQ = cons([GQelim,GPelim,D1],NQ);
                    637:        }
                    638:        return [Q,NQ];
                    639: }
                    640:
                    641: def complete_decomp(GD,V,Mod)
                    642: {
                    643:        G = GD[0]; D = GD[1];
                    644:        W = vars(G);
                    645:        PV = setminus(W,V);
                    646:        N = length(V); PN = length(PV);
                    647:        U = find_npos(GD,V,PV,Mod);
                    648:        NV = ttttt;
                    649:        M = gen_minipoly(cons(NV-U,G),cons(NV,V),PV,0,NV,Mod);
                    650:        M = ppart(M,NV,Mod);
                    651:        MF = Mod ? modfctr(M) : fctr(M);
                    652:        if ( length(MF) == 2 ) return [G];
                    653:        R = [];
                    654:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                    655:                Mt = subst(car(car(T)),NV,U);
                    656:                G1 = fast_gb(cons(Mt,G),W,Mod,0);
                    657:                if ( PV != [] ) G1 = elim_gb(G1,V,PV,Mod,[[0,N],[0,PN]]);
                    658:                R = cons(G1,R);
                    659:        }
                    660:        return R;
                    661: }
                    662:
                    663: def partial_decomp(B,V,Mod)
                    664: {
                    665:        Elim = (Elim=getopt(elim))&&type(Elim)!=-1 ? 1 : 0;
                    666:        N = length(V);
                    667:        W = vars(B);
                    668:        PV = setminus(W,V);
                    669:        NP = length(PV);
                    670:        W = append(V,PV);
                    671:        if ( Elim && PV != [] ) Ord = [[0,N],[0,NP]];
                    672:        else Ord = 0;
                    673:        if ( Mod )
                    674:                B = nd_f4(B,W,Mod,Ord);
                    675:        else
                    676:                B = nd_gr_trace(B,W,1,GBCheck,Ord);
                    677:        P = []; NP = [[B,vector(N+1)]];
                    678:        for ( I = length(V)-1; I >= 0; I-- ) {
                    679:                NP1 = [];
                    680:                for ( T = NP; T != []; T = cdr(T) ) {
                    681:                        L = partial_decomp0(car(T),V,PV,Ord,I,Mod);
                    682:                        P = append(L[0],P);
                    683:                        NP1 = append(L[1],NP1);
                    684:                }
                    685:                NP = NP1;
                    686:        }
                    687:        return [P,NP];
                    688: }
                    689:
                    690: def partial_decomp0(GD,V,PV,Ord,I,Mod)
                    691: {
                    692:        G = GD[0]; D = GD[1];
                    693:        N = length(V); PN = length(PV);
                    694:        W = append(V,PV);
                    695:        VI = V[I];
                    696:        M = gen_minipoly(G,V,PV,Ord,VI,Mod);
                    697:        M = ppart(M,VI,Mod);
                    698:        if ( Mod )
                    699:                MF = modfctr(M,Mod);
                    700:        else
                    701:                MF = fctr(M);
                    702:        if ( length(MF) == 2 && MF[1][1] == 1 ) {
                    703:                D1 = D*1;
                    704:                D1[I] = M;
                    705:                Gelim = elim_gb(G,V,PV,Mod,Ord);
                    706:                D1[N] = LD = ldim(Gelim,V);
                    707:                GD1 = [Gelim,D1];
                    708:                for ( J = 0; J < N; J++ )
                    709:                        if ( D1[J] && deg(D1[J],V[J]) == LD )
                    710:                                return [[GD1],[]];
                    711:                return [[],[GD1]];
                    712:        }
                    713:        P = []; NP = [];
                    714:        GI = elim_gb(G,V,PV,Mod,Ord);
                    715:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                    716:                Mt = car(car(T));
                    717:                D1 = D*1;
                    718:                D1[I] = Mt;
                    719:                GIt = map(p_nf,GI,[Mt],V,Ord);
                    720:                G1 = cons(Mt,GIt);
                    721:                Gelim = elim_gb(G1,V,PV,Mod,Ord);
                    722:                D1[N] = LD = ldim(Gelim,V);
                    723:                for ( J = 0; J < N; J++ )
                    724:                        if ( D1[J] && deg(D1[J],V[J]) == LD ) break;
                    725:                if ( J < N )
                    726:                        P = cons([Gelim,D1],P);
                    727:                else
                    728:                        NP = cons([Gelim,D1],NP);
                    729:        }
                    730:        return [P,NP];
                    731: }
                    732:
                    733: /* prime/primary components over rational function field */
                    734:
                    735: def zprimacomp(G,V) {
                    736:        L = zprimadec(G,V,0);
                    737:        R = [];
                    738:        dp_ord(0);
                    739:        for ( T = L; T != []; T = cdr(T) ) {
                    740:                S = car(T);
                    741:                UQ = contraction(S[0],V);
                    742:                UP = contraction(S[1],V);
                    743:                R = cons([UQ,UP],R);
                    744:        }
                    745:        return R;
                    746: }
                    747:
                    748: def zprimecomp(G,V,Indep) {
                    749:        W = maxindep(G,V,0);
                    750:        V0 = setminus(V,W);
                    751:        V1 = append(V0,W);
                    752: #if 0
                    753:        O1 = [[0,length(V0)],[0,length(W)]];
                    754:        G1 = nd_gr_trace(G,V1,1,GBCheck,O1);
                    755:        dp_ord(0);
                    756: #else
                    757:        G1 = G;
                    758: #endif
                    759:        PD = zprimedec(G1,V0,0);
                    760:        dp_ord(0);
                    761:        R = [];
                    762:        for ( T = PD; T != []; T = cdr(T) ) {
                    763:                U = contraction(car(T),V0);
                    764:                R = cons(U,R);
                    765:        }
                    766:        if ( Indep ) return [R,W];
                    767:        else return R;
                    768: }
                    769:
                    770: def fast_gb(B,V,Mod,Ord)
                    771: {
                    772:        NoRA = (NoRA=getopt(nora))&&type(NoRA)!=-1 ? 1 : 0;
                    773:        if ( Mod )
                    774:                G = nd_f4(B,V,Mod,Ord|nora=NoRA);
                    775:        else {
                    776:                if ( F4 )
                    777:                        G = map(ptozp,f4_chrem(B,V,Ord));
                    778:                else
                    779:                        G = nd_gr_trace(B,V,1,GBCheck,Ord|nora=NoRA);
                    780:        }
                    781:        return G;
                    782: }
                    783:
                    784:
                    785: def elim_gb(G,V,PV,Mod,Ord)
                    786: {
                    787:        N = length(V); PN = length(PV);
                    788:        O1 = [[0,N],[0,PN]];
                    789:        if ( Ord == O1 )
                    790:                Ord = Ord[0][0];
                    791:        if ( Mod ) /* XXX */
                    792:                G = dp_gr_mod_main(G,V,0,Mod,Ord);
                    793:        else if ( Procs ) {
                    794:                Arg0 = ["nd_gr_trace",G,V,1,GBCheck,Ord];
                    795:                Arg1 = ["nd_gr_trace_rat",G,V,PV,1,GBCheck,O1,Ord];
                    796:                G = competitive_exec(Procs,Arg0,Arg1);
                    797:        } else
                    798:                G = nd_gr_trace(G,V,1,GBCheck,Ord);
                    799:        return G;
                    800: }
                    801:
                    802: def ldim(G,V)
                    803: {
                    804:        O0 = dp_ord(); dp_ord(0);
                    805:        D = length(dp_mbase(map(dp_ptod,G,V)));
                    806:        dp_ord(O0);
                    807:        return D;
                    808: }
                    809:
                    810: def make_mod_subst(GD,V,PV,HC)
                    811: {
                    812:        N = length(V);
                    813:        PN = length(PV);
                    814:        G = GD[0]; D = GD[1];
                    815:        for ( I = 0; ; I = (I+1)%100 ) {
                    816:                Mod = lprime(I);
                    817:                S = [];
                    818:                for ( J = PN-1; J >= 0; J-- )
                    819:                        S = append([PV[J],random()%Mod],S);
                    820:                for ( T = HC; T != []; T = cdr(T) )
                    821:                        if ( !(subst(car(T),S)%Mod) ) break;
                    822:                if ( T != [] ) continue;
                    823:                for ( J = 0; J < N; J++ ) {
                    824:                        M = subst(D[J],S);
                    825:                        F = modsqfr(M,Mod);
                    826:                        if ( length(F) != 2 || F[1][1] != 1 ) break;
                    827:                }
                    828:                if ( J < N ) continue;
                    829:                G0 = map(subst,G,S);
                    830:                return [G0,Mod];
                    831:        }
                    832: }
                    833:
                    834: def rsgn()
                    835: {
                    836:        return random()%2 ? 1 : -1;
                    837: }
                    838:
                    839: def find_npos(GD,V,PV,Mod)
                    840: {
                    841:        N = length(V); PN = length(PV);
                    842:        G = GD[0]; D = GD[1]; LD = D[N];
                    843:        Ord0 = dp_ord(); dp_ord(0);
                    844:        HC = map(dp_hc,map(dp_ptod,G,V));
                    845:        dp_ord(Ord0);
                    846:        if ( !Mod ) {
                    847:                W = append(V,PV);
                    848:                G1 = nd_gr_trace(G,W,1,GBCheck,[[0,N],[0,PN]]);
                    849:                L = make_mod_subst([G1,D],V,PV,HC);
                    850:                return find_npos([L[0],D],V,[],L[1]);
                    851:        }
                    852:        N = length(V);
                    853:        NV = ttttt;
                    854:        for ( B = 2; ; B++ ) {
                    855:                for ( J = N-2; J >= 0; J-- ) {
                    856:                        for ( U = 0, K = J; K < N; K++ )
                    857:                                U += rsgn()*((random()%B+1))*V[K];
                    858:                        M = minipolym(G,V,0,U,NV,Mod);
                    859:                        if ( deg(M,NV) == LD ) return U;
                    860:                }
                    861:        }
                    862: }
                    863:
                    864: def gen_minipoly(G,V,PV,Ord,VI,Mod)
                    865: {
                    866:        if ( PV == [] ) {
                    867:                NV = ttttt;
                    868:                if ( Mod )
                    869:                        M = minipolym(G,V,Ord,VI,NV,Mod);
                    870:                else
                    871:                        M = minipoly(G,V,Ord,VI,NV);
                    872:                return subst(M,NV,VI);
                    873:        }
                    874:        W = setminus(V,[VI]);
                    875:        PV1 = cons(VI,PV);
                    876: #if 0
                    877:        while ( 1 ) {
                    878:                V1 = append(W,PV1);
                    879:                if ( Mod )
                    880:                        G = nd_gr(G,V1,Mod,[[0,1],[0,length(V1)-1]]|nora=1);
                    881:                else
                    882:                        G = nd_gr_trace(G,V1,1,GBCheck,[[0,1],[0,length(V1)-1]]|nora=1);
                    883:                if ( W == [] ) break;
                    884:                else {
                    885:                        W = cdr(W);
                    886:                        G = elimination(G,cdr(V1));
                    887:                }
                    888:        }
                    889: #elif 1
                    890:        if ( Mod ) {
                    891:                G = nd_gr(G,V1,Mod,[[0,length(W)],[0,length(PV1)]]|nora=1);
                    892:                G = elimination(G,PV1);
                    893:        } else {
                    894:                PV2 = setminus(PV1,[PV1[length(PV1)-1]]);
                    895:                V2 = append(W,PV2);
                    896:                G = nd_gr_trace(G,V2,1,GBCheck,[[0,length(W)],[0,length(PV2)]]|nora=1);
                    897:                G = elimination(G,PV1);
                    898:        }
                    899: #else
                    900:        V1 = append(W,PV1);
                    901:        if ( Mod )
                    902:                G = nd_gr(G,V1,Mod,[[0,length(W)],[0,length(PV1)]]|nora=1);
                    903:        else
                    904:                G = nd_gr_trace(G,V1,1,GBCheck,[[0,length(W)],[0,length(PV1)]]|nora=1);
                    905:        G = elimination(G,PV1);
                    906: #endif
                    907:        if ( Mod )
                    908:                G = nd_gr(G,PV1,Mod,[[0,1],[0,length(PV)]]|nora=1);
                    909:        else
                    910:                G = nd_gr_trace(G,PV1,1,GBCheck,[[0,1],[0,length(PV)]]|nora=1);
                    911:        for ( M = car(G), T = cdr(G); T != []; T = cdr(T) )
                    912:                if ( deg(car(T),VI) < deg(M,VI) ) M = car(T);
                    913:        return M;
                    914: }
                    915:
                    916: def indepset(V,H)
                    917: {
                    918:        if ( H == [] ) return V;
                    919:        N = -1;
                    920:        for ( T = V; T != []; T = cdr(T) ) {
                    921:                VI = car(T);
                    922:                HI = [];
                    923:                for ( S = H; S != []; S = cdr(S) )
                    924:                        if ( !tdiv(car(S),VI) ) HI = cons(car(S),HI);
                    925:                RI = indepset(setminus(V,[VI]),HI);
                    926:                if ( length(RI) > N ) {
                    927:                        R = RI; N = length(RI);
                    928:                }
                    929:        }
                    930:        return R;
                    931: }
                    932:
                    933: def maxindep(B,V,O)
                    934: {
                    935:        G = nd_gr_trace(B,V,1,GBCheck,O);
                    936:        Old = dp_ord();
                    937:        dp_ord(O);
                    938:        H = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
                    939:        H = dp_mono_raddec(H,V);
                    940:        N = length(V);
                    941:        Dep = [];
                    942:        for ( T = H, Len = N+1; T != []; T = cdr(T) ) {
                    943:                M = length(car(T));
                    944:                if ( M < Len ) {
                    945:                        Dep = [car(T)];
                    946:                        Len = M;
                    947:                } else if ( M == Len )
                    948:                        Dep = cons(car(T),Dep);
                    949:        }
                    950:        R = setminus(V,Dep[0]);
                    951:        dp_ord(Old);
                    952:        return R;
                    953: }
                    954:
                    955: /* ideal operations */
                    956: def contraction(G,V)
                    957: {
                    958:        C = [];
                    959:        for ( T = G; T != []; T = cdr(T) ) {
                    960:                C1 = dp_hc(dp_ptod(car(T),V));
                    961:                S = fctr(C1);
                    962:                for ( S = cdr(S); S != []; S = cdr(S) )
                    963:                        if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
                    964:        }
                    965:        W = vars(G);
                    966:        PV = setminus(W,V);
                    967:        W = append(V,PV);
                    968:        NV = ttttt;
                    969:        for ( T = C, S = 1; T != []; T = cdr(T) )
                    970:                S *= car(T);
                    971:        G = saturation([G,NV],S,W);
                    972:        return G;
                    973: }
                    974:
                    975: def ideal_list_intersection(L,V,Ord)
                    976: {
                    977:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    978:        N = length(L);
                    979:        if ( N == 0 ) return [1];
                    980:        if ( N == 1 ) return fast_gb(L[0],V,Mod,Ord);
                    981:        N2 = idiv(N,2);
                    982:        for ( L1 = [], I = 0; I < N2; I++ ) L1 = cons(L[I],L1);
                    983:        for ( L2 = []; I < N; I++ ) L2 = cons(L[I],L2);
                    984:        I1 = ideal_list_intersection(L1,V,Ord|mod=Mod);
                    985:        I2 = ideal_list_intersection(L2,V,Ord|mod=Mod);
                    986:        return ideal_intersection(I1,I2,V,Ord|mod=Mod,
                    987:                gbblock=[[0,length(I1)],[length(I1),length(I2)]]);
                    988: }
                    989:
                    990: def ideal_intersection(A,B,V,Ord)
                    991: {
                    992:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    993:        if ( type(Block=getopt(gbblock)) == -1 ) Block = 0;
                    994:        T = ttttt;
                    995:        if ( Mod )
                    996:                G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                    997:                        cons(T,V),Mod,[[0,1],[Ord,length(V)]]);
                    998:        else
                    999:        if ( Procs ) {
                   1000:                Arg0 = ["nd_gr",
                   1001:                        append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1002:                        cons(T,V),0,[[0,1],[Ord,length(V)]]];
                   1003:                Arg1 = ["nd_gr_trace",
                   1004:                        append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1005:                        cons(T,V),1,GBCheck,[[0,1],[Ord,length(V)]]];
                   1006:                G = competitive_exec(Procs,Arg0,Arg1);
                   1007:        } else {
                   1008:                if ( Block )
                   1009:                        G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1010:                                cons(T,V),0,[[0,1],[Ord,length(V)]]|gbblock=Block);
                   1011:                else
                   1012:                        G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1013:                                cons(T,V),0,[[0,1],[Ord,length(V)]]);
                   1014:        }
                   1015:        G0 = elimination(G,V);
                   1016:        return G0;
                   1017: }
                   1018:
                   1019: /* returns GB if F notin rad(G) */
                   1020:
                   1021: def radical_membership(F,G,V) {
                   1022:        F = p_nf(F,G,V,0);
                   1023:        if ( !F ) return 0;
                   1024:        NV = ttttt;
                   1025:        T = nd_gr_trace(cons(NV*F-1,G),cons(NV,V),1,GBCheck,0);
                   1026:        if ( type(car(T)) != 1 ) return [T,NV];
                   1027:        else return 0;
                   1028: }
                   1029:
                   1030: def quick_radical_membership(F,G,V) {
                   1031:        F = p_nf(F,G,V,0);
                   1032:        if ( !F ) return 1;
                   1033:        NV = ttttt;
                   1034:        T = nd_f4(cons(NV*F-1,G),cons(NV,V),lprime(0),0);
                   1035:        if ( type(car(T)) != 1 ) return 0;
                   1036:        else return 1;
                   1037: }
                   1038:
                   1039: def modular_radical_membership(F,G,V) {
                   1040:        F = p_nf(F,G,V,0);
                   1041:        if ( !F ) return 0;
                   1042:        NV = ttttt;
                   1043:        for ( J = 0; ; J++ ) {
                   1044:                Mod = lprime(J);
                   1045:                H = map(dp_hc,map(dp_ptod,G,V));
                   1046:                for ( ; H != []; H = cdr(H) ) if ( !(car(H)%Mod) ) break;
                   1047:                if ( H != [] ) continue;
                   1048:
                   1049:                T = nd_f4(cons(NV*F-1,G),cons(NV,V),Mod,0);
                   1050:                if ( type(car(T)) == 1 ) {
                   1051:                        I = radical_membership_rep(F,G,V,-1,0,Mod);
                   1052:                        I1 = radical_membership_rep(F,G,V,I,0,0);
                   1053:                        if ( I1 > 0 ) return 0;
                   1054:                }
                   1055:                return radical_membership(F,G,V);
                   1056:        }
                   1057: }
                   1058:
                   1059: def radical_membership_rep(F,G,V,Max,Ord,Mod) {
                   1060:        Ft = F;
                   1061:        I = 1;
                   1062:        while ( Max < 0 || I <= Max ) {
                   1063:                Ft = nd_nf(Ft,G,V,Ord,Mod);
                   1064:                if ( !Ft ) return I;
                   1065:                Ft *= F;
                   1066:                I++;
                   1067:        }
                   1068:        return -1;
                   1069: }
                   1070:
                   1071: def ideal_product(A,B,V)
                   1072: {
                   1073:        dp_ord(0);
                   1074:        DA = map(dp_ptod,A,V);
                   1075:        DB = map(dp_ptod,B,V);
                   1076:        DegA = map(dp_td,DA);
                   1077:        DegB = map(dp_td,DB);
                   1078:        for ( PA = [], T = A, DT = DegA; T != []; T = cdr(T), DT = cdr(DT) )
                   1079:                PA = cons([car(T),car(DT)],PA);
                   1080:        PA = reverse(PA);
                   1081:        for ( PB = [], T = B, DT = DegB; T != []; T = cdr(T), DT = cdr(DT) )
                   1082:                PB = cons([car(T),car(DT)],PB);
                   1083:        PB = reverse(PB);
                   1084:        R = [];
                   1085:        for ( T = PA; T != []; T = cdr(T) )
                   1086:                for ( S = PB; S != []; S = cdr(S) )
                   1087:                        R = cons([car(T)[0]*car(S)[0],car(T)[1]+car(S)[1]],R);
                   1088:        T = qsort(R,comp_by_second);
                   1089:        T = map(first_element,T);
                   1090:        Len = length(A)>length(B)?length(A):length(B);
                   1091:        Len *= 2;
                   1092:        L = sep_list(T,Len); B0 = L[0]; B1 = L[1];
                   1093:        R = nd_gr_trace(B0,V,0,-1,0);
                   1094:        while ( B1 != [] ) {
                   1095:                print(length(B1));
                   1096:                L = sep_list(B1,Len);
                   1097:                B0 = L[0]; B1 = L[1];
                   1098:                R = nd_gr_trace(append(R,B0),V,0,-1,0|gbblock=[[0,length(R)]],nora=1);
                   1099:        }
                   1100:        return R;
                   1101: }
                   1102:
                   1103: def saturation(GNV,F,V)
                   1104: {
                   1105:        G = GNV[0]; NV = GNV[1];
                   1106:        if ( Procs ) {
                   1107:                Arg0 = ["nd_gr_trace",
                   1108:                cons(NV*F-1,G),cons(NV,V),0,GBCheck,[[0,1],[0,length(V)]]];
                   1109:                Arg1 = ["nd_gr_trace",
                   1110:                cons(NV*F-1,G),cons(NV,V),1,GBCheck,[[0,1],[0,length(V)]]];
                   1111:                G1 = competitive_exec(Procs,Arg0,Arg1);
                   1112:        } else
                   1113:                G1 = nd_gr_trace(cons(NV*F-1,G),cons(NV,V),SatHomo,GBCheck,[[0,1],[0,length(V)]]);
                   1114:        return elimination(G1,V);
                   1115: }
                   1116:
                   1117: def sat(G,F,V)
                   1118: {
                   1119:        NV = ttttt;
                   1120:        if ( Procs ) {
                   1121:                Arg0 = ["nd_gr_trace",
                   1122:                cons(NV*F-1,G),cons(NV,V),0,GBCheck,[[0,1],[0,length(V)]]];
                   1123:                Arg1 = ["nd_gr_trace",
                   1124:                cons(NV*F-1,G),cons(NV,V),1,GBCheck,[[0,1],[0,length(V)]]];
                   1125:                G1 = competitive_exec(Procs,Arg0,Arg1);
                   1126:        } else
                   1127:                G1 = nd_gr_trace(cons(NV*F-1,G),cons(NV,V),SatHomo,GBCheck,[[0,1],[0,length(V)]]);
                   1128:        return elimination(G1,V);
                   1129: }
                   1130:
                   1131: def satind(G,F,V)
                   1132: {
                   1133:        NV = ttttt;
                   1134:        N = length(V);
                   1135:        B = append(G,[NV*F-1]);
                   1136:        V1 = cons(NV,V);
                   1137:        D = nd_gr_trace(B,V1,1,GBCheck,[[0,1],[0,N]]
                   1138:                |nora=1,gentrace=1,gbblock=[[0,length(G)]]);
                   1139:        G1 = D[0];
                   1140:        Len = length(G1);
                   1141:        Deg = compute_deg(B,V1,NV,D);
                   1142:        D1 = 0;
                   1143:        R = [];
                   1144:        M = length(B);
                   1145:        for ( I = 0; I < Len; I++ ) {
                   1146:                if ( !member(NV,vars(G1[I])) ) {
                   1147:                        for ( J = 1; J < M; J++ )
                   1148:                                D1 = MAX(D1,Deg[I][J]);
                   1149:                        R = cons(G1[I],R);
                   1150:                }
                   1151:        }
                   1152:        return [reverse(R),D1];
                   1153: }
                   1154:
                   1155: def sat_ind(G,F,V)
                   1156: {
                   1157:        NV = ttttt;
                   1158:        F = p_nf(F,G,V,0);
                   1159:        for ( I = 0, GI = G; ; I++ ) {
                   1160:                G1 = colon(GI,F,V);
                   1161:                if ( ideal_inclusion(G1,GI,V,0) )  {
                   1162:                        return [GI,I];
                   1163:                }
                   1164:                else GI = G1;
                   1165:        }
                   1166: }
                   1167:
                   1168: def colon(G,F,V)
                   1169: {
                   1170:        F = p_nf(F,G,V,0);
                   1171:        if ( !F ) return [1];
                   1172:        NV = ttttt;
                   1173:        V1 = cons(NV,V);
                   1174:        T = nd_gr_trace(append(vtol(NV*ltov(G)),[(1-NV)*F]),V1,1,GBCheck,
                   1175:                [[0,1],[0,length(V)]]|gbblock=[[0,length(G)]],nora=1);
                   1176:        T = elimination(T,V);
                   1177:        return map(ptozp,map(sdiv,T,F));
                   1178: }
                   1179:
                   1180: def ideal_colon(G,F,V)
                   1181: {
                   1182:        G = nd_gr(G,V,0,0);
                   1183:        L = mapat(colon,1,G,F,V);
                   1184:        return ideal_list_intersection(L,V,0);
                   1185: }
                   1186:
1.2     ! noro     1187: def ideal_sat(G,F,V)
        !          1188: {
        !          1189:        G = nd_gr(G,V,0,0);
        !          1190:        L = mapat(sat,1,G,F,V);
        !          1191:        return ideal_list_intersection(L,V,0);
        !          1192: }
        !          1193:
1.1       noro     1194: def ideal_inclusion(F,G,V,O)
                   1195: {
                   1196:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1197:        if ( Mod ) {
                   1198:                for ( T = F; T != []; T = cdr(T) )
                   1199:                        if ( p_nf_mod(car(T),G,V,O,Mod) ) return 0;
                   1200:        } else {
                   1201:                for ( T = F; T != []; T = cdr(T) )
                   1202:                        if ( p_nf(car(T),G,V,O) ) return 0;
                   1203:        }
                   1204:        return 1;
                   1205: }
                   1206:
                   1207: /* remove redundant components */
                   1208:
                   1209: def qd_simp_comp(QP,V)
                   1210: {
                   1211:        R = ltov(QP);
                   1212:        N = length(R);
                   1213:        for ( I = 0; I < N; I++ ) {
                   1214:                if ( R[I] ) {
                   1215:                        QI = R[I][0]; PI = R[I][1];
                   1216:                        for ( J = I+1; J < N; J++ )
                   1217:                                if ( R[J] && gb_comp(PI,R[J][1]) ) {
                   1218:                                        QI = ideal_intersection(QI,R[J][0],V,0);
                   1219:                                        R[J] = 0;
                   1220:                                }
                   1221:                        R[I] = [QI,PI];
                   1222:                }
                   1223:        }
                   1224:        for ( I = N-1, S = []; I >= 0; I-- )
                   1225:                if ( R[I] ) S = cons(R[I],S);
                   1226:        return S;
                   1227: }
                   1228:
                   1229: def qd_remove_redundant_comp(G,Iso,Emb,V,Ord)
                   1230: {
                   1231:        IsoInt = ideal_list_intersection(map(first_element,Iso),V,Ord);
                   1232:        Emb = qd_simp_comp(Emb,V);
                   1233:        Emb = qsort(Emb);
                   1234:        A = ltov(Emb);
                   1235:        N = length(A);
                   1236:        for ( I = 0; I < N; I++ ) {
                   1237:                if ( !A[I] ) continue;
                   1238:                for ( T = [], J = 0; J < N; J++ )
                   1239:                        if ( J != I && A[J] ) T = cons(A[J][0],T);
                   1240:                Int = ideal_list_intersection(T,V,Ord);
                   1241:                Int = ideal_intersection(IsoInt,Int,V,Ord);
                   1242:                if ( gb_comp(Int,G) ) A[I] = 0;
                   1243:        }
                   1244:        for ( T = [], I = 0; I < N; I++ )
                   1245:                if ( A[I] ) T = cons(A[I],T);
                   1246:        return reverse(T);
                   1247: }
                   1248:
                   1249: def remove_redundant_comp(G,Iso,Emb,V,Ord)
                   1250: {
                   1251:        IsoInt = ideal_list_intersection(Iso,V,Ord);
                   1252:
                   1253:        A = ltov(Emb);
                   1254:        N = length(A);
                   1255:        for ( I = 0; I < N; I++ ) {
                   1256:                if ( !A[I] ) continue;
                   1257:                for ( J = I+1; J < N; J++ )
                   1258:                        if ( A[J] && gb_comp(A[I],A[J]) ) A[J] = 0;
                   1259:        }
                   1260:        for ( I = 0; I < N; I++ ) {
                   1261:                if ( !A[I] ) continue;
                   1262:                for ( T = [], J = 0; J < N; J++ )
                   1263:                        if ( J != I && A[J] ) T = cons(A[J],T);
                   1264:                Int = ideal_list_intersection(cons(IsoInt,T),V,Ord);
                   1265:                if ( gb_comp(Int,G) ) A[I] = 0;
                   1266:        }
                   1267:        for ( T = [], I = 0; I < N; I++ )
                   1268:                if ( A[I] ) T = cons(A[I],T);
                   1269:        return reverse(T);
                   1270: }
                   1271:
                   1272: def remove_redundant_comp_first(G,P,V,Ord)
                   1273: {
                   1274:        A = ltov(P);
                   1275:        N = length(A);
                   1276:        for ( I = 0; I < N; I++ ) {
                   1277:                if ( !A[I] ) continue;
                   1278:                for ( J = I+1; J < N; J++ )
                   1279:                        if ( A[J] && gb_comp(A[I][0],A[J][0]) ) A[J] = 0;
                   1280:        }
                   1281:        for ( I = 0; I < N; I++ ) {
                   1282:                if ( !A[I] ) continue;
                   1283:                for ( T = [], J = 0; J < N; J++ )
                   1284:                        if ( J != I && A[J] ) T = cons(A[J][0],T);
                   1285:                Int = ideal_list_intersection(T,V,Ord);
                   1286:                if ( gb_comp(Int,G) ) A[I] = 0;
                   1287:        }
                   1288:        for ( T = [], I = 0; I < N; I++ )
                   1289:                if ( A[I] ) T = cons(A[I],T);
                   1290:        return reverse(T);
                   1291: }
                   1292:
                   1293: /* polynomial operations */
                   1294:
                   1295: def ppart(F,V,Mod)
                   1296: {
                   1297:        if ( !Mod )
                   1298:                G = nd_gr([F],[V],0,0);
                   1299:        else
                   1300:                G = dp_gr_mod_main([F],[V],0,Mod,0);
                   1301:        return G[0];
                   1302: }
                   1303:
                   1304:
                   1305: def sq(F)
                   1306: {
                   1307:        if ( !F ) return 0;
                   1308:        A = cdr(fctr(F));
                   1309:        for ( R = 1; A != []; A = cdr(A) )
                   1310:                R *= car(car(A));
                   1311:        return R;
                   1312: }
                   1313:
                   1314: def lcfactor(G,V,O)
                   1315: {
                   1316:        O0 = dp_ord(); dp_ord(O);
                   1317:        C = [];
                   1318:        for ( T = G; T != []; T = cdr(T) ) {
                   1319:                C1 = dp_hc(dp_ptod(car(T),V));
                   1320:                S = fctr(C1);
                   1321:                for ( S = cdr(S); S != []; S = cdr(S) )
                   1322:                        if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
                   1323:        }
                   1324:        dp_ord(O0);
                   1325:        return C;
                   1326: }
                   1327:
                   1328: /* Ti = [D,I,M,C] */
                   1329:
                   1330: def compute_deg0(Ti,P,V,TV)
                   1331: {
                   1332:        N = length(P[0]);
                   1333:        Num = vector(N);
                   1334:        for ( I = 0; I < N; I++ ) Num[I] = -1;
                   1335:        for ( ; Ti != []; Ti = cdr(Ti) ) {
                   1336:                Sj = car(Ti);
                   1337:                Dj = Sj[0];
                   1338:                Ij =Sj[1];
                   1339:                Mj = deg(type(Sj[2])==9?dp_dtop(Sj[2],V):Sj[2],TV);
                   1340:                Pj = P[Ij];
                   1341:                if ( Dj )
                   1342:                        for ( I = 0; I < N; I++ )
                   1343:                                if ( Pj[I] >= 0 ) {
                   1344:                                        T = Mj+Pj[I];
                   1345:                                        Num[I] = MAX(Num[I],T);
                   1346:                                }
                   1347:        }
                   1348:        return Num;
                   1349: }
                   1350:
                   1351: def compute_deg(B,V,TV,Data)
                   1352: {
                   1353:        GB = Data[0];
                   1354:        Homo = Data[1];
                   1355:        Trace = Data[2];
                   1356:        IntRed = Data[3];
                   1357:        Ind = Data[4];
                   1358:        DB = map(dp_ptod,B,V);
                   1359:        if ( Homo ) {
                   1360:                DB = map(dp_homo,DB);
                   1361:                V0 = append(V,[hhh]);
                   1362:        } else
                   1363:                V0 = V;
                   1364:        Perm = Trace[0]; Trace = cdr(Trace);
                   1365:        for ( I = length(Perm)-1, T = Trace; T != []; T = cdr(T) )
                   1366:                if ( (J=car(T)[0]) > I ) I = J;
                   1367:        N = I+1;
                   1368:        N0 = length(B);
                   1369:        P = vector(N);
                   1370:        for ( T = Perm, I = 0; T != []; T = cdr(T), I++ ) {
                   1371:                Pi = car(T);
                   1372:                C = vector(N0);
                   1373:                for ( J = 0; J < N0; J++ ) C[J] = -1;
                   1374:                C[Pi[1]] = 0;
                   1375:                P[Pi[0]] = C;
                   1376:        }
                   1377:        for ( T = Trace; T != []; T = cdr(T) ) {
                   1378:                Ti = car(T); P[Ti[0]] = compute_deg0(Ti[1],P,V0,TV);
                   1379:        }
                   1380:        M = length(Ind);
                   1381:        for ( T = IntRed; T != []; T = cdr(T) ) {
                   1382:                Ti = car(T); P[Ti[0]] = compute_deg0(Ti[1],P,V,TV);
                   1383:        }
                   1384:        R = [];
                   1385:        for ( J = 0; J < M; J++ ) {
                   1386:                U = P[Ind[J]];
                   1387:                R = cons(U,R);
                   1388:        }
                   1389:        return reverse(R);
                   1390: }
                   1391:
                   1392: /* set theoretic functions */
                   1393:
                   1394: def member(A,S)
                   1395: {
                   1396:        for ( ; S != []; S = cdr(S) )
                   1397:                if ( car(S) == A ) return 1;
                   1398:        return 0;
                   1399: }
                   1400:
                   1401: def elimination(G,V) {
                   1402:        for ( R = [], T = G; T != []; T = cdr(T) )
                   1403:                if ( setminus(vars(car(T)),V) == [] ) R =cons(car(T),R);
                   1404:        return R;
                   1405: }
                   1406:
                   1407: def setintersection(A,B)
                   1408: {
                   1409:        for ( L = []; A != []; A = cdr(A) )
                   1410:                if ( member(car(A),B) )
                   1411:                        L = cons(car(A),L);
                   1412:        return L;
                   1413: }
                   1414:
                   1415: def setminus(A,B) {
                   1416:        for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
                   1417:                for ( S = B, M = car(T); S != []; S = cdr(S) )
                   1418:                        if ( car(S) == M ) break;
                   1419:                if ( S == [] ) R = cons(M,R);
                   1420:        }
                   1421:        return R;
                   1422: }
                   1423:
                   1424: def sep_list(L,N)
                   1425: {
                   1426:        if ( length(L) <= N ) return [L,[]];
                   1427:        R = [];
                   1428:        for ( T = L, I = 0; I < N; I++, T = cdr(T) )
                   1429:                R = cons(car(T),R);
                   1430:        return [reverse(R),T];
                   1431: }
                   1432:
                   1433: def first_element(L)
                   1434: {
                   1435:        return L[0];
                   1436: }
                   1437:
                   1438: def comp_tdeg(A,B)
                   1439: {
                   1440:        DA = tdeg(A);
                   1441:        DB = tdeg(B);
                   1442:        if ( DA > DB ) return 1;
                   1443:        else if ( DA < DB ) return -1;
                   1444:        else return 0;
                   1445: }
                   1446:
                   1447: def tdeg(P)
                   1448: {
                   1449:        dp_ord(0);
                   1450:        return dp_td(dp_ptod(P,vars(P)));
                   1451: }
                   1452:
                   1453: def comp_by_ord(A,B)
                   1454: {
                   1455:        if ( dp_ht(A) > dp_ht(B) ) return 1;
                   1456:        else if ( dp_ht(A) < dp_ht(B) ) return -1;
                   1457:        else return 0;
                   1458: }
                   1459:
                   1460: def comp_by_second(A,B)
                   1461: {
                   1462:        if ( A[1] > B[1] ) return 1;
                   1463:        else if ( A[1] < B[1] ) return -1;
                   1464:        else return 0;
                   1465: }
                   1466: endmodule$
                   1467: end$

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