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

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

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