[BACK]Return to noro_pd.rr CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / packages / src

Annotation of OpenXM/src/asir-contrib/packages/src/noro_pd.rr, Revision 1.5

1.5     ! noro        1: /* $OpenXM: OpenXM/src/asir-contrib/packages/src/noro_pd.rr,v 1.4 2017/02/22 06:20:00 noro Exp $ */
1.1       noro        2: import("gr")$
                      3: module noro_pd$
1.5     ! noro        4: static GBF4,GBTrace,GBCheck,Chrem,EProcs,Procs,SatHomo,GBRat,SuccSat,RepColon,IntGB,UsualInt$
1.1       noro        5:
                      6: localf radical_membership_sat$
                      7: localf witness$
                      8: localf get_lc,tomonic,aa,ideal_intersection_m,redbase$
                      9: localf para_exec,nd_gr_rat,competitive_exec,call_func,call_func_serial$
                     10: localf call_ideal_list_intersection$
                     11: localf call_colon,call_prime_dec$
                     12: localf prime_dec2, prime_dec_main2$
                     13: localf first_second$
                     14: localf third$
                     15: localf locsat,iso_comp_para,extract_qj,colon_prime_dec,extract_comp$
                     16: localf separator$
1.3       noro       17: localf member,subset,mingen,compute_gbsyz,redcoef,recompute_trace,dtop,topnum$
1.1       noro       18: localf prepost$
                     19: localf monodec0,monodec,prod$
                     20: localf extract_qd,primary_check$
                     21: localf second$
                     22: localf gbrat,succsat,repcolon,comp_third_tdeg,comp_tord$
                     23: localf power$
                     24:
1.2       noro       25: localf syci_dec$
1.1       noro       26:
                     27: localf find_si0,find_si1,find_si2$
                     28: localf find_ssi0,find_ssi1,find_ssi2$
                     29:
                     30: localf init_pprocs, init_eprocs, init_procs, kill_procs$
                     31:
1.2       noro       32: localf iso_comp$
1.1       noro       33:
                     34: localf prime_dec, prime_dec_main, lex_predec1, zprimedec, zprimadec$
                     35: localf complete_qdecomp, partial_qdecomp, partial_qdecomp0, complete_decomp$
                     36: localf partial_decomp, partial_decomp0, zprimacomp, zprimecomp$
                     37: localf fast_gb, incremental_gb, elim_gb, ldim, make_mod_subst$
                     38: localf rsgn, find_npos, gen_minipoly, indepset$
                     39: localf maxindep, maxindep2, contraction, contraction_m, ideal_list_intersection, ideal_intersection$
                     40: localf radical_membership, modular_radical_membership$
                     41: localf radical_membership_rep, ideal_product, saturation$
1.3       noro       42: localf sat, sat_ind, sat_ind_var, colon, isat$
1.1       noro       43: localf ideal_colon, ideal_sat, ideal_inclusion, qd_simp_comp, qd_remove_redundant_comp$
                     44: localf pd_simp_comp, remove_identical_comp$
                     45: localf pd_remove_redundant_comp, ppart, sq, gen_fctr, gen_nf, gen_gb_comp$
                     46: localf gen_mptop, lcfactor, compute_deg0, compute_deg, member$
                     47: localf elimination, setintersection, setminus, sep_list$
                     48: localf first, comp_tdeg, comp_tdeg_first, tdeg, comp_by_ord, comp_by_second$
1.5     ! noro       49: localf gbcheck,chrem,sathomo,qd_check,qdb_check,intgb,gbtrace,gbf4,usualint$
        !            50: localf zero_dim$
1.1       noro       51:
                     52: SatHomo=0$
                     53: GBCheck=1$
                     54: GBRat=0$
                     55: SuccSat=0$
                     56: RepColon=0$
1.3       noro       57: IntGB=0$
1.5     ! noro       58: UsualInt=0$
1.3       noro       59: GBTrace=0$
                     60: GBF4=0$
1.1       noro       61:
                     62: #define MAX(a,b) ((a)>(b)?(a):(b))
                     63: #define ACCUM_TIME(C,R) {T1 = time(); C += (T1[0]-T0[0])+(T1[1]-T0[1]); R += (T1[3]-T0[3]); }
                     64:
1.5     ! noro       65: def usualint(A)
        !            66: {
        !            67:        UsualInt = A;
        !            68: }
        !            69:
1.3       noro       70: def gbf4(A)
                     71: {
                     72:        GBF4 = A;
                     73: }
                     74:
                     75: def gbtrace(A)
                     76: {
                     77:        GBTrace = A;
                     78: }
                     79:
1.1       noro       80: def gbrat(A)
                     81: {
                     82:        GBRat = A;
                     83: }
                     84:
                     85: def succsat(A)
                     86: {
                     87:        SuccSat = A;
                     88: }
                     89:
                     90: def repcolon(A)
                     91: {
                     92:        RepColon = A;
                     93: }
                     94:
                     95: def gbcheck(A)
                     96: {
                     97:        if ( A ) GBCheck = 1;
                     98:        else GBCheck = -1;
                     99: }
                    100:
1.3       noro      101: def chrem(A)
1.1       noro      102: {
1.3       noro      103:        if ( A ) Chrem = 1;
                    104:        else Chrem = 0;
1.1       noro      105: }
                    106:
                    107: def sathomo(A)
                    108: {
                    109:        if ( A ) SatHomo = 1;
                    110:        else SatHomo = 0;
                    111: }
                    112:
1.3       noro      113: def intgb(A)
                    114: {
                    115:        if ( A ) IntGB = 1;
                    116:        else IntGB = 0;
                    117: }
                    118:
                    119:
