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

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

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