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

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

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