1.1       noro      120: def init_eprocs()
                    121: {
                    122:        if ( type(NoX=getopt(nox)) == -1 ) NoX = 0;
                    123:        if ( !EProcs ) {
                    124:                if ( NoX ) {
                    125:                        P0 = ox_launch_nox();
                    126:                        P1 = ox_launch_nox();
                    127:                } else {
                    128:                        P0 = ox_launch();
                    129:                        P1 = ox_launch();
                    130:                }
                    131:                EProcs = [P0,P1];
                    132:        }
                    133: }
                    134:
                    135: def init_pprocs(N)
                    136: {
                    137:        if ( type(NoX=getopt(nox)) == -1 ) NoX = 0;
                    138:        for ( R = [], I = 0; I < N; I++ ) {
                    139:                P = NoX ? ox_launch_nox() : ox_launch();
                    140:                R = cons(P,R);
                    141:        }
                    142:        return reverse(R);
                    143: }
                    144:
                    145: def init_procs()
                    146: {
                    147:        if ( type(NoX=getopt(nox)) == -1 ) NoX = 0;
                    148:        if ( !Procs ) {
                    149:                if ( NoX ) {
                    150:                        P0 = ox_launch_nox();
                    151:                        P1 = ox_launch_nox();
                    152:                } else {
                    153:                        P0 = ox_launch();
                    154:                        P1 = ox_launch();
                    155:                }
                    156:                Procs = [P0,P1];
                    157:        }
                    158: }
                    159:
                    160: def kill_procs()
                    161: {
                    162:        if ( Procs ) {
                    163:                ox_shutdown(Procs[0]);
                    164:                ox_shutdown(Procs[1]);
                    165:                Procs = 0;
                    166:        }
                    167:        if ( EProcs ) {
                    168:                ox_shutdown(EProcs[0]);
                    169:                ox_shutdown(EProcs[1]);
                    170:                EProcs = 0;
                    171:        }
                    172: }
                    173:
                    174: def qd_check(B,V,QD)
                    175: {
                    176:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
1.3       noro      177:        G = fast_gb(B,V,Mod,0);
1.1       noro      178:        Iso = ideal_list_intersection(map(first,QD[0]),V,0|mod=Mod);
                    179:        Emb = ideal_list_intersection(map(first,QD[1]),V,0|mod=Mod);
                    180:        GG = ideal_intersection(Iso,Emb,V,0|mod=Mod);
                    181:        return gen_gb_comp(G,GG,Mod);
                    182: }
                    183:
                    184: def primary_check(B,V)
                    185: {
                    186:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
1.3       noro      187:        G = fast_gb(B,V,Mod,0);
1.1       noro      188:        PL = prime_dec(G,V|indep=1,mod=Mod);
                    189:        if ( length(PL) > 1 ) return 0;
                    190:        P = PL[0][0]; Y = PL[0][1];
                    191:        Z = setminus(V,Y);
                    192:        H = elim_gb(G,Z,Y,Mod,[[0,length(Z)],[0,length(Y)]]);
                    193:        H = contraction(H,Z|mod=Mod);
1.3       noro      194:        H = fast_gb(H,V,Mod,0);
                    195:        if ( gen_gb_comp(G,H,Mod) ) return fast_gb(P,V,Mod,0);
1.1       noro      196:        else return 0;
                    197: }
                    198:
                    199: def qdb_check(B,V,QD)
                    200: {
                    201:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
1.3       noro      202:        G = fast_gb(B,V,Mod,0);
1.1       noro      203:        N = length(QD);
                    204:        for ( I = 0, Q = [1]; I < N; I++ )
                    205:                for ( J = 0, QL = map(first,QD[I]), L = length(QL);
                    206:                        J < L; J++ )
                    207:                        Q = ideal_intersection_m(Q,QL[J],V,0|mod=Mod);
1.3       noro      208:        Q = fast_gb(Q,V,0,0);
1.1       noro      209:        if ( !gen_gb_comp(G,Q,Mod) )
                    210:                return 0;
                    211:        for ( I = 0; I < N; I++ ) {
                    212:                T = QD[I];
                    213:                M = length(T);
                    214:                for ( J = 0; J < M; J++ ) {
                    215:                        P = primary_check(T[J][0],V|mod=Mod);
                    216:                        if ( !P ) return 0;
1.3       noro      217:                        PP = fast_gb(T[J][1],V,Mod,0);
1.1       noro      218:                        if ( !gen_gb_comp(P,PP,Mod) ) return 0;
                    219:                }
                    220:        }
                    221:        return 1;
                    222: }
                    223:
                    224: def extract_qd(QD,V,Ind)
                    225: {
                    226:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    227:        N = length(Ind);
                    228:        for ( I = 0, Q = [1]; I < N; I++ )
                    229:                for ( J = 0, QL = map(first,QD[Ind[I]]), L = length(QL);
                    230:                        J < L; J++ )
                    231:                        Q = ideal_intersection(Q,QL[J],V,0|mod=Mod);
                    232:        return Q;
                    233: }
                    234:
                    235: static Tint2, RTint2$
                    236:
                    237: def syci_dec(B,V)
                    238: {
                    239:        if ( type(SI=getopt(si)) == -1 ) SI = 1;
                    240:        if ( SI<0 || SI>2 ) error("sycb_assdec : si should be 0,1,2");
                    241:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    242:        if ( type(Lexdec=getopt(lexdec)) == -1 ) Lexdec = 0;
                    243:        if ( type(Time=getopt(time)) == -1 ) Time = 0;
                    244:        if ( type(Iso=getopt(iso)) == -1 ) Iso = 0;
                    245:        if ( type(Ass=getopt(ass)) == -1 ) Ass = 0;
                    246:        if ( type(Colon=getopt(colon)) == -1 ) Colon = 0;
                    247:        if ( type(Para=getopt(para)) == -1 ) Para = 0;
1.4       noro      248:     /* global options */
1.3       noro      249:        if ( type(GBF4=getopt(f4)) == -1 ) GBF4 = 0;
                    250:        if ( type(GBTrace=getopt(trace)) == -1 ) GBTrace = 0;
1.4       noro      251:        if ( type(IntGB=getopt(intgb)) == -1 ) IntGB = 0;
1.3       noro      252:
1.1       noro      253:        Ord = 0;
                    254:        Tiso = Tint = Tpd = Text = Tint2 = 0;
                    255:        RTiso = RTint = RTpd = RText = RTint2 = 0;
                    256: T00 = time();
1.3       noro      257:        G = fast_gb(B,V,Mod,Ord);
1.1       noro      258:        IntQ = [1]; QL = RL = []; First = 1;
                    259:        for ( Level = 0; ; Level++ ) {
                    260: T0 = time();
                    261:                if ( !Level ) {
                    262:                        PtR = prime_dec(G,V|indep=1,lexdec=Lexdec,mod=Mod,radical=1);
                    263: ACCUM_TIME(Tfpd,RTfpd)
                    264:                        Pt = PtR[0]; IntPt = PtR[1]; Rad = IntPt;
                    265:                        if ( gen_gb_comp(G,Rad,Mod) ) {
                    266:                                /* Gt is radical and Gt = cap Pt */
                    267:                                for ( T = Pt, Qt = []; T != []; T = cdr(T) )
                    268:                                        Qt = cons([car(T)[0],car(T)[0],car(T)[1]],Qt);
                    269:                                return [reverse(Qt)];
                    270:                        }
                    271:                } else
                    272:                        Pt = colon_prime_dec(G,IntQ,V|lexdec=Lexdec,mod=Mod,para=Para);
                    273: ACCUM_TIME(Tpd,RTpd)
                    274: T0 = time();
                    275:                Rt = iso_comp(G,Pt,V,Ord|mod=Mod,iso=Iso,para=Para,intq=IntQ);
                    276: ACCUM_TIME(Tiso,RTiso)
                    277:                if ( !Level ) {
                    278:                        if ( Iso == 3 ) {
                    279:                                NI = length(Rt);
                    280:                                Q = IntQ;
                    281: T0 = time();
                    282:                                if ( Para ) {
                    283:                                        for ( J = 0, Task = []; J < NI; J++ ) {
                    284:                                                T = ["noro_pd.extract_qj",Q,V,Rt[J],Rad,Mod,SI,Colon,-1];
                    285:                                                Task = cons(T,Task);
                    286:                                        }
                    287:                                        Task = reverse(Task);
                    288:                                        print("comps:",2); print(length(Task),2);
                    289:                                        Rt = para_exec(Para,Task);
                    290:                                } else {
                    291:                                        for ( J = 0, T = []; J < NI; J++ ) {
                    292:                                                TJ = extract_qj(Q,V,Rt[J],Rad,Mod,SI,Colon,-1);
                    293:                                                T = cons(TJ,T);
                    294:                                        }
                    295:                                        Rt = reverse(T);
                    296:                                }
                    297: ACCUM_TIME(Text,RText)
                    298:                        }
                    299:                        print("");
                    300: T0 = time();
                    301:                        Int = Rad;
                    302:                        for ( T = Rt; T != []; T = cdr(T) )
                    303:                                if ( !gb_comp(car(T)[0],car(T)[1]) )
                    304:                                        Int = ideal_intersection_m(Int,car(T)[0],V,Ord|mod=Mod);
1.3       noro      305:                        IntQ = fast_gb(Int,V,Mod,Ord);
1.1       noro      306: ACCUM_TIME(Tint,RTint)
                    307:                        RL = append(RL,[Rt]);
                    308:                } else if ( Iso != 3 ) {
                    309: T0 = time();
                    310:                        IntQ = ideal_list_intersection(map(first,Rt),V,Ord|mod=Mod,isgb=1);
                    311:                        RL = append(RL,[Rt]);
                    312: ACCUM_TIME(Tint,RTint)
                    313:                } else {
                    314:                        NI = length(Rt);
                    315:                        Q = IntQ;
                    316:                        if ( Para ) {
                    317:                                for ( J = 0, Task = []; J < NI; J++ ) {
                    318:                                        T = ["noro_pd.extract_qj",Q,V,Rt[J],Rad,Mod,SI,Colon,-1];
                    319:                                        Task = cons(T,Task);
                    320:                                }
                    321:                                Task = reverse(Task);
                    322:                                print("comps:",2); print(length(Task),2);
                    323: T0 = time();
                    324:                                R = para_exec(Para,Task);
                    325: ACCUM_TIME(Text,RText)
                    326:                                print("");
                    327: T0 = time();
                    328:                                IntQ = ideal_list_intersection(cons(IntQ,map(first,R)),V,Ord|mod=Mod);
                    329: ACCUM_TIME(Tint,RTint)
                    330:                                RL = append(RL,[R]);
                    331:                        } else {
                    332:                                for ( J = 0, T = []; J < NI; J++ ) {
                    333: T0 = time();
                    334:                                        TJ = extract_qj(Q,V,Rt[J],Rad,Mod,SI,Colon,-1);
                    335: ACCUM_TIME(Text,RText)
                    336:                                        T = cons(TJ,T);
                    337: T0 = time();
                    338:                                        IntQ = ideal_intersection_m(IntQ,TJ[0],V,Ord|mod=Mod);
                    339: ACCUM_TIME(Tint,RTint)
                    340:                                }
                    341:                                print("");
                    342: T0 = time();
1.3       noro      343:                                IntQ = fast_gb(IntQ,V,Mod,Ord);
1.1       noro      344: ACCUM_TIME(Tint,RTint)
                    345:                                T = reverse(T); RL = append(RL,[T]);
                    346:                        }
                    347:                }
                    348:                QL = append(QL,[IntQ]);
                    349:                if ( gen_gb_comp(IntQ,G,Mod) ) break;
                    350:        }
                    351: T0 = time();
                    352:        if ( Iso != 3 && !Ass )
                    353:                RL = extract_comp(QL,RL,V,Rad|mod=Mod,para=Para,si=SI,colon=Colon,ass=Ass);
                    354: ACCUM_TIME(Text,RText)
                    355:        if ( Time ) {
                    356:                T1 = time();
                    357:                Tall = T1[0]-T00[0]+T1[1]-T00[1]; RTall += T1[3]-T00[3];
                    358:                Tass = Tall-Text; RTass = RTall-RText;
                    359:                print(["total",Tall,"ass",Tass,"pd",Tpd,"(fpd)",Tfpd,"iso",Tiso,"int",Tint,"ext",Text]);
                    360:                print(["elapsed",RTall,"ass",RTass,"pd",RTpd,"(fpd)",RTfpd,"iso",RTiso,"int",RTint,"ext",RText]);
                    361:        }
                    362:        return  RL;
                    363: }
                    364:
                    365: def extract_comp(QL,RL,V,Rad) {
                    366:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    367:        if ( type(Para=getopt(para)) == -1 ) Para = 0;
                    368:        if ( type(Colon=getopt(colon)) == -1 ) Colon = 0;
                    369:        if ( type(SI=getopt(si)) == -1 ) SI = 1;
                    370:        if ( type(Ass=getopt(ass)) == -1 ) Ass = 0;
                    371:
                    372:        L = length(QL);
                    373:        if ( Para ) {
                    374:                for ( Task = [], I = 1; I < L; I++ ) {
                    375:                        QI = QL[I-1]; RI = RL[I]; NI = length(RI);
                    376:                        for ( J = 0, TI = []; J < NI; J++ ) {
                    377:                                T = ["noro_pd.extract_qj",QI,V,RI[J],Rad,Mod,SI,Colon,I];
                    378:                                Task = cons(T,Task);
                    379:                        }
                    380:                }
                    381:                Task = reverse(Task);
                    382:                print("comps:",2); print(length(Task),2); print("");
                    383:                R = para_exec(Para,Task);
                    384:                S = vector(L);
                    385:                for ( I = 1; I < L; I++ ) S[I] = [];
                    386:                S[0] = RL[0];
                    387:                for ( T = R; T != []; T = cdr(T) ) {
                    388:                        U = car(T); Level = U[0]; Body = U[1];
                    389:                        S[Level] = cons(Body,S[Level]);
                    390:                }
                    391:                return vtol(S);
                    392:        } else {
                    393:                TL = [RL[0]];
                    394:                for ( I = 1; I < L; I++ ) {
                    395:                        print("level:",2); print(I,2);
                    396:                        print(" comps:",2); print(length(RL[I]),2); print("");
                    397:                        QI = QL[I-1]; RI = RL[I]; NI = length(RI);
                    398:                        for ( J = 0, TI = []; J < NI; J++ ) {
                    399:                                TIJ = extract_qj(QI,V,RI[J],Rad,Mod,SI,Colon,-1);
                    400:                                TI = cons(TIJ,TI);
                    401:                        }
                    402:                        TI = reverse(TI); TL = cons(TI,TL);
                    403:                }
                    404:                TL = reverse(TL);
                    405:        }
                    406:        return TL;
                    407: }
                    408:
                    409: def colon_prime_dec(G,IntQ,V) {
                    410:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    411:        if ( type(Lexdec=getopt(lexdec)) == -1 ) Lexdec = 0;
                    412:        if ( type(Para=getopt(para)) == -1 ) Para = 0;
                    413:        if ( !Para ) {
                    414:                print("colon_pd:",2); print(length(IntQ),2);
                    415:        }
                    416:        if ( !Mod ) M = mingen(IntQ,V);
                    417:        else M = IntQ;
                    418:        if ( Para ) {
                    419:                L = length(M);
                    420:                for ( Task = [], J = 0; J < L; J++ )
                    421:                        if ( gen_nf(M[J],G,V,Ord,Mod) ) {
                    422:                                T = ["noro_pd.call_colon",G,M[J],V,Mod,1];
                    423:                                Task = cons(T,Task);
                    424:                        }
                    425:                Task = reverse(Task);
                    426:                R = para_exec(Para,Task);
                    427:                R = pd_simp_comp(R,V|mod=Mod); L = length(R);
                    428:
                    429:                for ( Task = [], J = 0; J < L; J++ ) {
                    430:                        T = ["noro_pd.call_prime_dec",R[J],V,1,Lexdec,Mod];
                    431:                        Task = cons(T,Task);
                    432:                }
                    433:                Task = reverse(Task);
                    434:                R = para_exec(Para,Task);
                    435:
                    436:                for ( Pt = [], T = R; T != []; T = cdr(T) ) Pt = append(Pt,car(T));
                    437:        } else {
                    438:                for ( R = [], T = M; T != []; T = cdr(T) ) {
                    439:                        Ci = colon(G,car(T),V|isgb=1,mod=Mod);
                    440:                        R = cons(Ci,R);
                    441:                }
                    442:                print("->",2); print(length(M),2);
                    443:                R = pd_simp_comp(R,V|mod=Mod);
                    444:                print("->",2); print(length(R));
                    445: #if 1
                    446:                for ( Pt = [], T = R; T != []; T = cdr(T) ) {
                    447:                        Pi = prime_dec(car(T),V|indep=1,lexdec=Lexdec,mod=Mod);
                    448:                        Pt = append(Pt,Pi);
                    449:                }
                    450: #else
                    451:                J = ideal_list_intersection(R,V,0|mod=Mod);
                    452:                Pt = prime_dec(J,V|indep=1,lexdec=Lexdec,mod=Mod);
                    453: #endif
                    454:        }
                    455: #if 1
                    456:        Pt = pd_simp_comp(Pt,V|first=1,mod=Mod);
                    457: #endif
                    458:        return Pt;
                    459: }
                    460:
                    461: def call_colon(G,F,V,Mod,IsGB)
                    462: {
                    463:        return colon(G,F,V|isgb=1,mod=Mod);
                    464: }
                    465:
                    466: def call_prime_dec(G,V,Indep,Lexdec,Mod)
                    467: {
                    468:        if ( type(G[0]) != 1 )
                    469:                Pi = prime_dec(G,V|indep=Indep,lexdec=Lexdec,mod=Mod);
                    470:        else
                    471:                Pi = [];
                    472:        return Pi;
                    473: }
                    474:
                    475: def extract_qj(Q,V,QL,Rad,Mod,SI,Colon,Level)
                    476: {
                    477:        SIFList=[noro_pd.find_ssi0, noro_pd.find_ssi1, noro_pd.find_ssi2];
                    478:        SIF = SIFList[SI];
                    479:        G = QL[0]; P = QL[1]; PV = QL[2];
                    480:        if ( Q != [1] ) {
                    481:                C = Colon ? ideal_colon(G,Q,V|mod=Mod) : P;
                    482:                Ok = (*SIF)(C,G,Q,Rad,V,0|mod=Mod);
                    483:        } else
                    484:                Ok = [];
                    485:        V0 = setminus(V,PV);
                    486:        HJ = elim_gb(append(Ok,G),V0,PV,Mod,[[0,length(V0)],[0,length(PV)]]);
                    487:        HJ = contraction(HJ,V0|mod=Mod);
1.3       noro      488:        IJ = fast_gb(HJ,V,Mod,Ord);
1.1       noro      489:        return Level >= 0 ? [Level,[IJ,P]] : [IJ,P];
                    490: }
                    491:
                    492: def power(A,I) { return A^I; }
                    493:
                    494:
                    495: /* functions for computating a separing ideal  */
                    496: /* C=G:Q, Rad=rad(Q), return J s.t. Q cap (G+J) = G */
                    497:
                    498: def find_si0(C,G,Q,Rad,V,Ord) {
                    499:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    500:        for ( CI = C, I = 1; ; I++ ) {
                    501:                for ( T = CI, S = []; T != []; T = cdr(T) )
                    502:                        if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
                    503:                if ( S == [] )
                    504:                        error("find_si0 : cannot happen");
                    505:                G1 = append(S,G);
                    506:                Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    507:                /* check whether (Q cap (G+S)) = G */
                    508:                if ( gen_gb_comp(Int,G,Mod) ) { print([0]); return reverse(S); }
                    509:                CI = ideal_product(CI,C,V|mod=Mod);
                    510:        }
                    511: }
                    512:
                    513: def find_si1(C,G,Q,Rad,V,Ord) {
                    514:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    515:        for ( T = C, S = []; T != []; T = cdr(T) )
                    516:                if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
                    517:        if ( S == [] )
                    518:                error("find_si1 : cannot happen");
                    519:        G1 = append(S,G);
                    520:        Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    521:        /* check whether (Q cap (G+S)) = G */
                    522:        if ( gen_gb_comp(Int,G,Mod) ) { print([0]); return reverse(S); }
                    523:
                    524:        C = qsort(C,noro_pd.comp_tdeg);
                    525:
                    526:        Tmp = ttttt; TV = cons(Tmp,V); Ord1 = [[0,1],[Ord,length(V)]];
                    527:        Int0 = incremental_gb(append(vtol(ltov(G)*Tmp),vtol(ltov(Q)*(1-Tmp))),
                    528:                TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
                    529:        Dp = dp_gr_print(); dp_gr_print(0);
                    530:        for ( T = C, S = []; T != []; T = cdr(T) ) {
                    531:                if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
                    532:                Ui = U = car(T);
                    533:                for ( I = 1; ; I++ ) {
                    534:                        G1 = cons(Ui,G);
                    535:                        Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    536:                        if ( gen_gb_comp(Int,G,Mod) ) break;
                    537:                        else
                    538:                                Ui = gen_nf(Ui*U,G,V,Ord,Mod);
                    539:                }
                    540:                print([length(T),I],2);
                    541:                Int1 = incremental_gb(append(Int0,[Tmp*Ui]),TV,Ord1
                    542:                        |gbblock=[[0,length(Int0)]],mod=Mod);
                    543:                Int = elimination(Int1,V);
                    544:                if ( !gen_gb_comp(Int,G,Mod) ) {
                    545:                        break;
                    546:                } else {
                    547:                        Int0 = Int1;
                    548:                        S = cons(Ui,S);
                    549:                }
                    550:        }
                    551:        print("");
                    552:        dp_gr_print(Dp);
                    553:        return reverse(S);
                    554: }
                    555:
                    556: def find_si2(C,G,Q,Rad,V,Ord) {
                    557:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    558:        for ( T = C, S = []; T != []; T = cdr(T) )
                    559:                if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
                    560:        if ( S == [] )
                    561:                error("find_si2 : cannot happen");
                    562:        G1 = append(S,G);
                    563:        Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    564:        /* check whether (Q cap (G+S)) = G */
                    565:        if ( gen_gb_comp(Int,G,Mod) ) { print([0]); return reverse(S); }
                    566:
                    567:        C = qsort(C,noro_pd.comp_tdeg);
                    568:
                    569:        Dp = dp_gr_print(); dp_gr_print(0);
                    570:        Tmp = ttttt; TV = cons(Tmp,V); Ord1 = [[0,1],[Ord,length(V)]];
                    571:        Int0 = incremental_gb(append(vtol(ltov(G)*Tmp),vtol(ltov(Q)*(1-Tmp))),
                    572:                TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
                    573:        for ( T = C, S = []; T != []; T = cdr(T) ) {
                    574:                if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
                    575:                Ui = U = car(T);
                    576:                for ( I = 1; ; I++ ) {
                    577:                        Int1 = incremental_gb(append(Int0,[Tmp*Ui]),TV,Ord1
                    578:                                        |gbblock=[[0,length(Int0)]],mod=Mod);
                    579:                        Int = elimination(Int1,V);
                    580:                        if ( gen_gb_comp(Int,G,Mod) ) break;
                    581:                        else
                    582:                                Ui = gen_nf(Ui*U,G,V,Ord,Mod);
                    583:                }
                    584:                print([length(T),I],2);
                    585:                S = cons(Ui,S);
                    586:        }
                    587:        S = qsort(S,noro_pd.comp_tdeg);
                    588:        print("");
                    589:        End = Len = length(S);
                    590:
                    591:        Tmp = ttttt; TV = cons(Tmp,V); Ord1 = [[0,1],[Ord,length(V)]];
                    592:        Prev = 1;
                    593:        G1 = append(G,[S[0]]);
                    594:        Int0 = incremental_gb(append(vtol(ltov(G1)*Tmp),vtol(ltov(Q)*(1-Tmp))),
                    595:                TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
                    596:        if ( End > 1 ) {
                    597:                Cur = 2;
                    598:                while ( Prev < Cur ) {
                    599:                        for ( St = [], I = Prev; I < Cur; I++ ) St = cons(Tmp*S[I],St);
                    600:                        Int1 = incremental_gb(append(Int0,St),TV,Ord1
                    601:                                |gbblock=[[0,length(Int0)]],mod=Mod);
                    602:                        Int = elimination(Int1,V);
                    603:                        if ( gen_gb_comp(Int,G,Mod) ) {
                    604:                                print([Cur],2);
                    605:                                Prev = Cur;
                    606:                                Cur = Cur+idiv(End-Cur+1,2);
                    607:                                Int0 = Int1;
                    608:                        } else {
                    609:                                End = Cur;
                    610:                                Cur = Prev + idiv(Cur-Prev,2);
                    611:                        }
                    612:                }
                    613:                for ( St = [], I = 0; I < Prev; I++ ) St = cons(S[I],St);
                    614:        } else
                    615:                St = [S[0]];
                    616:        print("");
                    617:        for ( I = Prev; I < Len; I++ ) {
                    618:                Int1 = incremental_gb(append(Int0,[Tmp*S[I]]),TV,Ord1
                    619:                        |gbblock=[[0,length(Int0)]],mod=Mod);
                    620:                Int = elimination(Int1,V);
                    621:                if ( gen_gb_comp(Int,G,Mod) ) {
                    622:                        print([I],2);
                    623:                        St = cons(S[I],St);
                    624:                        Int0 = Int1;
                    625:                }
                    626:        }
                    627:        Ok = reverse(St);
                    628:        print("");
                    629:        print([length(S),length(Ok)]);
                    630:        dp_gr_print(Dp);
                    631:        return Ok;
                    632: }
                    633:
                    634: /* functions for computing a saturated separating ideal */
                    635:
                    636: def find_ssi0(C,G,Q,Rad,V,Ord) {
                    637:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    638:        if ( type(Reduce=getopt(red)) == -1 ) Reduce = 0;
                    639:        for ( T = C, S = []; T != []; T = cdr(T) )
                    640:                if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
                    641:        if ( S == [] )
                    642:                error("find_ssi0 : cannot happen");
                    643:        G1 = append(S,G);
                    644:        Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    645:        /* check whether (Q cap (G+S)) = G */
                    646:        if ( gen_gb_comp(Int,G,Mod) ) { print([0]); return reverse(S); }
                    647:
                    648:        if ( Reduce ) {
                    649:                for ( T = C, U = []; T != []; T = cdr(T) )
                    650:                        if ( gen_nf(car(T),Rad,V,Ord,Mod) ) U = cons(car(T),U);
                    651:                U = reverse(U);
                    652:        } else
                    653:                U = C;
                    654:
                    655:        for ( I = 1; ; I++ ) {
                    656:                print([I],2);
                    657:                S = map(power,U,I);
                    658:                G1 = append(S,G);
                    659:                Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    660:                /* check whether (Q cap (G+S)) = G */
                    661:                if ( gen_gb_comp(Int,G,Mod) ) { print(""); return reverse(S); }
                    662:        }
                    663: }
                    664:
                    665: def find_ssi1(C,G,Q,Rad,V,Ord) {
                    666:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    667:        if ( type(Reduce=getopt(red)) == -1 ) Reduce = 0;
                    668:        for ( T = C, S = []; T != []; T = cdr(T) )
                    669:                if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
                    670:        if ( S == [] )
                    671:                error("find_ssi1 : cannot happen");
                    672:        G1 = append(S,G);
                    673:        Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    674:        /* check whether (Q cap (G+S)) = G */
                    675:        if ( gen_gb_comp(Int,G,Mod) ) { print([0]); return reverse(S); }
                    676:
                    677:        dp_ord(Ord); DC = map(dp_ptod,C,V);
                    678:        DC = qsort(DC,noro_pd.comp_tord); C = map(dp_dtop,DC,V);
                    679:        print(length(C),2);
                    680:        if ( Reduce ) {
                    681:                SC = map(sq,C,Mod);
                    682:                SC = reverse(SC); C = reverse(C);
                    683:                for ( T = C, C1 = [],  R1 = append(SC,Rad); T != []; T = cdr(T) ) {
                    684:                        R0 = car(R1); R1 = cdr(R1);
                    685:                        if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
                    686:                        if ( radical_membership(R0,R1,V|mod=Mod)  ) {
                    687:                                C1 = cons(car(T),C1);
                    688:                                R1 = append(R1,[R0]);
                    689:                        }
                    690:                }
                    691:                print("->",0); print(length(C1),2);
                    692:                C = C1;
                    693:        }
                    694:        print(" ",2);
                    695:
                    696:        Tmp = ttttt; TV = cons(Tmp,V); Ord1 = [[0,1],[Ord,length(V)]];
                    697:        Int0 = incremental_gb(append(vtol(ltov(G)*Tmp),vtol(ltov(Q)*(1-Tmp))),
                    698:                TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
                    699:        Dp = dp_gr_print(); dp_gr_print(0);
                    700:        for ( J = 0, T = C, S = [], GS = G; T != []; T = cdr(T), J++ ) {
                    701:                if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
                    702:                Ui = U = car(T);
                    703:                for ( I = 1; ; I++ ) {
                    704:                        Int1 = nd_gr(append(Int0,[Tmp*Ui]),TV,Mod,Ord1
                    705:                                |gbblock=[[0,length(Int0)]],newelim=1);
                    706:                        if ( Int1 ) {
                    707:                                Int = elimination(Int1,V);
                    708:                                if ( gen_gb_comp(Int,G,Mod) ) break;
                    709:                        }
                    710:                        print("x",2);
                    711:                        Ui = gen_nf(Ui*U,G,V,Ord,Mod);
                    712:                }
                    713:                print(J,2);
                    714:                Int0 = Int1;
                    715:                S = cons(Ui,S);
                    716:        }
                    717:        print("");
                    718:        dp_gr_print(Dp);
                    719:        return reverse(S);
                    720: }
                    721:
                    722: def find_ssi2(C,G,Q,Rad,V,Ord) {
                    723:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    724:        if ( type(Reduce=getopt(red)) == -1 ) Reduce = 0;
                    725:        for ( T = C, S = []; T != []; T = cdr(T) )
                    726:                if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
                    727:        if ( S == [] )
                    728:                error("find_ssi2 : cannot happen");
                    729:        G1 = append(S,G);
                    730:        Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
                    731:        /* check whether (Q cap (G+S)) = G */
                    732:        if ( gen_gb_comp(Int,G,Mod) ) { print([0]); return reverse(S); }
                    733:
                    734: #if 0
                    735:        dp_ord(Ord); DC = map(dp_ptod,C,V);
                    736:        DC = qsort(DC,noro_pd.comp_tord); C = map(dp_dtop,DC,V);
                    737: #else
                    738:        C = qsort(C,noro_pd.comp_tdeg);
                    739: #endif
                    740:        if ( Reduce ) {
                    741:                for ( T = C, C1 = [], R1 = Rad; T != []; T = cdr(T) ) {
                    742:                        if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
                    743:                        if ( radical_membership(car(T),R1,V|mod=Mod)  ) {
                    744:                                C1 = cons(car(T),C1);
                    745:                                R1 = cons(sq(car(T),Mod),R1);
                    746:                        }
                    747:                }
                    748:                print(["C",length(C),"->",length(C1)]);
                    749:                C = reverse(C1);
                    750:        }
                    751:        for ( T = C, S = []; T != []; T = cdr(T) ) {
                    752:                if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
                    753:                Ui = U = car(T);
                    754:                S = cons([Ui,U],S);
                    755:        }
                    756:        S = qsort(S,noro_pd.comp_tdeg_first);
                    757:        print("");
                    758:
                    759:        Dp = dp_gr_print(); dp_gr_print(0);
                    760:        Tmp = ttttt; TV = cons(Tmp,V); Ord1 = [[0,1],[Ord,length(V)]];
                    761:        Int0 = incremental_gb(append(vtol(ltov(G)*Tmp),vtol(ltov(Q)*(1-Tmp))),
                    762:                TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
                    763:        OK = [];
                    764:        while ( S != [] ) {
                    765:                Len = length(S); print("remaining gens : ",0); print(Len);
                    766:                S1 = [];
                    767:                for ( Start = Prev = 0; Start < Len; Start = Prev ) {
                    768:                        Cur = Start+1;
                    769:                        print(Start,2);
                    770:                        while ( Prev < Len ) {
                    771:                                for ( St = [], I = Prev; I < Cur; I++ ) St = cons(Tmp*S[I][0],St);
                    772:                                Int1 = nd_gr(append(Int0,St),TV,Mod,Ord1|gbblock=[[0,length(Int0)]],newelim=1);
                    773:                                if ( !Int1 ) {
                    774:                                        print("x",0); break;
                    775:                                }
                    776:                                Int = elimination(Int1,V);
                    777:                                if ( gen_gb_comp(Int,G,Mod) ) {
                    778:                                        print([Prev,Cur-1],2);
                    779:                                        Prev = Cur;
                    780:                                        Cur += (Prev-Start)+1;
                    781:                                        if ( Cur > Len ) Cur = Len;
                    782:                                        Int0 = Int1;
                    783:                                } else
                    784:                                        break;
                    785:                        }
                    786:                        for ( I = Start; I < Prev; I++ ) OK = cons(S[I][0],OK);
                    787:                        if ( Prev == Start ) {
                    788:                                Ui = S[I][0]; U = S[I][1];
                    789:                                Ui = gen_nf(Ui*U,G,V,Ord,Mod);
                    790:                                S1 = cons([Ui,U],S1);
                    791:                                Prev++;
                    792:                        }
                    793:                }
                    794:                S = reverse(S1);
                    795:                print("");
                    796:        }
                    797:        print("");
                    798:        OK = reverse(OK);
                    799:        dp_gr_print(Dp);
                    800:        return OK;
                    801: }
                    802:
                    803: def iso_comp(G,L,V,Ord)
                    804: {
                    805:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    806:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                    807:        if ( type(Iso=getopt(iso)) == -1 ) Iso = 0;
                    808:        if ( type(Para=getopt(para)) == -1 ) Para = 0;
                    809:        if ( type(Q=getopt(intq)) == -1 ) Q = 0;
                    810:        if ( type(S=getopt(sep)) == -1 ) S = 0;
                    811:
                    812:        if ( !S ) S = separator(L,V|mod=Mod);
                    813:        N = length(L);
                    814:        print("comps : ",2); print(N); print("",2);
                    815:        if ( Para ) {
                    816:                Task = [];
                    817:                for ( I = 0; I < N; I++ ) {
                    818:                        T = ["noro_pd.locsat",G,V,L[I],S[I],Mod,IsGB,Iso,Q];
                    819:                        Task = cons(T,Task);
                    820:                }
                    821:                Task = reverse(Task);
                    822:                R = para_exec(Para,Task);
                    823:                return R;
                    824:        } else {
                    825:                for ( I = 0, R = []; I < N; I++ ) {
                    826:                        QI = locsat(G,V,L[I],S[I],Mod,IsGB,Iso,Q);
                    827:                        if ( type(QI[0][0])==1 )
                    828:                                error("iso_comp : cannot happen");
                    829:                        print(".",2);
                    830:                        R = cons(QI,R);
                    831:                }
                    832:                print("");
                    833:                return reverse(R);
                    834:        }
                    835: }
                    836:
                    837: def locsat(G,V,L,S,Mod,IsGB,Iso,Q)
                    838: {
                    839:        P = L[0]; PV = L[1]; V0 = setminus(V,PV);
                    840:        if ( Iso==1 ) {
                    841:                QI = sat(G,S,V|isgb=IsGB,mod=Mod);
                    842:                GI = elim_gb(QI,V0,PV,Mod,[[0,length(V0)],[0,length(PV)]]);
1.3       noro      843:                GI = fast_gb(contraction(GI,V0|mod=Mod,allv=V),V,Mod,0);
1.1       noro      844:        } else if ( Iso==0 ) {
                    845:                HI = elim_gb(G,V0,PV,Mod,[[0,length(V0)],[0,length(PV)]]);
1.3       noro      846:                GI = fast_gb(contraction(HI,V0|mod=Mod,allv=V),V,Mod,0);
1.1       noro      847:                GI = sat(GI,S,V|isgb=IsGB,mod=Mod);
                    848:        } else if ( Iso==2 ) {
                    849:                HI = elim_gb(G,V0,PV,Mod,[[0,length(V0)],[0,length(PV)]]);
                    850:                TV = ttttt;
                    851:                if ( Mod )
1.3       noro      852:                        GI = fast_gb(cons(TV*S-1,HI),cons(TV,V0),Mod,[[0,1],[0,length(V0)]]);
1.1       noro      853:                else
1.3       noro      854:                        GI = nd_f4_trace(append(HI,[TV*S-1]),cons(TV,V0),
1.1       noro      855:                                1,1,[[0,1],[0,length(V0)]]|gbblock=[[0,length(HI)]]);
                    856:                GI = elimination(GI,V);
1.3       noro      857:                GI = fast_gb(contraction(GI,V0|mod=Mod,allv=V),V,Mod,0);
1.1       noro      858:        } else if ( Iso==3 ) {
                    859:                GI = sat(G,S,V|isgb=IsGB,mod=Mod);
                    860:        }
                    861:        if ( Q )
                    862:                GI = ideal_intersection(Q,GI,V,0|mod=Mod);
                    863:        return [GI,P,PV];
                    864: }
                    865:
                    866: /* SL prime decomposition */
                    867:
                    868: def prime_dec(B,V)
                    869: {
                    870:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    871:        if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
                    872:        if ( type(LexDec=getopt(lexdec)) == -1 ) LexDec = 0;
                    873:        if ( type(Rad=getopt(radical)) == -1 ) Rad = 0;
                    874:        B = map(sq,B,Mod);
                    875:        if ( LexDec )
                    876:                PD = lex_predec1(B,V|mod=Mod);
                    877:        else
                    878:                PD = [B];
                    879:        if ( length(PD) > 1 ) {
                    880:                G = ideal_list_intersection(PD,V,0|mod=Mod);
                    881:                PD = pd_remove_redundant_comp(G,PD,V,0|mod=Mod);
                    882:        }
                    883:        R = []; RL = [];
                    884:        for ( T = PD; T != []; T = cdr(T) ) {
1.5     ! noro      885:                PDT = prime_dec_main(car(T),V|indep=Indep,mod=Mod,radical=Rad);
        !           886:         if ( Rad ) {
        !           887:                  R = append(R,PDT[0]);
        !           888:                  GT = fast_gb(PDT[1],V,Mod,0);
        !           889:                  RL = append(RL,[GT]);
        !           890:         } else {
        !           891:           R = append(R,PDT);
        !           892:         }
1.1       noro      893:        }
                    894:        if ( LexDec ) R = pd_simp_comp(R,V|first=Indep,mod=Mod);
                    895:        if ( Rad ) {
                    896:                G = ideal_list_intersection(RL,V,0|mod=Mod);
                    897:                return [R,G];
                    898:        } else return R;
                    899: }
                    900:
                    901: def prime_dec2(B,V)
                    902: {
                    903:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    904:        if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
                    905:        if ( type(LexDec=getopt(lexdec)) == -1 ) LexDec = 0;
                    906:        if ( type(Rad=getopt(radical)) == -1 ) Rad = 0;
                    907:        if ( type(Para=getopt(para)) == -1 || type(Para) != 4 ) Para = [];
                    908:        B = map(sq,B,Mod);
                    909:        if ( LexDec )
                    910:                PD = lex_predec1(B,V|mod=Mod);
                    911:        else
                    912:                PD = [B];
                    913:        if ( length(PD) > 1 ) {
                    914:                G = ideal_list_intersection(PD,V,0|mod=Mod);
                    915:                PD = pd_remove_redundant_comp(G,PD,V,0|mod=Mod);
                    916:        }
                    917:        R = [];
                    918:        for ( T = PD; T != []; T = cdr(T) )
                    919:                R = append(prime_dec_main2(car(T),V|indep=Indep,mod=Mod,para=Para),R);
                    920:        if ( Indep ) {
                    921:                G = ideal_list_intersection(map(first,R),V,0|mod=Mod);
                    922:                R = pd_simp_comp(R,V|first=1,mod=Mod);
                    923:        } else {
                    924:                G = ideal_list_intersection(R,V,0|mod=Mod);
                    925:                R = pd_simp_comp(R,V|mod=Mod);
                    926:        }
                    927:        return Rad ? [R,G] : R;
                    928: }
                    929:
                    930: /* returns [PD,rad(I)] */
                    931:
                    932: def prime_dec_main(B,V)
                    933: {
                    934:        Tpint = RTpint = 0;
                    935:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    936:        if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
1.5     ! noro      937:        if ( type(Rad=getopt(radical)) == -1 ) Rad = 0;
1.1       noro      938:        G = fast_gb(B,V,Mod,0);
1.5     ! noro      939:     if ( zero_dim(G,V,0) ) {
        !           940:       PD = zprimecomp(G,V,Indep|mod=Mod);
        !           941:       if ( Rad ) {
        !           942:                Int = ideal_list_intersection(map(first,PD),V,0|mod=Mod);
        !           943:         return [PD,Int];
        !           944:       } else
        !           945:         return PD;
        !           946:     }
1.1       noro      947:        IntP = [1];
                    948:        PD = [];
                    949:        DG = ltov(map(dp_ptod,G,V));
                    950:        for ( Ind = [], I = length(G)-1; I >= 0; I-- ) Ind = cons(I,Ind);
                    951:        if ( Mod ) DG = map(dp_mod,DG,Mod,[]);
                    952:        while ( 1 ) {
                    953:                print([length(PD)],2);
                    954:                /* rad(G) subset IntP */
                    955:                /* check if IntP subset rad(G) */
                    956:                /* print([length(PD),length(IntP)],2); */
                    957:                for ( T = IntP; T != []; T = cdr(T) )
                    958:                        if ( (G0 = radical_membership_sat(car(T),G,V|mod=Mod,isgb=1,dg=[DG,Ind])) ) {
                    959:                                F = car(T);
                    960:                                break;
                    961:                        }
                    962:                if ( T == [] ) {
                    963:                        print(["pint",Tpint,"rpint",RTpint]);
1.5     ! noro      964:                        return Rad ? [PD,IntP] : PD;
1.1       noro      965:                }
                    966:                PD0 = zprimecomp(G0,V,Indep|mod=Mod);
                    967:                Int = ideal_list_intersection(Indep?map(first,PD0):PD0,V,0|mod=Mod);
                    968:                PD = append(PD,PD0);
                    969: #if 1
                    970: T0=time();
                    971:                IntP = ideal_intersection_m(IntP,Int,V,0|mod=Mod);
                    972:                dp_ord(0); DC = map(dp_ptod,IntP,V);
                    973:                DC = qsort(DC,noro_pd.comp_tord); IntP = map(dp_dtop,DC,V);
                    974: ACCUM_TIME(Tpint,RTpint)
                    975: #else
                    976:                IntP = ideal_intersection(IntP,Int,V,0|mod=Mod,gbblock=[[0,length(IntP)]]);
                    977: #endif
                    978:        }
                    979: }
                    980:
                    981: localf callsat,callzcomp;
                    982:
                    983: def callsat(F,G,V,Mod,DG)
                    984: {
                    985:        return radical_membership(F,G,V|mod=Mod,isgb=1,dg=DG,sat=1);
                    986: }
                    987:
                    988: def callzcomp(F,V,Indep,Mod)
                    989: {
                    990:        PD0 = zprimecomp(F,V,Indep|mod=Mod);
                    991:        Int = ideal_list_intersection(Indep?map(first,PD0):PD0,V,0|mod=Mod);
                    992:        return [PD0,Int];
                    993: }
                    994:
                    995: def prime_dec_main2(B,V)
                    996: {
                    997:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                    998:        if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
                    999:        if ( type(Para=getopt(para)) == -1 || type(Para) != 4 ) Para = [];
                   1000:        NPara = length(Para);
                   1001:
                   1002:        G = fast_gb(B,V,Mod,0);
                   1003:        IntP = [1];
                   1004:        PD = [];
                   1005:        DG = ltov(map(dp_ptod,G,V));
                   1006:        for ( Ind = [], I = length(G)-1; I >= 0; I-- ) Ind = cons(I,Ind);
                   1007:        if ( Mod ) DG = map(dp_mod,DG,Mod,[]);
                   1008:        if ( NPara )
                   1009:                while ( 1 ) {
                   1010:                        IntPM = mingen(IntP,V);
                   1011:                        for ( T = IntPM, CallSat = []; T != []; T = cdr(T) )
                   1012:                                CallSat = cons(["noro_pd.callsat",car(T),G,V,Mod,[DG,Ind]],CallSat);
                   1013:                        CallSat = reverse(CallSat);
                   1014:                        /* SatL = [[..],0,[...],...] */
                   1015:                        SatL = para_exec(Para,CallSat);
                   1016:                        for ( T = SatL, Sat = []; T != []; T = cdr(T) ) if ( car(T) ) Sat = cons(car(T),Sat);
                   1017:                        if ( Sat == [] ) return PD;
                   1018:                        print(length(Sat),2); print("->",2);
                   1019:                        Sat = remove_identical_comp(Sat|mod=Mod);
                   1020:                        print(length(Sat));
                   1021:                        for ( T = Sat, CallComp = []; T != []; T = cdr(T) )
                   1022:                                CallComp = cons(["noro_pd.callzcomp",car(T),V,Indep,Mod],CallComp);
                   1023:                        CallComp = reverse(CallComp);
                   1024:                        /* PDL = [[PD0,Int],...] */
                   1025:                        PDL = para_exec(Para,CallComp);
                   1026:                        for ( T = PDL; T != []; T = cdr(T) ) PD = append(PD,car(T)[0]);
                   1027:                        Int = ideal_list_intersection(map(second,PDL),V,0|mod=Mod);
                   1028:                        IntP = ideal_intersection(IntP,Int,V,0|mod=Mod,gbblock=[[0,length(IntP)]]);
                   1029:                }
                   1030:        else
                   1031:                while ( 1 ) {
                   1032:                        /* rad(G) subset IntP */
                   1033:                        /* check if IntP subset rad(G) */
                   1034:                        /* print([length(PD),length(IntP)],2); */
                   1035:                        Sat = [];
                   1036:                        IntPM = mingen(IntP,V);
                   1037:                        for ( T = IntPM; T != [] && length(Sat) < 16; T = cdr(T) )
                   1038:                                if ( G0 = radical_membership(car(T),G,V|mod=Mod,isgb=1,dg=[DG,Ind],sat=1) )
                   1039:                                        Sat = cons(G0,Sat);
                   1040:                        if ( Sat == [] ) return PD;
                   1041:                        print(length(Sat),2); print("->",2);
                   1042:                        Sat = remove_identical_comp(Sat|mod=Mod);
                   1043:                        print(length(Sat));
                   1044:                        for ( T = Sat; T != []; T = cdr(T) ) {
                   1045:                                PD0 = zprimecomp(car(T),V,Indep|mod=Mod); PD = append(PD,PD0);
                   1046:                                Int = ideal_list_intersection(Indep?map(first,PD0):PD0,V,0|mod=Mod);
                   1047:                                IntP = ideal_intersection(IntP,Int,V,0|mod=Mod,gbblock=[[0,length(IntP)]]);
                   1048:                        }
                   1049:                }
                   1050: }
                   1051:
                   1052: /* pre-decomposition */
                   1053:
                   1054: def lex_predec1(B,V)
                   1055: {
                   1056:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1057:        G = fast_gb(B,V,Mod,2);
                   1058:        for ( T = G; T != []; T = cdr(T) ) {
                   1059:                F = gen_fctr(car(T),Mod);
                   1060:                if ( length(F) > 2 || length(F) == 2 && F[1][1] > 1 ) {
                   1061:                        for ( R = [], S = cdr(F); S != []; S = cdr(S) ) {
                   1062:                                Ft = car(S)[0];
                   1063:                                Gt = map(ptozp,map(gen_nf,G,[Ft],V,0,Mod));
                   1064:                                R1 = fast_gb(cons(Ft,Gt),V,Mod,0);
                   1065:                                R = cons(R1,R);
                   1066:                        }
                   1067:                        return R;
                   1068:                }
                   1069:        }
                   1070:        return [G];
                   1071: }
                   1072:
                   1073: /* zero-dimensional prime/primary decomosition */
                   1074:
                   1075: def zprimedec(B,V,Mod)
                   1076: {
                   1077:        L = partial_decomp(B,V,Mod);
                   1078:        P = L[0]; NP = L[1];
                   1079:        R = [];
                   1080:        for ( ; P != []; P = cdr(P) ) R = cons(car(car(P)),R);
                   1081:        for ( T = NP; T != []; T = cdr(T) ) {
                   1082:                R1 = complete_decomp(car(T),V,Mod);
                   1083:                R = append(R1,R);
                   1084:        }
                   1085:        return R;
                   1086: }
                   1087:
                   1088: def zprimadec(B,V,Mod)
                   1089: {
                   1090:        L = partial_qdecomp(B,V,Mod);
                   1091:        Q = L[0]; NQ = L[1];
                   1092:        R = [];
                   1093:        for ( ; Q != []; Q = cdr(Q) ) {
                   1094:                T = car(Q); R = cons([T[0],T[1]],R);
                   1095:        }
                   1096:        for ( T = NQ; T != []; T = cdr(T) ) {
                   1097:                R1 = complete_qdecomp(car(T),V,Mod);
                   1098:                R = append(R1,R);
                   1099:        }
                   1100:        return R;
                   1101: }
                   1102:
                   1103: def complete_qdecomp(GD,V,Mod)
                   1104: {
                   1105:        GQ = GD[0]; GP = GD[1]; D = GD[2];
                   1106:        W = vars(GP);
                   1107:        PV = setminus(W,V);
                   1108:        N = length(V); PN = length(PV);
                   1109:        U = find_npos([GP,D],V,PV,Mod);
                   1110:        NV = ttttt;
                   1111:        M = gen_minipoly(cons(NV-U,GQ),cons(NV,V),PV,0,NV,Mod);
                   1112:        M = ppart(M,NV,Mod);
                   1113:        MF = Mod ? modfctr(M,Mod) : fctr(M);
                   1114:        R = [];
                   1115:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                   1116:                S = car(T);
                   1117:                Mt = subst(S[0],NV,U);
                   1118:                GP1 = fast_gb(cons(Mt,GP),W,Mod,0);
                   1119:                GQ1 = fast_gb(cons(Mt^S[1],GQ),W,Mod,0);
                   1120:                if ( PV != [] ) {
                   1121:                        GP1 = elim_gb(GP1,V,PV,Mod,[[0,N],[0,PN]]);
                   1122:                        GQ1 = elim_gb(GQ1,V,PV,Mod,[[0,N],[0,PN]]);
                   1123:                }
                   1124:                R = cons([GQ1,GP1],R);
                   1125:        }
                   1126:        return R;
                   1127: }
                   1128:
                   1129: def partial_qdecomp(B,V,Mod)
                   1130: {
                   1131:        Elim = (Elim=getopt(elim))&&type(Elim)!=-1 ? 1 : 0;
                   1132:        N = length(V);
                   1133:        W = vars(B);
                   1134:        PV = setminus(W,V);
                   1135:        NP = length(PV);
                   1136:        W = append(V,PV);
                   1137:        if ( Elim && PV != [] ) Ord = [[0,N],[0,NP]];
                   1138:        else Ord = 0;
                   1139:        if ( Mod )
                   1140:                B = nd_f4(B,W,Mod,Ord);
                   1141:        else
                   1142:                B = nd_gr_trace(B,W,1,GBCheck,Ord);
                   1143:        Q = []; NQ = [[B,B,vector(N+1)]];
                   1144:        for ( I = length(V)-1; I >= 0; I-- ) {
                   1145:                NQ1 = [];
                   1146:                for ( T = NQ; T != []; T = cdr(T) ) {
                   1147:                        L = partial_qdecomp0(car(T),V,PV,Ord,I,Mod);
                   1148:                        Q = append(L[0],Q);
                   1149:                        NQ1 = append(L[1],NQ1);
                   1150:                }
                   1151:                NQ = NQ1;
                   1152:        }
                   1153:        return [Q,NQ];
                   1154: }
                   1155:
                   1156: def partial_qdecomp0(GD,V,PV,Ord,I,Mod)
                   1157: {
                   1158:        GQ = GD[0]; GP = GD[1]; D = GD[2];
                   1159:        N = length(V); PN = length(PV);
                   1160:        W = append(V,PV);
                   1161:        VI = V[I];
                   1162:        M = gen_minipoly(GQ,V,PV,Ord,VI,Mod);
                   1163:        M = ppart(M,VI,Mod);
                   1164:        if ( Mod )
                   1165:                MF = modfctr(M,Mod);
                   1166:        else
                   1167:                MF = fctr(M);
                   1168:        Q = []; NQ = [];
                   1169:        if ( length(MF) == 2 && MF[1][1] == 1 ) {
                   1170:                D1 = D*1; D1[I] = M;
                   1171:                GQelim = elim_gb(GQ,V,PV,Mod,Ord);
                   1172:                GPelim = elim_gb(GP,V,PV,Mod,Ord);
                   1173:                LD = ldim(GQelim,V);
                   1174:                if ( deg(M,VI) == LD )
                   1175:                        Q = cons([GQelim,GPelim,D1],Q);
                   1176:                else
                   1177:                        NQ = cons([GQelim,GPelim,D1],NQ);
                   1178:                return [Q,NQ];
                   1179:        }
                   1180:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                   1181:                S = car(T); Mt = S[0]; D1 = D*1; D1[I] = Mt;
                   1182:
                   1183:                GQ1 = fast_gb(cons(Mt^S[1],GQ),W,Mod,Ord);
                   1184:                GQelim = elim_gb(GQ1,V,PV,Mod,Ord);
                   1185:                GP1 = fast_gb(cons(Mt,GP),W,Mod,Ord);
                   1186:                GPelim = elim_gb(GP1,V,PV,Mod,Ord);
                   1187:
                   1188:                D1[N] = LD = ldim(GPelim,V);
                   1189:
                   1190:                for ( J = 0; J < N; J++ )
                   1191:                        if ( D1[J] && deg(D1[J],V[J]) == LD ) break;
                   1192:                if ( J < N )
                   1193:                        Q = cons([GQelim,GPelim,D1],Q);
                   1194:                else
                   1195:                        NQ = cons([GQelim,GPelim,D1],NQ);
                   1196:        }
                   1197:        return [Q,NQ];
                   1198: }
                   1199:
                   1200: def complete_decomp(GD,V,Mod)
                   1201: {
                   1202:        G = GD[0]; D = GD[1];
                   1203:        W = vars(G);
                   1204:        PV = setminus(W,V);
                   1205:        N = length(V); PN = length(PV);
                   1206:        U = find_npos(GD,V,PV,Mod);
                   1207:        NV = ttttt;
                   1208:        M = gen_minipoly(cons(NV-U,G),cons(NV,V),PV,0,NV,Mod);
                   1209:        M = ppart(M,NV,Mod);
                   1210:        MF = Mod ? modfctr(M,Mod) : fctr(M);
                   1211:        if ( length(MF) == 2 ) return [G];
                   1212:        R = [];
                   1213:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                   1214:                Mt = subst(car(car(T)),NV,U);
                   1215:                G1 = fast_gb(cons(Mt,G),W,Mod,0);
                   1216:                if ( PV != [] ) G1 = elim_gb(G1,V,PV,Mod,[[0,N],[0,PN]]);
                   1217:                R = cons(G1,R);
                   1218:        }
                   1219:        return R;
                   1220: }
                   1221:
                   1222: def partial_decomp(B,V,Mod)
                   1223: {
                   1224:        Elim = (Elim=getopt(elim))&&type(Elim)!=-1 ? 1 : 0;
                   1225:        N = length(V);
                   1226:        W = vars(B);
                   1227:        PV = setminus(W,V);
                   1228:        NP = length(PV);
                   1229:        W = append(V,PV);
                   1230:        if ( Elim && PV != [] ) Ord = [[0,N],[0,NP]];
                   1231:        else Ord = 0;
                   1232:        if ( Mod )
                   1233:                B = nd_f4(B,W,Mod,Ord);
                   1234:        else
1.3       noro     1235:                B = nd_f4_trace(B,W,1,GBCheck,Ord);
1.1       noro     1236:        P = []; NP = [[B,vector(N+1)]];
                   1237:        for ( I = length(V)-1; I >= 0; I-- ) {
                   1238:                NP1 = [];
                   1239:                for ( T = NP; T != []; T = cdr(T) ) {
                   1240:                        L = partial_decomp0(car(T),V,PV,Ord,I,Mod);
                   1241:                        P = append(L[0],P);
                   1242:                        NP1 = append(L[1],NP1);
                   1243:                }
                   1244:                NP = NP1;
                   1245:        }
                   1246:        return [P,NP];
                   1247: }
                   1248:
                   1249: def partial_decomp0(GD,V,PV,Ord,I,Mod)
                   1250: {
                   1251:        G = GD[0]; D = GD[1];
                   1252:        N = length(V); PN = length(PV);
                   1253:        W = append(V,PV);
                   1254:        VI = V[I];
                   1255:        M = gen_minipoly(G,V,PV,Ord,VI,Mod);
                   1256:        M = ppart(M,VI,Mod);
                   1257:        if ( Mod )
                   1258:                MF = modfctr(M,Mod);
                   1259:        else
                   1260:                MF = fctr(M);
                   1261:        if ( length(MF) == 2 && MF[1][1] == 1 ) {
                   1262:                D1 = D*1;
                   1263:                D1[I] = M;
                   1264:                Gelim = elim_gb(G,V,PV,Mod,Ord);
                   1265:                D1[N] = LD = ldim(Gelim,V);
                   1266:                GD1 = [Gelim,D1];
                   1267:                for ( J = 0; J < N; J++ )
                   1268:                        if ( D1[J] && deg(D1[J],V[J]) == LD )
                   1269:                                return [[GD1],[]];
                   1270:                return [[],[GD1]];
                   1271:        }
                   1272:        P = []; NP = [];
                   1273:        GI = elim_gb(G,V,PV,Mod,Ord);
                   1274:        for ( T = cdr(MF); T != []; T = cdr(T) ) {
                   1275:                Mt = car(car(T));
                   1276:                D1 = D*1;
                   1277:                D1[I] = Mt;
                   1278:                GIt = map(gen_nf,GI,[Mt],V,Ord,Mod);
                   1279:                G1 = cons(Mt,GIt);
                   1280:                Gelim = elim_gb(G1,V,PV,Mod,Ord);
                   1281:                D1[N] = LD = ldim(Gelim,V);
                   1282:                for ( J = 0; J < N; J++ )
                   1283:                        if ( D1[J] && deg(D1[J],V[J]) == LD ) break;
                   1284:                if ( J < N )
                   1285:                        P = cons([Gelim,D1],P);
                   1286:                else
                   1287:                        NP = cons([Gelim,D1],NP);
                   1288:        }
                   1289:        return [P,NP];
                   1290: }
                   1291:
                   1292: /* prime/primary components over rational function field */
                   1293:
                   1294: def zprimacomp(G,V) {
                   1295:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1296:        L = zprimadec(G,V,0|mod=Mod);
                   1297:        R = [];
                   1298:        dp_ord(0);
                   1299:        for ( T = L; T != []; T = cdr(T) ) {
                   1300:                S = car(T);
                   1301:                UQ = contraction(S[0],V|mod=Mod);
                   1302:                UP = contraction(S[1],V|mod=Mod);
                   1303:                R = cons([UQ,UP],R);
                   1304:        }
                   1305:        return R;
                   1306: }
                   1307:
                   1308: def zprimecomp(G,V,Indep) {
                   1309:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1310:        W = maxindep(G,V,0|mod=Mod);
                   1311:        V0 = setminus(V,W);
                   1312:        V1 = append(V0,W);
                   1313: #if 0
                   1314:        O1 = [[0,length(V0)],[0,length(W)]];
                   1315:        G1 = fast_gb(G,V1,Mod,O1);
                   1316:        dp_ord(0);
                   1317: #else
                   1318:        G1 = G;
                   1319: #endif
                   1320:        PD = zprimedec(G1,V0,Mod);
                   1321:        dp_ord(0);
                   1322:        R = [];
                   1323:        for ( T = PD; T != []; T = cdr(T) ) {
                   1324:                U = contraction(car(T),V0|mod=Mod);
                   1325:                U = nd_gr(U,V,Mod,0);
                   1326:                R = cons(Indep?[U,W]:U,R);
                   1327:        }
                   1328:        return R;
                   1329: }
                   1330:
                   1331: def fast_gb(B,V,Mod,Ord)
                   1332: {
1.3       noro     1333:        if ( type(Block=getopt(gbblock)) == -1 ) Block = [];
1.1       noro     1334:        if ( type(NoRA=getopt(nora)) == -1 ) NoRA = 0;
                   1335:        if ( Mod )
                   1336:                G = nd_f4(B,V,Mod,Ord|nora=NoRA);
1.3       noro     1337:        else if ( Chrem )
1.1       noro     1338:                G = map(ptozp,f4_chrem(B,V,Ord));
1.3       noro     1339:        else if ( GBTrace ) {
                   1340:         if ( GBF4 )
                   1341:                  G = nd_f4_trace(B,V,1,1,Ord|nora=NoRA,gbblock=Block);
                   1342:         else
                   1343:                  G = nd_gr_trace(B,V,1,1,Ord|nora=NoRA,gbblock=Block);
1.1       noro     1344:        } else {
1.3       noro     1345:         if ( GBF4 )
                   1346:              G = nd_f4(B,V,0,Ord|nora=NoRA,gbblock=Block);
                   1347:                else
                   1348:           G = nd_gr(B,V,0,Ord|nora=NoRA,gbblock=Block);
1.1       noro     1349:        }
                   1350:        return G;
                   1351: }
                   1352:
                   1353: def incremental_gb(A,V,Ord)
                   1354: {
                   1355:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
1.3       noro     1356:        if ( type(Block=getopt(gbblock)) == -1 ) Block = [];
1.1       noro     1357:        if ( Mod ) {
1.3       noro     1358:          G = fast_gb(A,V,Mod,Ord|gbblock=Block);
1.1       noro     1359:        } else if ( Procs ) {
1.3       noro     1360:                Arg0 = ["nd_f4",A,V,0,Ord];
                   1361:                Arg1 = ["nd_f4_trace",A,V,1,GBCheck,Ord];
1.1       noro     1362:                G = competitive_exec(Procs,Arg0,Arg1);
1.3       noro     1363:        } else
                   1364:          G = fast_gb(A,V,0,Ord|gbblock=Block);
1.1       noro     1365:        return G;
                   1366: }
                   1367:
                   1368: def elim_gb(G,V,PV,Mod,Ord)
                   1369: {
                   1370:        N = length(V); PN = length(PV);
                   1371:        O1 = [[0,N],[0,PN]];
                   1372:        if ( Ord == O1 )
                   1373:                Ord = Ord[0][0];
                   1374:        if ( Mod ) /* XXX */ {
                   1375:                for ( T = G, H = []; T != []; T = cdr(T) )
                   1376:                        if ( car(T) ) H = cons(car(T),H);
                   1377:                G = reverse(H);
1.5     ! noro     1378:         if ( PV == [] )
        !          1379:                  G = nd_f4(G,V,Mod,Ord);
        !          1380:         else
        !          1381:                  G = dp_gr_mod_main(G,V,0,Mod,Ord);
1.1       noro     1382:        } else if ( EProcs ) {
                   1383: #if 1
                   1384:                Arg0 = ["dp_gr_main",G,V,0,0,Ord];
                   1385: #else
                   1386:                Arg0 = ["nd_gr",G,V,0,Ord];
                   1387: #endif
                   1388:                Arg1 = ["noro_pd.nd_gr_rat",G,V,PV,O1,Ord];
                   1389:                G = competitive_exec(EProcs,Arg0,Arg1);
                   1390:        } else if ( GBRat ) {
                   1391:                G1 = nd_gr(G,append(V,PV),0,O1);
                   1392:                if ( GBRat == 1 )
                   1393:                G1 = nd_gr_postproc(G1,V,0,Ord,0|nora=1);
                   1394:                return G1;
                   1395:        } else
                   1396: #if 1
                   1397: #if 0
                   1398:                G = dp_gr_main(G,V,0,0,Ord);
                   1399: #else
1.3       noro     1400:     if ( subset(vars(G),V) )
                   1401:                G = nd_f4_trace(G,V,1,1,Ord);
                   1402:     else
1.1       noro     1403:                G = nd_gr_trace(G,V,1,1,Ord);
                   1404: #endif
                   1405: #else
                   1406:                G = nd_gr(G,V,0,Ord);
                   1407: #endif
                   1408:        return G;
                   1409: }
                   1410:
                   1411: def ldim(G,V)
                   1412: {
                   1413:        O0 = dp_ord(); dp_ord(0);
                   1414:        D = length(dp_mbase(map(dp_ptod,G,V)));
                   1415:        dp_ord(O0);
                   1416:        return D;
                   1417: }
                   1418:
                   1419: /* over Q only */
                   1420:
                   1421: def make_mod_subst(GD,V,PV,HC)
                   1422: {
                   1423:        N = length(V);
                   1424:        PN = length(PV);
                   1425:        G = GD[0]; D = GD[1];
                   1426:        for ( I = 0; ; I = (I+1)%100 ) {
                   1427:                Mod = lprime(I);
                   1428:                S = [];
                   1429:                for ( J = PN-1; J >= 0; J-- )
                   1430:                        S = append([PV[J],random()%Mod],S);
                   1431:                for ( T = HC; T != []; T = cdr(T) )
                   1432:                        if ( !(subst(car(T),S)%Mod) ) break;
                   1433:                if ( T != [] ) continue;
                   1434:                for ( J = 0; J < N; J++ ) {
                   1435:                        M = subst(D[J],S);
                   1436:                        F = modsqfr(M,Mod);
                   1437:                        if ( length(F) != 2 || F[1][1] != 1 ) break;
                   1438:                }
                   1439:                if ( J < N ) continue;
                   1440:                G0 = map(subst,G,S);
                   1441:                return [G0,Mod];
                   1442:        }
                   1443: }
                   1444:
                   1445: def rsgn()
                   1446: {
                   1447:        return random()%2 ? 1 : -1;
                   1448: }
                   1449:
                   1450: def find_npos(GD,V,PV,Mod)
                   1451: {
                   1452:        N = length(V); PN = length(PV);
                   1453:        G = GD[0]; D = GD[1]; LD = D[N];
                   1454:        DH = map(dp_dtop,map(dp_ht,map(dp_ptod,D,V)),V);
                   1455:        Ord0 = dp_ord(); dp_ord(0);
                   1456:        HC = map(dp_hc,map(dp_ptod,G,V));
                   1457:        dp_ord(Ord0);
                   1458:        if ( !Mod ) {
                   1459:                W = append(V,PV);
                   1460:                G1 = nd_gr_trace(G,W,1,GBCheck,[[0,N],[0,PN]]);
                   1461:                L = make_mod_subst([G1,D],V,PV,HC);
                   1462:                return find_npos([L[0],D],V,[],L[1]);
                   1463:        }
                   1464:        N = length(V);
                   1465:        NV = ttttt;
                   1466:        for ( B = 2; ; B++ ) {
                   1467:                for ( J = N-2; J >= 0; J-- ) {
                   1468:                        for ( U = 0, K = J; K < N; K++ ) {
                   1469:                                if ( DH[K] == V[K] ) continue;
                   1470:                                U += rsgn()*((random()%B+1))*V[K];
                   1471:                        }
                   1472: #if 0
                   1473:                        M = minipolym(G,V,0,U,NV,Mod);
                   1474: #else
                   1475:                        M = gen_minipoly(cons(NV-U,G),cons(NV,V),PV,0,NV,Mod);
                   1476: #endif
                   1477:                        if ( deg(M,NV) == LD ) return U;
                   1478:                }
                   1479:        }
                   1480: }
                   1481:
                   1482: def gen_minipoly(G,V,PV,Ord,VI,Mod)
                   1483: {
                   1484:        O0 = dp_ord();
                   1485:        if ( PV == [] ) {
                   1486:                NV = sssss;
                   1487:                if ( Mod )
                   1488:                        M = minipolym(G,V,Ord,VI,NV,Mod);
                   1489:                else
                   1490:                        M = minipoly(G,V,Ord,VI,NV);
                   1491:                dp_ord(O0);
                   1492:                return subst(M,NV,VI);
                   1493:        }
                   1494:        W = setminus(V,[VI]);
                   1495:        PV1 = cons(VI,PV);
                   1496: #if 0
                   1497:        while ( 1 ) {
                   1498:                V1 = append(W,PV1);
                   1499:                if ( Mod )
                   1500:                        G = nd_gr(G,V1,Mod,[[0,1],[0,length(V1)-1]]|nora=1);
                   1501:                else
                   1502:                        G = nd_gr_trace(G,V1,1,GBCheck,[[0,1],[0,length(V1)-1]]|nora=1);
                   1503:                if ( W == [] ) break;
                   1504:                else {
                   1505:                        W = cdr(W);
                   1506:                        G = elimination(G,cdr(V1));
                   1507:                }
                   1508:        }
                   1509: #elif 1
                   1510:        if ( Mod ) {
                   1511:                V1 = append(W,PV1);
                   1512:                G = nd_gr(G,V1,Mod,[[0,length(W)],[0,length(PV1)]]);
                   1513:                G = elimination(G,PV1);
                   1514:        } else {
                   1515:                PV2 = setminus(PV1,[PV1[length(PV1)-1]]);
                   1516:                V2 = append(W,PV2);
                   1517:                G = nd_gr_trace(G,V2,1,GBCheck,[[0,length(W)],[0,length(PV2)]]|nora=1);
                   1518:                G = elimination(G,PV1);
                   1519:        }
                   1520: #else
                   1521:        V1 = append(W,PV1);
                   1522:        if ( Mod )
                   1523:                G = nd_gr(G,V1,Mod,[[0,length(W)],[0,length(PV1)]]|nora=1);
                   1524:        else
                   1525:                G = nd_gr_trace(G,V1,1,GBCheck,[[0,length(W)],[0,length(PV1)]]|nora=1);
                   1526:        G = elimination(G,PV1);
                   1527: #endif
                   1528:        if ( Mod )
                   1529:                G = nd_gr(G,PV1,Mod,[[0,1],[0,length(PV)]]|nora=1);
                   1530:        else
                   1531:                G = nd_gr_trace(G,PV1,1,GBCheck,[[0,1],[0,length(PV)]]|nora=1);
                   1532:        for ( M = car(G), T = cdr(G); T != []; T = cdr(T) )
                   1533:                if ( deg(car(T),VI) < deg(M,VI) ) M = car(T);
                   1534:        dp_ord(O0);
                   1535:        return M;
                   1536: }
                   1537:
                   1538: def indepset(V,H)
                   1539: {
                   1540:        if ( H == [] ) return V;
                   1541:        N = -1;
                   1542:        for ( T = V; T != []; T = cdr(T) ) {
                   1543:                VI = car(T);
                   1544:                HI = [];
                   1545:                for ( S = H; S != []; S = cdr(S) )
                   1546:                        if ( !tdiv(car(S),VI) ) HI = cons(car(S),HI);
                   1547:                RI = indepset(setminus(V,[VI]),HI);
                   1548:                if ( length(RI) > N ) {
                   1549:                        R = RI; N = length(RI);
                   1550:                }
                   1551:        }
                   1552:        return R;
                   1553: }
                   1554:
                   1555: def maxindep(B,V,O)
                   1556: {
                   1557:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1558:        G = fast_gb(B,V,Mod,O);
                   1559:        Old = dp_ord();
                   1560:        dp_ord(O);
                   1561:        H = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
                   1562:        H = map(sq,H,0);
                   1563:        H = nd_gr(H,V,0,0);
                   1564:        H = monodec0(H,V);
                   1565:        N = length(V);
                   1566:        Dep = [];
                   1567:        for ( T = H, Len = N+1; T != []; T = cdr(T) ) {
                   1568:                M = length(car(T));
                   1569:                if ( M < Len ) {
                   1570:                        Dep = [car(T)];
                   1571:                        Len = M;
                   1572:                } else if ( M == Len )
                   1573:                        Dep = cons(car(T),Dep);
                   1574:        }
                   1575:        R = setminus(V,Dep[0]);
                   1576:        dp_ord(Old);
                   1577:        return R;
                   1578: }
                   1579:
                   1580: def maxindep2(B,V,O)
                   1581: {
                   1582:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1583:        G = fast_gb(B,V,Mod,O);
                   1584:        Old = dp_ord();
                   1585:        dp_ord(O);
                   1586:        H = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
                   1587:        H = map(sq,H,0);
                   1588:        H = nd_gr(H,V,0,0);
                   1589:        H = monodec0(H,V);
                   1590:        N = length(V);
                   1591:        Dep = [];
                   1592:        for ( T = H, Len = N+1; T != []; T = cdr(T) ) {
                   1593:                M = length(car(T));
                   1594:                if ( M < Len ) {
                   1595:                        Dep = [car(T)];
                   1596:                        Len = M;
                   1597:                } else if ( M == Len )
                   1598:                        Dep = cons(car(T),Dep);
                   1599:        }
                   1600:        R = [];
                   1601:        for ( T = Dep; T != []; T = cdr(T) )
                   1602:                R = cons(setminus(V,car(T)),R);
                   1603:        dp_ord(Old);
                   1604:        return reverse(R);
                   1605: }
                   1606:
                   1607:
                   1608: /* ideal operations */
                   1609: def contraction(G,V)
                   1610: {
                   1611:        if ( type(AllV=getopt(allv)) == -1 ) AllV = 0;
                   1612:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1613:
                   1614:        if ( RepColon ) return contraction_m(G,V|allv=AllV,mod=Mod);
                   1615:
                   1616:        C = [];
                   1617:        for ( T = G; T != []; T = cdr(T) ) {
                   1618:                C1 = dp_hc(dp_ptod(car(T),V));
                   1619:                S = gen_fctr(C1,Mod);
                   1620:                for ( S = cdr(S); S != []; S = cdr(S) )
                   1621:                        if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
                   1622:        }
                   1623:        W = vars(G);
                   1624:        PV = setminus(W,V);
                   1625:        if ( AllV ) W = AllV;
                   1626:        else W = append(V,PV);
                   1627:        NV = ttttt;
                   1628:        if ( SuccSat ) {
                   1629:                W1 = cons(NV,W);
                   1630:                O1 = [[0,1],[0,length(W)]];
                   1631:                Block = [];
                   1632:                for ( T = C; T != []; T = cdr(T) ) {
                   1633:                        G1 = nd_gr(append(G,[NV*car(T)-1]),W1,Mod,O1|gbblock=Block);
                   1634:                        G = elimination(G1,W);
                   1635:                        Block = [[0,length(G)]];
                   1636:                }
                   1637:        } else {
                   1638:                for ( T = C, S = 1; T != []; T = cdr(T) )
                   1639:                        S *= car(T);
                   1640:                G = saturation([G,NV],S,W|mod=Mod);
                   1641:        }
                   1642:        return G;
                   1643: }
                   1644:
                   1645: def contraction_m(G,V)
                   1646: {
                   1647:        if ( type(AllV=getopt(allv)) == -1 ) AllV = 0;
                   1648:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1649:        C = [];
                   1650:        for ( T = G; T != []; T = cdr(T) ) {
                   1651:                C1 = dp_hc(dp_ptod(car(T),V));
                   1652:                S = gen_fctr(C1,Mod);
                   1653:                for ( S = cdr(S); S != []; S = cdr(S) )
                   1654:                        if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
                   1655:        }
                   1656:        W = vars(G);
                   1657:        PV = setminus(W,V);
                   1658:        if ( AllV ) W = AllV;
                   1659:        else W = append(V,PV);
                   1660:        H = H0 = G;
                   1661:        while ( 1 ) {
                   1662:                for ( T = C; T != []; T = cdr(T) )
                   1663:                        H = map(sdiv,ideal_intersection_m([car(T)],H,W,0),car(T));
                   1664:                H = nd_gr(H,W,0,0);
                   1665:                if ( gb_comp(H0,H) ) break;
                   1666:                else H0 = H;
                   1667:        }
                   1668:        return H;
                   1669: }
                   1670:
                   1671: def ideal_list_intersection(L,V,Ord)
                   1672: {
                   1673:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1674:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                   1675:        N = length(L);
                   1676:        if ( N == 0 ) return [1];
                   1677:        if ( N == 1 )
                   1678:                return IsGB ? L[0] : fast_gb(L[0],V,Mod,Ord);
                   1679:        else {
                   1680:                for ( I = 0, T = [1]; I < N; I++ )
                   1681:                        T = ideal_intersection_m(T,L[I],V,Ord|mod=Mod);
                   1682:                T = nd_gr(T,V,Mod,Ord);
                   1683:                return T;
                   1684:        }
                   1685: }
                   1686:
                   1687: def call_ideal_list_intersection(L,V,Mod,Ord,IsGB)
                   1688: {
                   1689:        return ideal_list_intersection(L,V,Ord|mod=Mod,isgb=IsGB);
                   1690: }
                   1691:
                   1692: def ideal_intersection(A,B,V,Ord)
                   1693: {
                   1694:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
1.3       noro     1695:        if ( type(Block=getopt(gbblock)) == -1 ) Block = [];
1.1       noro     1696:        T = ttttt;
                   1697:        if ( Mod ) {
1.3       noro     1698:          G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1699:                cons(T,V),Mod,[[0,1],[Ord,length(V)]]|gbblock=Block,nora=0);
1.1       noro     1700:        } else
                   1701:        if ( Procs ) {
                   1702:                Arg0 = ["nd_gr",
                   1703:                        append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1704:                        cons(T,V),0,[[0,1],[Ord,length(V)]]];
                   1705:                Arg1 = ["nd_gr_trace",
                   1706:                        append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1707:                        cons(T,V),1,GBCheck,[[0,1],[Ord,length(V)]]];
                   1708:                G = competitive_exec(Procs,Arg0,Arg1);
                   1709:        } else {
1.3       noro     1710:          G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
                   1711:                cons(T,V),0,[[0,1],[Ord,length(V)]]|gbblock=Block,nora=0);
1.1       noro     1712:        }
                   1713:        G0 = elimination(G,V);
                   1714:        if ( 0 && !Procs )
                   1715:                G0 = nd_gr_postproc(G0,V,Mod,Ord,0);
                   1716:        return G0;
                   1717: }
                   1718:
                   1719:
                   1720: def aa(A) { return [A,A]; }
                   1721:
                   1722: def ideal_intersection_m(A,B,V,Ord)
                   1723: {
                   1724:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1725:
                   1726:        dp_ord(Ord);
1.5     ! noro     1727:
        !          1728:     if ( UsualInt )
        !          1729:       return ideal_intersection(A,B,V,Ord|mod=Mod);
        !          1730:
1.1       noro     1731:        DA = map(dp_ptod,A,V); DB = ltov(map(dp_ptod,B,V));
                   1732:        if ( Mod ) {
                   1733:                DA = map(dp_mod,DA,Mod,[]); DB = map(dp_mod,DB,Mod,[]);
                   1734:                setmod(Mod);
                   1735:        }
                   1736:        N = length(B);
                   1737:        for ( Ind = [], I = N-1; I >= 0; I-- ) Ind = cons(I,Ind);
                   1738:        for ( T = DA, C = []; T != []; T = cdr(T) ) {
1.3       noro     1739:                R = Mod?dp_nf_mod(Ind,car(T),DB,1,Mod):dp_nf(Ind,car(T),DB,1);
                   1740:         if ( !R )
                   1741:                  C = cons([0,dp_dtop(car(T),V)],C);
                   1742:         else
                   1743:                  C = cons([dp_dtop(car(T),V),0],C);
                   1744:        }
                   1745:     if ( IntGB ) {
                   1746:          G = nd_gr(append(C,map(aa,B)),V,Mod,[1,Ord]);
                   1747:       for ( T = G, G1 = []; T != []; T = cdr(T) )
                   1748:         if ( !car(T)[0] ) G1 = cons(car(T)[1],G1);
                   1749:       G = reverse(G1);
                   1750:     } else {
                   1751:          G = nd_gr(append(C,map(aa,B)),V,Mod,[1,Ord]|intersect=1);
                   1752:          G = map(second,G);
                   1753:     }
1.1       noro     1754:        return G;
                   1755: }
                   1756:
                   1757: /* returns GB if F notin rad(G) */
                   1758:
                   1759: def radical_membership(F,G,V) {
                   1760:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1761:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                   1762:        if ( type(L=getopt(dg)) == -1 ) L = 0;
                   1763:        if ( type(Sat=getopt(sat)) == -1 ) Sat = 0;
                   1764:        dp_ord(0);
                   1765:        if ( L ) { DG = L[0]; Ind = L[1]; }
                   1766:        else {
                   1767:                DG = ltov(map(dp_ptod,G,V));
                   1768:                if ( Mod ) DG = map(dp_mod,DG,Mod,[]);
                   1769:                for ( Ind = [], I = length(G)-1; I >= 0; I-- ) Ind = cons(I,Ind);
                   1770:        }
                   1771:        DF = dp_ptod(F,V); DFI = dp_ptod(1,V);
                   1772:        if ( Mod ) {
                   1773:                DF = dp_mod(DF,Mod,[]); DFI = dp_mod(DFI,Mod,[]);
                   1774:                setmod(Mod);
                   1775:        }
                   1776:        for ( I = 0; I < 3; I++ ) {
                   1777:                DFI = Mod?dp_nf_mod(Ind,DF*DFI,DG,0,Mod):dp_nf(Ind,DF*DFI,DG,0);
                   1778:                if ( !DFI ) return 0;
                   1779:        }
                   1780:        NV = ttttt;
                   1781:        if ( IsGB )
                   1782:                T = nd_gr(append(G,[NV*F-1]),cons(NV,V),Mod,0
                   1783:                        |gbblock=[[0,length(G)]]);
                   1784:        else
                   1785:                T = nd_gr(append(G,[NV*F-1]),cons(NV,V),Mod,0);
                   1786:        if ( type(car(T)) == 1 ) return 0;
                   1787:        else if ( Sat ) {
                   1788:                G1 = fast_gb(T,cons(NV,V),Mod,[[0,1],[0,length(V)]]);
                   1789:                G0 = elimination(G1,V);
                   1790:                return G0;
                   1791:        } else return [T,NV];
                   1792: }
                   1793:
                   1794: def radical_membership_sat(F,G,V) {
                   1795:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1796:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                   1797:        if ( type(L=getopt(dg)) == -1 ) L = 0;
                   1798:        dp_ord(0);
                   1799:        if ( L ) { DG = L[0]; Ind = L[1]; }
                   1800:        else {
                   1801:                DG = ltov(map(dp_ptod,G,V));
                   1802:                if ( Mod ) DG = map(dp_mod,DG,Mod,[]);
                   1803:                for ( Ind = [], I = length(G)-1; I >= 0; I-- ) Ind = cons(I,Ind);
                   1804:        }
                   1805:        DF = dp_ptod(F,V); DFI = dp_ptod(1,V);
                   1806:        if ( Mod ) {
                   1807:                DF = dp_mod(DF,Mod,[]); DFI = dp_mod(DFI,Mod,[]);
                   1808:                setmod(Mod);
                   1809:        }
                   1810:        for ( I = 0; I < 3; I++ ) {
                   1811:                DFI = Mod?dp_nf_mod(Ind,DF*DFI,DG,0,Mod):dp_nf(Ind,DF*DFI,DG,0);
                   1812:                if ( !DFI ) return 0;
                   1813:        }
                   1814:        NV = ttttt;
                   1815:        if ( IsGB )
                   1816:                T = nd_gr(append(G,[NV*F-1]),cons(NV,V),Mod,[[0,1],[0,length(V)]]
                   1817:                        |gbblock=[[0,length(G)]]);
                   1818:        else
                   1819:                T = nd_gr(append(G,[NV*F-1]),cons(NV,V),Mod,[[0,1],[0,length(V)]]);
                   1820:        if ( type(car(T)) == 1 ) return 0;
                   1821:        G0 = elimination(T,V);
                   1822:        return G0;
                   1823: }
                   1824:
                   1825: def modular_radical_membership(F,G,V) {
                   1826:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1827:        if ( Mod )
                   1828:                return radical_membership(F,G,V|mod=Mod);
                   1829:
                   1830:        F = gen_nf(F,G,V,0,0);
                   1831:        if ( !F ) return 0;
                   1832:        NV = ttttt;
                   1833:        for ( J = 0; ; J++ ) {
                   1834:                Mod = lprime(J);
                   1835:                H = map(dp_hc,map(dp_ptod,G,V));
                   1836:                for ( ; H != []; H = cdr(H) ) if ( !(car(H)%Mod) ) break;
                   1837:                if ( H != [] ) continue;
                   1838:
                   1839:                T = nd_f4(cons(NV*F-1,G),cons(NV,V),Mod,0);
                   1840:                if ( type(car(T)) == 1 ) {
                   1841:                        I = radical_membership_rep(F,G,V,-1,0,Mod);
                   1842:                        I1 = radical_membership_rep(F,G,V,I,0,0);
                   1843:                        if ( I1 > 0 ) return 0;
                   1844:                }
                   1845:                return radical_membership(F,G,V);
                   1846:        }
                   1847: }
                   1848:
                   1849: def radical_membership_rep(F,G,V,Max,Ord,Mod) {
                   1850:        Ft = F;
                   1851:        I = 1;
                   1852:        while ( Max < 0 || I <= Max ) {
                   1853:                Ft = gen_nf(Ft,G,V,Ord,Mod);
                   1854:                if ( !Ft ) return I;
                   1855:                Ft *= F;
                   1856:                I++;
                   1857:        }
                   1858:        return -1;
                   1859: }
                   1860:
                   1861: def ideal_product(A,B,V)
                   1862: {
                   1863:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1864:        dp_ord(0);
                   1865:        DA = map(dp_ptod,A,V);
                   1866:        DB = map(dp_ptod,B,V);
                   1867:        DegA = map(dp_td,DA);
                   1868:        DegB = map(dp_td,DB);
                   1869:        for ( PA = [], T = A, DT = DegA; T != []; T = cdr(T), DT = cdr(DT) )
                   1870:                PA = cons([car(T),car(DT)],PA);
                   1871:        PA = reverse(PA);
                   1872:        for ( PB = [], T = B, DT = DegB; T != []; T = cdr(T), DT = cdr(DT) )
                   1873:                PB = cons([car(T),car(DT)],PB);
                   1874:        PB = reverse(PB);
                   1875:        R = [];
                   1876:        for ( T = PA; T != []; T = cdr(T) )
                   1877:                for ( S = PB; S != []; S = cdr(S) )
                   1878:                        R = cons([car(T)[0]*car(S)[0],car(T)[1]+car(S)[1]],R);
                   1879:        T = qsort(R,noro_pd.comp_by_second);
                   1880:        T = map(first,T);
                   1881:        Len = length(A)>length(B)?length(A):length(B);
                   1882:        Len *= 2;
                   1883:        L = sep_list(T,Len); B0 = L[0]; B1 = L[1];
                   1884:        R = fast_gb(B0,V,Mod,0);
                   1885:        while ( B1 != [] ) {
                   1886:                print(length(B1));
                   1887:                L = sep_list(B1,Len);
                   1888:                B0 = L[0]; B1 = L[1];
                   1889:                R = fast_gb(append(R,B0),V,Mod,0|gbblock=[[0,length(R)]],nora=1);
                   1890:        }
                   1891:        return R;
                   1892: }
                   1893:
                   1894: def saturation(GNV,F,V)
                   1895: {
                   1896:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1897:        G = GNV[0]; NV = GNV[1];
                   1898:        if ( Mod )
1.3       noro     1899:                G1 = fast_gb(cons(NV*F-1,G),cons(NV,V),Mod,[[0,1],[0,length(V)]]);
1.1       noro     1900:        else if ( Procs ) {
1.3       noro     1901:                Arg0 = ["nd_f4_trace",
1.1       noro     1902:                cons(NV*F-1,G),cons(NV,V),0,GBCheck,[[0,1],[0,length(V)]]];
1.3       noro     1903:                Arg1 = ["nd_f4_trace",
1.1       noro     1904:                cons(NV*F-1,G),cons(NV,V),1,GBCheck,[[0,1],[0,length(V)]]];
                   1905:                G1 = competitive_exec(Procs,Arg0,Arg1);
                   1906:        } else
1.3       noro     1907:                G1 = fast_gb(cons(NV*F-1,G),cons(NV,V),0,[[0,1],[0,length(V)]]);
1.1       noro     1908:        return elimination(G1,V);
                   1909: }
                   1910:
                   1911: def sat(G,F,V)
                   1912: {
                   1913:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1914:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                   1915:        NV = ttttt;
                   1916:        if ( Mod )
1.3       noro     1917:                G1 = fast_gb(cons(NV*F-1,G),cons(NV,V),Mod,[[0,1],[0,length(V)]]);
1.1       noro     1918:        else if ( Procs ) {
1.3       noro     1919:                Arg0 = ["nd_f4_trace",
1.1       noro     1920:                cons(NV*F-1,G),cons(NV,V),0,GBCheck,[[0,1],[0,length(V)]]];
1.3       noro     1921:                Arg1 = ["nd_f4_trace",
1.1       noro     1922:                cons(NV*F-1,G),cons(NV,V),1,GBCheck,[[0,1],[0,length(V)]]];
                   1923:                G1 = competitive_exec(Procs,Arg0,Arg1);
                   1924:        } else {
                   1925:                B1 = append(G,[NV*F-1]);
                   1926:                V1 = cons(NV,V);
                   1927:                Ord1 = [[0,1],[0,length(V)]];
                   1928:                if ( IsGB )
1.3       noro     1929:                        G1 = fast_gb(B1,V1,0,Ord1|gbblock=[[0,length(G)]]);
1.1       noro     1930:                else
1.3       noro     1931:                        G1 = fast_gb(B1,V1,0,Ord1);
1.1       noro     1932:        }
                   1933:        return elimination(G1,V);
                   1934: }
                   1935:
                   1936: def isat(B,S,V)
                   1937: {
                   1938:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1939:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                   1940:        F = cdr(fctr(S));
                   1941:        R = B;
                   1942:        for ( T = F; T != []; T = cdr(T) )
                   1943:                R = sat(R,car(T)[0],V|mod=Mod,isgb=IsGB);
                   1944:        return R;
                   1945: }
                   1946:
                   1947: /* homogeneous case only */
                   1948:
                   1949: def sat_ind_var(G,F,V)
                   1950: {
                   1951:        if ( type(Ord=getopt(ord)) == -1 ) Ord = 0;
                   1952:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1953:   V0 = append(setminus(V,[F]),[F]);
                   1954:   G0 = nd_gr(G,V0,Mod,0);
                   1955:   M = 0;
                   1956:   for ( G1 = [], T = G0; T != []; T = cdr(T) ) {
                   1957:     S = car(T);
                   1958:     M1 = mindeg(S,F);
                   1959:     S = sdiv(S,F^M1);
                   1960:     G1 = cons(S,G1);
                   1961:     if ( M1 > M ) M = M1;
                   1962:   }
                   1963:   G1 = nd_gr(G1,V,Mod,Ord);
                   1964:   return [G1,M];
                   1965: }
                   1966:
                   1967: def sat_ind(G,F,V)
                   1968: {
                   1969:        if ( type(Ord=getopt(ord)) == -1 ) Ord = 0;
                   1970:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1971:        NV = ttttt;
                   1972:        F = gen_nf(F,G,V,Ord,Mod);
                   1973:        for ( I = 0, GI = G; ; I++ ) {
                   1974:                G1 = colon(GI,F,V|mod=Mod,ord=Ord);
                   1975:                if ( ideal_inclusion(G1,GI,V,Ord|mod=Mod) )  {
                   1976:                        return [GI,I];
                   1977:                }
                   1978:                else GI = G1;
                   1979:        }
                   1980: }
                   1981:
                   1982: def colon(G,F,V)
                   1983: {
                   1984:        if ( type(Ord=getopt(ord)) == -1 ) Ord = 0;
                   1985:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   1986:        if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
                   1987:        F = gen_nf(F,G,V,Ord,Mod);
                   1988:        if ( !F ) return [1];
                   1989:        if ( IsGB )
                   1990:                T = ideal_intersection(G,[F],V,Ord|gbblock=[[0,length(G)]],mod=Mod);
                   1991:        else
                   1992:                T = ideal_intersection(G,[F],V,Ord|mod=Mod);
                   1993:        Gen = Mod?map(sdivm,T,F,Mod):map(ptozp,map(sdiv,T,F));
                   1994:        return nd_gr(Gen,V,Mod,Ord);
                   1995: }
                   1996:
                   1997: #if 1
                   1998: def ideal_colon(G,F,V)
                   1999: {
                   2000:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2001:        G = nd_gr(G,V,Mod,0);
                   2002:        C = [1];
                   2003:        TV = ttttt;
                   2004:        F = qsort(F,noro_pd.comp_tdeg);
                   2005:        for ( T = F; T != []; T = cdr(T) ) {
                   2006:                S = colon(G,car(T),V|isgb=1,mod=Mod);
                   2007:                if ( type(S[0])!= 1 ) {
                   2008:                        C = nd_gr(append(vtol(ltov(C)*TV),vtol(ltov(S)*(1-TV))),
                   2009:                                cons(TV,V),Mod,[[0,1],[Ord,length(V)]]|gbblock=[[0,length(C)]]);
                   2010:                        C = elimination(C,V);
                   2011:                }
                   2012:        }
                   2013:        return C;
                   2014: }
                   2015: #else
                   2016: def ideal_colon(G,F,V)
                   2017: {
                   2018:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2019:        G = nd_gr(G,V,Mod,0);
                   2020:        for ( T = F, L = []; T != []; T = cdr(T) ) {
                   2021:                C = colon(G,car(T),V|isgb=1,mod=Mod);
                   2022:                if ( type(C[0]) != 1 ) L = cons(C,L);
                   2023:        }
                   2024:        L = reverse(L);
                   2025:        return ideal_list_intersection(L,V,0|mod=Mod);
                   2026: }
                   2027:
                   2028: #endif
                   2029:
                   2030: def member(A,L)
                   2031: {
                   2032:        for ( ; L != []; L = cdr(L) )
                   2033:                if ( car(L) == A ) return 1;
                   2034:        return 0;
                   2035: }
                   2036:
1.3       noro     2037: /* return 1 of A is a subset of B */
                   2038:
                   2039: def subset(A,B)
                   2040: {
                   2041:   for ( ; A != []; A = cdr(A) )
                   2042:     if ( !member(car(A),B) ) return 0;
                   2043:   return 1;
                   2044: }
                   2045:
1.1       noro     2046: def mingen(B,V) {
                   2047:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2048:        Data = nd_gr(B,V,Mod,O|gentrace=1,gensyz=1);
                   2049:     G = Data[0]; STrace = Data[6];
                   2050:        N = length(G);
                   2051:        S = compute_gbsyz(N,V,STrace,Mod);
                   2052:        for ( T = S, R = []; T != []; T = cdr(T) ) {
                   2053:                for ( A = car(T); A1 = dp_rest(A); A = A1);
                   2054:                if ( type(dp_hc(A)) ==1 ) R = cons(dp_etov(A)[0],R);
                   2055:        }
                   2056:        for ( I = 0, U = []; I < N; I++ ) if ( !member(I,R) ) U = cons(G[I],U);
                   2057:        return U;
                   2058: }
                   2059:
                   2060: def compute_gbsyz(N,V,Trace,Mod)
                   2061: {
                   2062:        P = vector(N);
                   2063:        for ( I = 0; I < N; I++ ) P[I] = dp_ptod(x^I,[x]);
                   2064:        for ( U = [], T = Trace; T != []; T = cdr(T) ) {
                   2065:                Ti = car(T);
                   2066:                if ( Ti[0] != -1 ) error("Input is not a GB");
                   2067:                R = recompute_trace(Ti[1],P,V,Mod);
                   2068:                U = cons(R,U);
                   2069:        }
                   2070:        return reverse(U);
                   2071: }
                   2072:
                   2073: def recompute_trace(Ti,P,V,Mod)
                   2074: {
                   2075:   for ( Num = 0, Den = 1; Ti != []; Ti = cdr(Ti) ) {
                   2076:     Sj = car(Ti); Dj = Sj[0]; Ij =Sj[1]; Mj = dp_dtop(Sj[2],V); Cj = Sj[3];
                   2077:     /* Num/Den <- (Dj*Num+Den*Mj*P[Ij])/(Den*Cj) */
                   2078:     if ( Dj ) Num = (Dj*Num+Den*Mj*P[Ij]);
                   2079:     Den *= Cj;
                   2080:   }
                   2081:   return Num;
                   2082: }
                   2083:
                   2084: def ideal_sat(G,F,V)
                   2085: {
                   2086:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2087:        G = nd_gr(G,V,Mod,0);
                   2088:        for ( T = F, L = []; T != []; T = cdr(T) )
                   2089:                L = cons(sat(G,car(T),V|mod=Mod),L);
                   2090:        L = reverse(L);
                   2091:        return ideal_list_intersection(L,V,0|mod=Mod);
                   2092: }
                   2093:
                   2094: def ideal_inclusion(F,G,V,O)
                   2095: {
                   2096:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2097:        for ( T = F; T != []; T = cdr(T) )
                   2098:                if ( gen_nf(car(T),G,V,O,Mod) ) return 0;
                   2099:        return 1;
                   2100: }
                   2101:
                   2102: /* remove redundant components */
                   2103:
                   2104: def qd_simp_comp(QP,V)
                   2105: {
                   2106:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2107:        R = ltov(QP);
                   2108:        N = length(R);
                   2109:        for ( I = 0; I < N; I++ ) {
                   2110:                if ( R[I] ) {
                   2111:                        QI = R[I][0]; PI = R[I][1];
                   2112:                        for ( J = I+1; J < N; J++ )
                   2113:                                if ( R[J] && gen_gb_comp(PI,R[J][1],Mod) ) {
                   2114:                                        QI = ideal_intersection(QI,R[J][0],V,0|mod=Mod);
                   2115:                                        R[J] = 0;
                   2116:                                }
                   2117:                        R[I] = [QI,PI];
                   2118:                }
                   2119:        }
                   2120:        for ( I = N-1, S = []; I >= 0; I-- )
                   2121:                if ( R[I] ) S = cons(R[I],S);
                   2122:        return S;
                   2123: }
                   2124:
                   2125: def qd_remove_redundant_comp(G,Iso,Emb,V,Ord)
                   2126: {
                   2127:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2128:        IsoInt = ideal_list_intersection(map(first,Iso),V,Ord|mod=Mod);
                   2129:        Emb = qd_simp_comp(Emb,V|mod=Mod);
                   2130:        Emb = reverse(qsort(Emb));
                   2131:        A = ltov(Emb); N = length(A);
                   2132:        Pre = IsoInt; Post = vector(N+1);
                   2133:        for ( Post[N] = IsoInt, I = N-1; I >= 1; I-- )
                   2134:                Post[I] = ideal_intersection(Post[I+1],A[I][0],V,Ord|mod=Mod);
                   2135:        for ( I = 0; I < N; I++ ) {
                   2136:                print(".",2);
                   2137:                Int = ideal_intersection(Pre,Post[I+1],V,Ord|mod=Mod);
                   2138:                if ( gen_gb_comp(Int,G,Mod) ) A[I] = 0;
                   2139:                else
                   2140:                        Pre = ideal_intersection(Pre,A[I][0],V,Ord|mod=Mod);
                   2141:        }
                   2142:        for ( T = [], I = 0; I < N; I++ )
                   2143:                if ( A[I] ) T = cons(A[I],T);
                   2144:        return reverse(T);
                   2145: }
                   2146:
                   2147: def pd_simp_comp(PL,V)
                   2148: {
                   2149:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2150:        if ( type(First=getopt(first)) == -1 ) First = 0;
                   2151:        A = ltov(PL); N = length(A);
                   2152:        if ( N == 1 )  return PL;
                   2153:        for ( I = 0; I < N; I++ ) {
                   2154:                if ( !A[I] ) continue;
                   2155:                AI = First?A[I][0]:A[I];
                   2156:                for ( J = 0; J < N; J++ ) {
                   2157:                        if ( J == I || !A[J] ) continue;
                   2158:                        AJ = First?A[J][0]:A[J];
                   2159:                        if ( gen_gb_comp(AI,AJ,Mod) || ideal_inclusion(AI,AJ,V,Ord|mod=Mod) )
                   2160:                                A[J] = 0;
                   2161:                }
                   2162:        }
                   2163:        for ( I = 0, T = []; I < N; I++ ) if ( A[I] ) T = cons(A[I],T);
                   2164:        return reverse(T);
                   2165: }
                   2166:
                   2167: def pd_remove_redundant_comp(G,P,V,Ord)
                   2168: {
                   2169:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2170:        if ( type(First=getopt(first)) == -1 ) First = 0;
                   2171:        if ( length(P) == 1 )  return P;
                   2172:
                   2173:        A = ltov(P); N = length(A);
                   2174:        for ( I = 0; I < N; I++ ) {
                   2175:                if ( !A[I] ) continue;
                   2176:                for ( J = I+1; J < N; J++ )
                   2177:                        if ( A[J] &&
                   2178:                                gen_gb_comp(First?A[I][0]:A[I],First?A[J][0]:A[J],Mod) ) A[J] = 0;
                   2179:        }
                   2180:        for ( I = 0, T = []; I < N; I++ ) if ( A[I] ) T = cons(A[I],T);
                   2181:        A = ltov(reverse(T)); N = length(A);
                   2182:        Pre = [1]; Post = vector(N+1);
                   2183:        for ( Post[N] = [1], I = N-1; I >= 1; I-- )
                   2184:                Post[I] = ideal_intersection(Post[I+1],First?A[I][0]:A[I],V,Ord|mod=Mod);
                   2185:        for ( I = 0; I < N; I++ ) {
                   2186:                Int = ideal_intersection(Pre,Post[I+1],V,Ord|mod=Mod);
                   2187:                if ( gen_gb_comp(Int,G,Mod) ) A[I] = 0;
                   2188:                else
                   2189:                        Pre = ideal_intersection(Pre,First?A[I][0]:A[I],V,Ord|mod=Mod);
                   2190:        }
                   2191:        for ( T = [], I = 0; I < N; I++ ) if ( A[I] ) T = cons(A[I],T);
                   2192:        return reverse(T);
                   2193: }
                   2194:
                   2195: def remove_identical_comp(L)
                   2196: {
                   2197:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2198:        if ( length(L) == 1 )  return L;
                   2199:
                   2200:        A = ltov(L); N = length(A);
                   2201:        for ( I = 0; I < N; I++ ) {
                   2202:                if ( !A[I] ) continue;
                   2203:                for ( J = I+1; J < N; J++ )
                   2204:                        if ( A[J] &&
                   2205:                                gen_gb_comp(A[I],A[J],Mod) ) A[J] = 0;
                   2206:        }
                   2207:        for ( I = 0, T = []; I < N; I++ ) if ( A[I] ) T = cons(A[I],T);
                   2208:        return reverse(T);
                   2209: }
                   2210:
                   2211: /* polynomial operations */
                   2212:
                   2213: def ppart(F,V,Mod)
                   2214: {
                   2215:        if ( !Mod )
                   2216:                G = nd_gr([F],[V],0,0);
                   2217:        else
                   2218:                G = dp_gr_mod_main([F],[V],0,Mod,0);
                   2219:        return G[0];
                   2220: }
                   2221:
                   2222:
                   2223: def sq(F,Mod)
                   2224: {
                   2225:        if ( !F ) return 0;
                   2226:        A = cdr(gen_fctr(F,Mod));
                   2227:        for ( R = 1; A != []; A = cdr(A) )
                   2228:                R *= car(car(A));
                   2229:        return R;
                   2230: }
                   2231:
                   2232: def lcfactor(G,V,O,Mod)
                   2233: {
                   2234:        O0 = dp_ord(); dp_ord(O);
                   2235:        C = [];
                   2236:        for ( T = G; T != []; T = cdr(T) ) {
                   2237:                C1 = dp_hc(dp_ptod(car(T),V));
                   2238:                S = gen_fctr(C1,Mod);
                   2239:                for ( S = cdr(S); S != []; S = cdr(S) )
                   2240:                        if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
                   2241:        }
                   2242:        dp_ord(O0);
                   2243:        return C;
                   2244: }
                   2245:
                   2246: def gen_fctr(F,Mod)
                   2247: {
                   2248:        if ( Mod ) return modfctr(F,Mod);
                   2249:        else return fctr(F);
                   2250: }
                   2251:
                   2252: def gen_mptop(F)
                   2253: {
                   2254:        if ( !F ) return F;
                   2255:        else if ( type(F)==1 )
                   2256:                if ( ntype(F)==5 ) return mptop(F);
                   2257:                else return F;
                   2258:        else {
                   2259:                V = var(F);
                   2260:                D = deg(F,V);
                   2261:                for ( R = 0, I = 0; I <= D; I++ )
                   2262:                        if ( C = coef(F,I,V) ) R += gen_mptop(C)*V^I;
                   2263:                return R;
                   2264:        }
                   2265: }
                   2266:
                   2267: def gen_nf(F,G,V,Ord,Mod)
                   2268: {
                   2269:        if ( !Mod ) return p_nf(F,G,V,Ord);
                   2270:
                   2271:        setmod(Mod);
                   2272:        dp_ord(Ord); DF = dp_mod(dp_ptod(F,V),Mod,[]);
                   2273:        N = length(G); DG = newvect(N);
                   2274:        for ( I = N-1, IL = []; I >= 0; I-- ) {
                   2275:                DG[I] = dp_mod(dp_ptod(G[I],V),Mod,[]);
                   2276:                IL = cons(I,IL);
                   2277:        }
                   2278:        T = dp_nf_mod(IL,DF,DG,1,Mod);
                   2279:        for ( R = 0; T; T = dp_rest(T) )
                   2280:                R += gen_mptop(dp_hc(T))*dp_dtop(dp_ht(T),V);
                   2281:        return R;
                   2282: }
                   2283:
                   2284: /* Ti = [D,I,M,C] */
                   2285:
                   2286: def compute_deg0(Ti,P,V,TV)
                   2287: {
                   2288:        N = length(P[0]);
                   2289:        Num = vector(N);
                   2290:        for ( I = 0; I < N; I++ ) Num[I] = -1;
                   2291:        for ( ; Ti != []; Ti = cdr(Ti) ) {
                   2292:                Sj = car(Ti);
                   2293:                Dj = Sj[0];
                   2294:                Ij =Sj[1];
                   2295:                Mj = deg(type(Sj[2])==9?dp_dtop(Sj[2],V):Sj[2],TV);
                   2296:                Pj = P[Ij];
                   2297:                if ( Dj )
                   2298:                        for ( I = 0; I < N; I++ )
                   2299:                                if ( Pj[I] >= 0 ) {
                   2300:                                        T = Mj+Pj[I];
                   2301:                                        Num[I] = MAX(Num[I],T);
                   2302:                                }
                   2303:        }
                   2304:        return Num;
                   2305: }
                   2306:
                   2307: def compute_deg(B,V,TV,Data)
                   2308: {
                   2309:        GB = Data[0];
                   2310:        Homo = Data[1];
                   2311:        Trace = Data[2];
                   2312:        IntRed = Data[3];
                   2313:        Ind = Data[4];
                   2314:        DB = map(dp_ptod,B,V);
                   2315:        if ( Homo ) {
                   2316:                DB = map(dp_homo,DB);
                   2317:                V0 = append(V,[hhh]);
                   2318:        } else
                   2319:                V0 = V;
                   2320:        Perm = Trace[0]; Trace = cdr(Trace);
                   2321:        for ( I = length(Perm)-1, T = Trace; T != []; T = cdr(T) )
                   2322:                if ( (J=car(T)[0]) > I ) I = J;
                   2323:        N = I+1;
                   2324:        N0 = length(B);
                   2325:        P = vector(N);
                   2326:        for ( T = Perm, I = 0; T != []; T = cdr(T), I++ ) {
                   2327:                Pi = car(T);
                   2328:                C = vector(N0);
                   2329:                for ( J = 0; J < N0; J++ ) C[J] = -1;
                   2330:                C[Pi[1]] = 0;
                   2331:                P[Pi[0]] = C;
                   2332:        }
                   2333:        for ( T = Trace; T != []; T = cdr(T) ) {
                   2334:                Ti = car(T); P[Ti[0]] = compute_deg0(Ti[1],P,V0,TV);
                   2335:        }
                   2336:        M = length(Ind);
                   2337:        for ( T = IntRed; T != []; T = cdr(T) ) {
                   2338:                Ti = car(T); P[Ti[0]] = compute_deg0(Ti[1],P,V,TV);
                   2339:        }
                   2340:        R = [];
                   2341:        for ( J = 0; J < M; J++ ) {
                   2342:                U = P[Ind[J]];
                   2343:                R = cons(U,R);
                   2344:        }
                   2345:        return reverse(R);
                   2346: }
                   2347:
                   2348: /* set theoretic functions */
                   2349:
                   2350: def member(A,S)
                   2351: {
                   2352:        for ( ; S != []; S = cdr(S) )
                   2353:                if ( car(S) == A ) return 1;
                   2354:        return 0;
                   2355: }
                   2356:
                   2357: def elimination(G,V) {
                   2358:        for ( R = [], T = G; T != []; T = cdr(T) )
                   2359:                if ( setminus(vars(car(T)),V) == [] ) R =cons(car(T),R);
                   2360:        return R;
                   2361: }
                   2362:
                   2363: def setintersection(A,B)
                   2364: {
                   2365:        for ( L = []; A != []; A = cdr(A) )
                   2366:                if ( member(car(A),B) )
                   2367:                        L = cons(car(A),L);
                   2368:        return L;
                   2369: }
                   2370:
                   2371: def setminus(A,B) {
                   2372:        for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
                   2373:                for ( S = B, M = car(T); S != []; S = cdr(S) )
                   2374:                        if ( car(S) == M ) break;
                   2375:                if ( S == [] ) R = cons(M,R);
                   2376:        }
                   2377:        return R;
                   2378: }
                   2379:
                   2380: def sep_list(L,N)
                   2381: {
                   2382:        if ( length(L) <= N ) return [L,[]];
                   2383:        R = [];
                   2384:        for ( T = L, I = 0; I < N; I++, T = cdr(T) )
                   2385:                R = cons(car(T),R);
                   2386:        return [reverse(R),T];
                   2387: }
                   2388:
                   2389: def first(L)
                   2390: {
                   2391:        return L[0];
                   2392: }
                   2393:
                   2394: def second(L)
                   2395: {
                   2396:        return L[1];
                   2397: }
                   2398:
                   2399: def third(L)
                   2400: {
                   2401:        return L[2];
                   2402: }
                   2403:
                   2404: def first_second(L)
                   2405: {
                   2406:        return [L[0],L[1]];
                   2407: }
                   2408:
                   2409: def comp_tord(A,B)
                   2410: {
                   2411:        DA = dp_ht(A);
                   2412:        DB = dp_ht(B);
                   2413:        if ( DA > DB ) return 1;
                   2414:        else if ( DA < DB ) return -1;
                   2415:        else return 0;
                   2416: }
                   2417:
                   2418: def comp_tdeg(A,B)
                   2419: {
                   2420:        DA = tdeg(A);
                   2421:        DB = tdeg(B);
                   2422:        if ( DA > DB ) return 1;
                   2423:        else if ( DA < DB ) return -1;
                   2424:        else return 0;
                   2425: }
                   2426:
                   2427: def comp_tdeg_first(A,B)
                   2428: {
                   2429:        DA = tdeg(A[0]);
                   2430:        DB = tdeg(B[0]);
                   2431:        if ( DA > DB ) return 1;
                   2432:        else if ( DA < DB ) return -1;
                   2433:        else return 0;
                   2434: }
                   2435:
                   2436: def comp_third_tdeg(A,B)
                   2437: {
                   2438:        if ( A[2] > B[2] ) return 1;
                   2439:        if ( A[2] < B[2] ) return -1;
                   2440:        DA = tdeg(A[0]);
                   2441:        DB = tdeg(B[0]);
                   2442:        if ( DA > DB ) return 1;
                   2443:        else if ( DA < DB ) return -1;
                   2444:        else return 0;
                   2445: }
                   2446:
                   2447: def tdeg(P)
                   2448: {
                   2449:        dp_ord(0);
                   2450:        return dp_td(dp_ptod(P,vars(P)));
                   2451: }
                   2452:
                   2453: def comp_by_ord(A,B)
                   2454: {
                   2455:        if ( dp_ht(A) > dp_ht(B) ) return 1;
                   2456:        else if ( dp_ht(A) < dp_ht(B) ) return -1;
                   2457:        else return 0;
                   2458: }
                   2459:
                   2460: def comp_by_second(A,B)
                   2461: {
                   2462:        if ( A[1] > B[1] ) return 1;
                   2463:        else if ( A[1] < B[1] ) return -1;
                   2464:        else return 0;
                   2465: }
                   2466:
                   2467: def get_lc(F)
                   2468: {
                   2469:        if ( type(F)==1 ) return F;
                   2470:        V = var(F);
                   2471:        D = deg(F,V);
                   2472:        return get_lc(coef(F,D,V));
                   2473: }
                   2474:
                   2475: def tomonic(F,Mod)
                   2476: {
                   2477:        C = get_lc(F);
                   2478:        IC = inv(C,Mod);
                   2479:        return (IC*F)%Mod;
                   2480: }
                   2481:
                   2482: def gen_gb_comp(A,B,Mod)
                   2483: {
                   2484:        if ( !Mod ) return gb_comp(A,B);
                   2485:        LA = length(A); LB = length(B);
                   2486:        if ( LA != LB ) return 0;
                   2487:        A = map(tomonic,A,Mod);
                   2488:        B = map(tomonic,B,Mod);
                   2489:        A = qsort(A); B = qsort(B);
                   2490:        if ( A != B ) return 0;
                   2491:        return 1;
                   2492: }
                   2493:
                   2494: def prod(L)
                   2495: {
                   2496:        for ( R = 1; L != []; L = cdr(L) )
                   2497:                R *= car(L);
                   2498:        return R;
                   2499: }
                   2500:
                   2501: def monodec0(B,V)
                   2502: {
                   2503:        M = monodec(B,V);
                   2504:        return map(vars,M);
                   2505: }
                   2506:
                   2507: def monodec(B,V)
                   2508: {
                   2509:        B = map(sq,B,0);
                   2510:        G = nd_gr_postproc(B,V,0,0,0);
                   2511:        V = vars(G);
                   2512:        N = length(V);
                   2513:        if ( N == 0 ) return G == [] ? [[]] : [];
                   2514:        if ( N == 1 ) return G;
                   2515:        if ( N < 20 ) {
                   2516:                T = dp_mono_raddec(G,V);
                   2517:                return map(prod,T);
                   2518:        }
                   2519:        X = car(V); W = cdr(V);
                   2520:        D0 = monodec(map(subst,B,X,0),W);
                   2521:        T0 = map(dp_ptod,D0,W);
                   2522:        D1 = monodec(map(subst,B,X,1),W);
                   2523:        T1 = map(dp_ptod,D1,W);
                   2524: #if 0
                   2525:        for ( T = T1; T != []; T = cdr(T) ) {
                   2526:                for ( M = car(T), S1 = [], S = T0; S != []; S = cdr(S) )
                   2527:                        if ( !dp_redble(car(S),M) ) S1= cons(car(S),S1);
                   2528:                T0 = S1;
                   2529:        }
                   2530: #else
                   2531:        T0 = dp_mono_reduce(T0,T1);
                   2532: #endif
                   2533:        D0 = map(dp_dtop,T0,W);
                   2534:        D0 = vtol(X*ltov(D0));
                   2535:        return append(D0,D1);
                   2536: }
                   2537:
                   2538: def separator(P,V)
                   2539: {
                   2540:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2541:        N = length(P);
                   2542:        dp_ord(0);
                   2543:        DP = vector(N);
                   2544:        for ( I = 0; I < N; I++ ) DP[I] = qsort(ltov(map(dp_ptod,P[I][0],V)),noro_pd.comp_tord);
                   2545:        if ( Mod )
                   2546:                for ( I = 0; I < N; I++ ) DP[I] = map(dp_mod,DP[I],Mod,[]);
                   2547:        Ind = vector(N);
                   2548:        for ( I = 0; I < N; I++ ) {
                   2549:                for ( K = [], J = length(DP[I])-1; J >= 0; J-- ) K = cons(J,K);
                   2550:                Ind[I] = K;
                   2551:        }
                   2552:        S = vector(N);
                   2553:        for ( I = 0; I < N; I++ ) S[I] = 1;
                   2554:        for ( I = 0; I < N; I++ ) {
                   2555:                print(".",2);
                   2556:                for ( J = 0; J < N; J++ ) {
                   2557:                        if ( J == I ) continue;
                   2558:                        T = DP[I]; L = length(T);
                   2559:                        if ( Mod ) {
                   2560:                                for ( K = 0; K < L; K++ )
                   2561:                                        if ( dp_nf_mod(Ind[J],T[K],DP[J],0,Mod) ) break;
                   2562:                        } else {
                   2563:                                for ( K = 0; K < L; K++ )
                   2564:                                        if ( dp_nf(Ind[J],T[K],DP[J],0) ) break;
                   2565:                        }
                   2566:                        S[J] = lcm(S[J],dp_dtop(T[K],V));
                   2567:                }
                   2568:        }
                   2569:        print("");
                   2570:        return S;
                   2571: }
                   2572:
                   2573: def prepost(PL,V)
                   2574: {
                   2575:        if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
                   2576:        A = ltov(PL); N = length(A);
                   2577:        Pre = vector(N);
                   2578:        Post = vector(N);
                   2579:        R = vector(N);
                   2580:        Pre[0] = [1];
                   2581:        print("pre ",2);
                   2582:        for ( I = 1; I < N; I++, print(".",2) )
                   2583:                Pre[I] = ideal_intersection_m(Pre[I-1],A[I-1],V,0|mod=Mod);
                   2584:        print("done");
                   2585:        print("post ",2);
                   2586:        Post[N-1] = [1];
                   2587:        for ( I = N-2; I >= 0; I--, print(".",2) )
                   2588:                Post[I] = ideal_intersection_m(Post[I+1],A[I+1],V,0|mod=Mod);
                   2589:        print("done");
                   2590:        print("int ",2);
                   2591:        for ( I = 0; I < N; I++, print(".",2) )
                   2592:                R[I] = ideal_intersection_m(Pre[I],Post[I],V,0|mod=Mod);
                   2593:        print("done");
                   2594:        return R;
                   2595: }
                   2596:
                   2597: /* XXX */
                   2598:
                   2599: def call_func(Arg)
                   2600: {
                   2601:        F = car(Arg);
                   2602:        return call(strtov(F),cdr(Arg));
                   2603: }
                   2604:
                   2605: def call_func_serial(Arg,Serial)
                   2606: {
                   2607:        F = car(Arg);
                   2608:        return [call(strtov(F),cdr(Arg)),Serial];
                   2609: }
                   2610:
                   2611: def competitive_exec(P,Arg0,Arg1)
                   2612: {
                   2613:        P0 = P[0]; P1 = P[1];
                   2614:        ox_cmo_rpc(P0,"noro_pd.call_func",Arg0|sync=1);
                   2615:        ox_cmo_rpc(P1,"noro_pd.call_func",Arg1|sync=1);
                   2616:        F = ox_select(P);
                   2617:        R = ox_get(F[0]);
                   2618:        if ( length(F) == 2 ) {
                   2619:                ox_get(F[1]);
                   2620:        } else {
                   2621:                U = setminus(P,F);
                   2622:                ox_reset(U[0]);
                   2623:        }
                   2624:        return R;
                   2625: }
                   2626:
                   2627:
                   2628: def nd_gr_rat(B,V,PV,Ord1,Ord)
                   2629: {
                   2630:        G = nd_gr(B,append(V,PV),0,Ord1);
                   2631:        G1 = nd_gr_postproc(G,V,0,Ord,0);
                   2632:        return G1;
                   2633: }
                   2634:
                   2635: /* Task[i] = [fname,[arg0,...,argn]] */
                   2636:
                   2637: def para_exec(Proc,Task) {
                   2638:        Free = Proc;
                   2639:        N = length(Task);
                   2640:        R = [];
                   2641:        print([N],2); print("->",2);
                   2642:        Serial = 0;
                   2643:        while ( N ) {
                   2644:                while ( Task != [] && Free != [] ) {
                   2645:                        T = car(Task); Task = cdr(Task);
                   2646:                        ox_rpc(car(Free),"noro_pd.call_func_serial",T,Serial++);
                   2647:                        ox_push_cmd(car(Free),258); Free = cdr(Free);
                   2648:                }
                   2649:                Finish0 = Finish = ox_select(Proc);
                   2650:                for ( ; Finish != []; Finish = cdr(Finish) ) {
                   2651:                        print(".",2);
                   2652:                        L = ox_get(car(Finish));
                   2653:                        R = cons(L,R);
                   2654:                        N--;
                   2655:                }
                   2656:                print([N],2);
                   2657:                Free = append(Free,Finish0);
                   2658:        }
                   2659:        print("");
                   2660:        R = qsort(R,noro_pd.comp_by_second);
                   2661:        R = map(first,R);
                   2662:        return R;
                   2663: }
                   2664:
                   2665: def redbase(B,V,Mod,Ord)
                   2666: {
                   2667:        M = nd_gr_postproc(B,V,Mod,Ord,0);
                   2668:        dp_ord(Ord);
                   2669:        DM = ltov(map(dp_ptod,M,V));
                   2670:        if ( Mod ) DM = map(dp_mod,DM,Mod,[]);
                   2671:        N = length(DM);
                   2672:        for ( Ind = [], I = N-1; I >= 0; I-- ) Ind = cons(I,Ind);
                   2673:        for ( T = B, R = vtol(DM); T != []; T = cdr(T) ) {
                   2674:                D = dp_ptod(car(T),V);
                   2675:                if ( Mod ) D = dp_mod(D,Mod,[]);
                   2676:                D = Mod?dp_nf_mod(Ind,D,DM,1,Mod):dp_nf(Ind,D,DM,1);
                   2677:                if ( D ) R = cons(D,R);
                   2678:        }
                   2679:        D = qsort(R,noro_pd.comp_tord);
                   2680:        return map(dp_dtop,D,V);
                   2681: }
                   2682:
                   2683: def witness(A,B,V)
                   2684: {
                   2685:        G = nd_gr(B,V,0,Mod);
                   2686:        L = length(A);
                   2687:        QL = []; PL = [];
                   2688:        for ( I = L-1; I >= 0; I-- ) {
                   2689:                QL = append(map(first,A[I]),QL);
                   2690:                PL = append(map(second,A[I]),PL);
                   2691:        }
                   2692:        N = length(QL);
                   2693:        Qhat = prepost(QL,V);
                   2694:        for ( I = 0, W = []; I < N; I++ ) {
                   2695:                for ( T = Qhat[I]; T != []; T = cdr(T) )
                   2696:                        if ( gen_nf(car(T),QL[I],V,0,Mod) ) break;
                   2697:                Ai = car(T);
                   2698:                Ji = colon(G,Ai,V|isgb=1,mod=Mod);
                   2699:                Ji = nd_gr(Ji,V,Mod,0);
                   2700:                if ( gen_gb_comp(Ji,PL[I],Mod) ) Bi = 1;
                   2701:                else {
                   2702:                        Ki = ideal_colon(Ji,PL[I],V|mod=Mod);
                   2703:                        for ( T = Ki; T != []; T = cdr(T) )
                   2704:                                if ( gen_nf(car(T),Ji,V,0,Mod) ) break;
                   2705:                        Bi = car(T);
                   2706:                }
                   2707:                W = cons(Ai*Bi,W);
                   2708:                Li = colon(G,W[0],V|isgb=1,mod=Mod);
                   2709:                Li = nd_gr(Li,V,Mod,0);
                   2710:                if ( !gen_gb_comp(Li,PL[I],Mod) )
                   2711:                        error("afo");
                   2712:        }
                   2713:        return reverse(W);
                   2714: }
1.5     ! noro     2715:
        !          2716: def zero_dim(G,V,O) {
        !          2717:        dp_ord(O);
        !          2718:        HL = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
        !          2719:        for ( L = []; HL != []; HL = cdr(HL) )
        !          2720:                if ( length(vars(car(HL))) == 1 )
        !          2721:                        L = cons(car(HL),L);
        !          2722:        return length(vars(L)) == length(V) ? 1 : 0;
        !          2723: }
1.1       noro     2724: endmodule$
                   2725: end$

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