![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to primdec CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib|
Return to primdec CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 noro 1: /* $OpenXM: OpenXM/src/asir99/lib/primdec,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
2: /* Primary decomposition & Radical decomposition program */
3: /* written by T.Shimoyama, Fujitsu Lab. Date: 1995.10.12 */
4:
5: /* comments about flags *
6: PRIMEORD term order of Groebner basis in primedec
7: PRIMAORD term order of Groebner basis in primadec
8: PRINTLOG print out log (1 (simple) -- 5 (pricise))
9: SHOWTIME if 1 : print out timings and data
10: ISOLATED if 1 : compute only isolated primary components
11: NOEMBEDDED if 1 : compute only pseudo-primary components
12: NOINITGB if 1 : no initial G-base computation
13: REDUNDANT if 1 : no redundancy check
14: COMPLETE if 1 : use complete criterion of redundancy
15: COMMONCHECK if 1 : redundancy check by intersection (in gr_fctr_sub)
16: SELECTFLAG selection strategy of separators (0 -- 3)
17: */
18:
19: if (vtype(minipoly) != 3) load("gr")$$
20:
21: #define GR(R,F,V,O) T2=newvect(4,time());R=dp_gr_main(F,V,0,0,O);GRTIME+=newvect(4,time())-T2;
22: #define HGRM(R,F,V,O) T2=newvect(4,time());R=dp_gr_main(F,V,1,1,O);GRTIME+=newvect(4,time())-T2;
23: #define NF(R,IN,F,G,O) T2=newvect(4,time());R=dp_nf(IN,F,G,O);NFTIME+=newvect(4,time())-T2;
24:
25: /* optional flags */
26: extern PRIMAORD,PRIMEORD,PRINTLOG,SHOWTIME,NOINITGB$
27: extern COMMONCHECK,ISOLATED,NOEMBEDDED,REDUNDANT,SELECTFLAG,COMPLETE$
28: extern EXTENDED,CONTINUE,BSAVE,MAXALGIND,NOSPECIALDEC$
29: def primflags() {
30: print("PRIMAORD,PRIMEORD,PRINTLOG,SHOWTIME,NOINITGB,ISOLATED,NOEMBEDDED,COMMONCHECK");
31: print("REDUNDANT,SELECTFLAG,COMPLETE,EXTENDED,CONTINUE,BSAVE,MAXALGIND,NOSPECIALDEC");}
32: PRIMAORD=0$ PRIMEORD=3$
33:
34: /* global variables */
35: extern DIRECTRY,COMMONIDEAL,NOMRADDEC,NOMISOLATE,RADDECTIME$
36: extern MISTIME,ISORADDEC,GRTIME,NFTIME$
37:
38:
39: /* primary ideal decomposition of ideal(F) */
40: /* return the list of [P,Q] where Q is P-primary component of ideal(F) */
41: def primadec(F,VL)
42: {
43: if ( !VL ) {
44: print("invalid variables"); return 0; }
45: if ( !F ) {
46: print("invalid argument"); return 0; }
47: if ( F == [] )
48: return [];
49: if ( length(F) == 1 && type(F[0]) == 1 )
50: return [[1],[1]];
51: NOMRADDEC = NOMISOLATE = 0; RADDECTIME = newvect(4); T0 = newvect(4,time());
52: GRTIME = newvect(4); NFTIME = newvect(4); MISTIME = newvect(4);
53:
54: R = primadec_main(F,[["begin",F,F,[1],[]]],[],VL);
55:
56: if ( PRINTLOG ) {
57: print(""); print("primary decomposition done.");
58: }
59: if ( PRINTLOG || SHOWTIME ) {
60: print("number of radical decompositions = ",0); print(NOMRADDEC);
61: print("number of primary components = ",0); print(length(R),0);
62: print(" ( isolated = ",0); print(NOMISOLATE,0); print(" )");
63: print("total time of m.i.s. computation = ",0);
64: print(MISTIME[0],0); GT = MISTIME[1];
65: if ( GT ) { print(" + gc ",0);print(GT); } else print("");
66: print("total time of n-form computation = ",0);
67: print(NFTIME[0],0); GT = NFTIME[1];
68: if ( GT ) { print(" + gc ",0);print(GT); } else print("");
69: print("total time of G-base computation = ",0);
70: print(GRTIME[0],0); GT = GRTIME[1];
71: if ( GT ) { print(" + gc ",0);print(GT); } else print("");
72: print("total time of radical decomposition = ",0);
73: print(RADDECTIME[0],0); GT = RADDECTIME[1];
74: if ( GT ) { print(" + gc ",0);print(GT,0); }
75: print(" (iso ",0); print(ISORADDEC[0],0); GT = ISORADDEC[1];
76: if ( GT ) { print(" +gc ",0);print(GT,0); } print(")");
77: print("total time of primary decomposition = ",0);
78: TT = newvect(4,time())-T0; print(TT[0],0);
79: if ( TT[1] ) { print(" + gc ",0);print(TT[1]); } else print("");
80: }
81: return R;
82: }
83:
84: /* prime ideal decomposition of radical(F) */
85: /* return the list of prime components of radical(F) */
86:
87: def primedec(F,VL)
88: {
89: if ( !VL ) {
90: print("invalid variables"); return 0; }
91: if ( !F ) {
92: print("invalid argument"); return 0; }
93: if ( F == [] )
94: return [];
95: GRTIME = newvect(4);
96: T0 = newvect(4,time());
97: if ( !NOINITGB ) {
98: if ( PRINTLOG ) {
99: print("G-base computation w.r.t ordering = ",0);
100: print(PRIMEORD);
101: }
102: T1 = newvect(4,time());
103: HGRM(F,F,VL,PRIMEORD);
104: Tg = newvect(4,time())-T1;
105: if ( PRINTLOG > 0 ) {
106: print(" G-base time = ",0); print(Tg[0],0);
107: if ( Tg[1] ) { print(" + gc : ",0); print(Tg[1]); }
108: else print("");
109: }
110: }
111:
112: R = primedec_main([F],VL);
113:
114: Ta = newvect(4,time())-T0;
115: if ( PRINTLOG || SHOWTIME ) {
116: print("number of prime components = ",0); print(length(R));
117: print("G-base time = ",0); print(GRTIME[0],0);
118: if ( GRTIME[1] ) { print(" + gc : ",0); print(GRTIME[1]); }
119: else print("");
120: print("total time = ",0); print(Ta[0],0);
121: if ( Ta[1] ) { print(" + gc : ",0); print(Ta[1]); } else print("");
122: }
123: return R;
124: }
125:
126: /* main procedure for primary decomposition. */
127: /* F : ideal, VL : variable list, REMS : remaining components */
128: def primadec_main(F,REMS,H,VL)
129: {
130: DEC = RES = [];
131: for (D = [], I = 0; I < length(REMS); I++) {
132: MARK = REMS[I][0]; WF = REMS[I][1]; WP = REMS[I][2]; SC = REMS[I][3];
133: if ( !NOINITGB || MARK != "begin" ) {
134: ORD = PRIMAORD;
135: if ( PRINTLOG > 1 ) {
136: if ( MARK != "begin" ) print("");
137: print("G-base computation w.r.t ordering = ",0);
138: print(ORD);
139: }
140: T1 = newvect(4,time());
141: /* G-base of ideal */
142: GR(GF,WF,VL,ORD);
143: if ( MARK != "begin" && ( COMPLETE || EXTENDED ) ) {
144: if ( SC!=[1] && SC!=[-1] ) {
145: LG = localize(WF,SC,VL,ORD); /* VR_s\cap R */
146: if ( EXTENDED ) { GF = LG; }
147: } else
148: LG = GF;
149: if ( idealinc(H,LG,VL,ORD) ) {
150: if ( PRINTLOG ) {
151: print("eliminate the remaining component ",0);
152: print("by complete criterion");
153: }
154: continue; /* complete criterion */
155: }
156: }
157: /* G-base of radical */
158: if ( MARK == "begin" ) {
159: RA = ["begin",GF];
160: } else if ( MARK == "ext" ) {
161: if ( WF==WP || idealinc(WP,GF,VL,ORD) )
162: RA = ["ext",GF];
163: else {
164: if ( EXTENDED ) {
165: GA = localize(WP,SC,VL,PRIMEORD);
166: } else {
167: GR(GA,WP,VL,PRIMEORD);
168: }
169: RA = ["ext",GA];
170: }
171: } else if ( MARK == "sep" ) {
172: for (U=[], T=WP,J=length(T)-1;J>=0;J--) {
173: if ( EXTENDED ) {
174: GA = localize(T[J],SC,VL,PRIMEORD);
175: } else {
176: GR(GA,T[J],VL,PRIMEORD);
177: }
178: if (GA != [1] && GA != [-1])
179: U = cons(GA,U);
180: }
181: RA = ["sep",U];
182: } else debug;
183: Tg = newvect(4,time())-T1;
184: if ( PRINTLOG > 1 ) {
185: print(" G-base time = ",0); print(Tg[0],0);
186: if ( Tg[1] ) { print(" + gc : ",0); print(Tg[1]); }
187: else print("");
188: }
189: } else {
190: GF = F; /* NOINITGB = 1 */
191: RA = ["begin",GF];
192: }
193: if ( PRINTLOG ) {
194: if ( MARK == "begin" ) {
195: print("primary decomposition of the ideal");
196: } else { /* at the begining */
197: print("");
198: print("primary decomposition of the remaining component");
199: }
200: if ( MARK == "begin" && PRINTLOG > 1 ) { /* at the begining */
201: print(" ideal = ",0);
202: print(WF);
203: } else if ( PRINTLOG >= 4 ) {
204: print(" remaining component = ",0);
205: print(GF);
206: }
207: }
208:
209: /* input: init, generator, G-base, radical, intersection, sep.cond.*/
210: /* output: primary comp. remaining comp. */
211: Y = isolated(F,WF,GF,RA,REMS[I][4],SC,VL);
212:
213: D = append(D,Y[0]);
214: if ( MARK == "begin" )
215: NOMISOLATE = length(Y[0]);
216: RES = append(RES,Y[1]);
217: }
218: if ( MARK == "begin" ) {
219: F = GF; /* input polynomial -> G-base of it */
220: }
221: DEC = append(DEC,D);
222: if ( PRINTLOG ) {
223: print("");
224: print("# number of remainig components = ",0); print(length(RES));
225: }
226: if ( !length(RES) )
227: return DEC;
228: if ( !REDUNDANT ) { /* check whether Id(F,RES) is redundant or not */
229: for(L = [H],I = length(D)-1; I>=0; I--)
230: L=append(L,[D[I][0]]);
231: H = idealsetcap(L,VL,ORD); /* new intersection H */
232: if ( idealinc(H,F,VL,ORD) ) {
233: if ( PRINTLOG ) {
234: print("");
235: print("all primary components are found!");
236: }
237: return DEC;
238: }
239: REMS = mksepext(RES,H,VL); /* new remainig comps. */
240: } else {
241: REMS = mksepext(RES,H,VL); /* new rmaining comps. */
242: }
243: D = primadec_main(F,REMS,H,VL);
244: return append(DEC,D);
245: }
246:
247: /* isolated components and embedded components */
248: /* GF is the G-base of F, RA is the radical of F */
249: def isolated(IP,F,GF,RA,H,SC,VL)
250: {
251: T0 = newvect(4,time());
252: if ( RA[0] == "begin" )
253: PD = primedec_main([RA[1]],VL);
254: else if ( RA[0] == "ext" || RA[0] == "sep" ) {
255: if ( RA[0] == "sep" )
256: T = prime_irred(idealsav(RA[1]),VL);
257: else
258: T = [RA[1]];
259: if ( !NOSPECIALDEC )
260: PD = primedec_special(T,VL);
261: else
262: PD = primedec_main(T,VL);
263: } else debug;
264: T1 = newvect(4,time())-T0;
265: if ( RA[0] == "begin" ) ISORADDEC = T1;
266: NOMRADDEC++; RADDECTIME += T1;
267: if ( PRINTLOG ) {
268: print("number of prime components = ",0); print(length(PD),0);
269: print(", time = ",0); print(T1[0],0);
270: if ( T1[1] ) { print(" + gc : ",0); print(T1[1]); } else print("");
271: if ( PRINTLOG > 0 ) {
272: print("dimensions of primes = ",0);
273: for (I = 0; I < length(PD); I++) {
274: print(length(VL)-length(minalgdep(PD[I],VL,PRIMEORD)),0);
275: print(", ",0);
276: }
277: print("");
278: }
279: }
280: if ( RA[0] == "begin" ) { /* isolated part of initial ideal */
281: if ( PRINTLOG ) {
282: print("check 'prime decomposition = primary decomposition?'");
283: }
284: CP = idealsetcap(PD,VL,PRIMAORD);
285: if ( idealinc(CP,GF,VL,PRIMAORD) ) {
286: if ( PRINTLOG ) {
287: print("lucky!");
288: }
289: for (L = [],I = length(PD)-1; I >= 0; I--)
290: L = cons([PD[I],PD[I]],L);
291: return [L,[]];
292: }
293: if ( PRINTLOG ) {
294: print("done.");
295: }
296: }
297:
298: R = pseudo_extract(IP,F,GF,PD,H,SC,VL);
299:
300: return R;
301: }
302:
303: def pseudo_extract(IP,F,GF,PD,H,SC,VL)
304: {
305: REMS = DEC = PDEC = SEL = RES = [];
306: ZERODIM = 1; CAP = H;
307: for (I = 0; I < length(PD); I++) {
308: P = PD[I];
309: if ( length(minalgdep(P,VL,PRIMEORD)) != length(VL) )
310: ZERODIM=0;
311: if ( length(PD) == 1 ) { /* pseudo primary ideal case */
312: if ( PRINTLOG >= 1 ) { print(""); print("pseudo primary ideal"); }
313: DD = GF; SEL = [SC];
314: } else {
315: T0 = time();
316: Y = pseudodec_main(F,P,PD,VL);
317: T1 = time();
318: DD = Y[0]; SS = Y[1]; SEL = append(SEL,[SS]);
319: PDEC = append(PDEC,[[DD,P]]);
320: if ( PRINTLOG >= 1 ) {
321: print(""); print("pseudo primary component of ",0);
322: print(I,0); print(", time = ",0); print(T1[0]-T0[0]);
323: if ( PRINTLOG >= 4 ) { print(" = ",0); print(DD); }
324: }
325: if ( NOEMBEDDED )
326: continue;
327: }
328: if ( !REDUNDANT && H != [] ) { /* elim. pseu-comp. */
329: if ( !sepcond_ps(P,SC,VL) )
330: continue;
331: LG = localize(DD,SC,VL,PRIMAORD);
332: if ( idealinc(H,LG,VL,PRIMAORD)) {
333: if ( PRINTLOG ) {
334: print("eliminate the pseudo primary component ",0);
335: print(I);
336: }
337: continue;
338: }
339: }
340: if ( PRINTLOG ) {
341: print("extraction of the pseudo primary component ",0);
342: print(I);
343: if ( PRINTLOG > 2 ) {
344: print(" associated prime of pseudo primary ideal = ",0);
345: print(P);
346: }
347: }
348: U = extraction(DD,P,VL);
349: if ( !REDUNDANT && H != [] && idealinc(H,U[0][0],VL,PRIMAORD) ) {
350: if ( PRINTLOG )
351: print("redundant primary component!");
352: } else {
353: DEC = append(DEC,[U[0]]);
354: H = idealcap(H,U[0][0],VL,PRIMAORD);
355: if ( idealeq(IP,H) ) {
356: if ( PRINTLOG ) {
357: print("");
358: print("all primary components are found!");
359: }
360: return [DEC,[]];
361: }
362: }
363: if ( !ISOLATED && U[1] != [] )
364: if ( sepcond_re(U[1],SC,VL) ) {
365: NSC = setunion([SS,SC]); /* new separating condition */
366: REM = cons(DD,append(U[1],[NSC,H]));
367: REMS = append(REMS,[REM]);
368: }
369: }
370: if ( NOEMBEDDED )
371: DEC = PDEC;
372: if ( length(PD) != 1 && !NOEMBEDDED && !ISOLATED && !ZERODIM ) {
373: for (QD=[],I=length(PDEC)-1;I>=0;I--)
374: QD = cons(PDEC[I][0],QD);
375: RES = ["sep",PD,QD,GF];
376: if ( sepcond_re(append(RES,[SEL]),SC,VL) ) {
377: REM = cons(F,append(RES,[SC,H]));
378: REMS = append(REMS,[REM]);
379: }
380: }
381: return [DEC,REMS];
382: }
383:
384: /* F:input set, PD:primes of radical, E:higher dimensional ideal */
385: /* length(PD) > 1 */
386: /* output : pseudo primary comp., remaining comp., selectors */
387: def pseudodec_main(F,P,PD,VL)
388: {
389: ZERODIM = 1;
390: S = selects(P,PD,VL,SELECTFLAG);
391: R = localize(F,S,VL,PRIMAORD);
392: if ( R[0] == 1 || R[0] == -1 ) {
393: print("error, R = ",0); print(R);
394: debug;
395: }
396: R = idealnormal(R);
397: return [R,S];
398: }
399:
400: /* Id(GF) is a pseudo primary ideal. (GF must be the G-base.) */
401: def extraction(GF,Pr,VL)
402: {
403: TMPORD1=TMPORD2=0;
404: V = minalgdep(Pr,VL,PRIMEORD);
405: U = listminus(VL,V);
406: V0 = append(V,U);
407: if ( V0 != VL ) {
408: ORD = [[TMPORD1,length(V)],[TMPORD2,length(U)]];
409: GR(G,GF,V0,ORD);
410: } else
411: G = GF;
412: dp_ord(TMPORD1);
413: for (LL = [],HC = 1,I = 0; I < length(G); I++)
414: LL = append(LL,cdr(fctr(dp_hc(dp_ptod(G[I],V)))));
415: for (L=[],I=0;I<length(LL);I++)
416: L = cons(LL[I][0],L);
417: L = setunion([L]);
418: for (S=1,SL=[],I=0;I<length(L);I++) {
419: S *= L[I];
420: if ( SELECTFLAG )
421: SL = cons(L[I],SL);
422: }
423: if ( !SELECTFLAG )
424: SL= [S];
425: if ( PRINTLOG > 1 ) {
426: print("extractor = ",0);
427: print(S);
428: }
429: T0 = time()[0];
430: Q = localize(GF,SL,VL,PRIMAORD);
431: Q = idealnormal(Q);
432: DEC = [Q,Pr];
433: if ( PRINTLOG ) {
434: print("find a primary component! time = ",0);
435: print(time()[0]-T0);
436: if (PRINTLOG >= 3){
437: print(" associated prime of primary component = ",0);
438: print(DEC[1]);
439: }
440: if (PRINTLOG >= 4){print(" primary component = ",0); print(Q);}
441: }
442: if ( !ISOLATED && !NOEMBEDDED && SL != [1]
443: && length(V) != length(VL) /* nonzerodim */
444: && (REDUNDANT || !idealinc(Q,GF,VL,PRIMAORD)) ) {
445: REM = ["ext",[Q,Pr],GF,S];
446: if ( PRINTLOG ) {
447: print("find the remaining component of the extraction");
448: }
449: } else {
450: REM = [];
451: if ( PRINTLOG ) {
452: print("eliminate the remaining component of the extraction");
453: }
454: }
455: return [DEC,REM];
456: }
457:
458: /* sep. cond. for pseudo-primary components */
459: def sepcond_ps(P,SC,VL)
460: {
461: for (J = 0; J < length(SC); J++) {
462: if ( idealinc([SC[J]],P,VL,PRIMEORD) )
463: break; /* separating condition */
464: }
465: if ( J != length(SC) ) {
466: if ( PRINTLOG ) {
467: print("");
468: print("elim. the pseudo primary comp. by separating cond.");
469: }
470: return 0;
471: }
472: return 1;
473: }
474:
475: /* sep. cond. for rem. components. */
476: /* REM = ["ext",[Q,Pr],GF,S] or ["sep",PD,QD,GF,SEL], SC : sep.cond. */
477: def sepcond_re(REM,SC,VL)
478: {
479: for (S=1,I=0;I<length(SC);I++)
480: S *= SC[I];
481: if (REM[0] == "ext") {
482: F = cons(REM[3],REM[1][1]);
483: L = localize(F,[S],VL,PRIMAORD);
484: if ( L != [1] )
485: return 1;
486: else
487: return 0;
488: } else if (REM[0] == "sep") {
489: PL = REM[1]; SEL = REM[4];
490: for (I=0;I<length(PL);I++) {
491: F = append(SEL[I],PL[I]);
492: L = localize(F,[S],VL,PRIMAORD);
493: if ( L != [1] )
494: return 1;
495: }
496: return 0;
497: }
498: }
499:
500: def minalgdep(Pr,VL,ORD)
501: {
502: T0=newvect(4,time());
503: if (MAXALGIND)
504: R = minalgdep1(Pr,VL,ORD); /* M.I.S. */
505: else
506: R = minalgdep0(Pr,VL,ORD); /* M.S.I.S. */
507: T1=newvect(4,time())-T0; MISTIME += T1;
508: if ( PRINTLOG >= 6 ) {
509: if (MAXALGIND) print("M.I.S. time = ",0);
510: else print("M.S.I.S. time = ",0);
511: print(T1[0]);
512: }
513: return R;
514: }
515:
516: def minalgdep0(Pr,VL,ORD)
517: {
518: dp_ord(ORD); N = length(Pr);
519: for (HD = [],I = N-1; I >= 0; I--)
520: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
521: for (R = X = [], I = length(VL)-1; I >= 0; I--) {
522: V = VL[I];
523: TX = cons(V,X);
524: for (J = 0; J < N; J++) {
525: if ( varinc(HD[J],TX) )
526: break;
527: }
528: if ( J == N )
529: X = TX;
530: else
531: R = cons(V,R);
532: }
533: return R;
534: }
535:
536: def minalgdep1(Pr,VL,ORD)
537: {
538: R = all_msis(Pr,VL,ORD);
539: for (E=I=0,D=length(R[0]);I<length(R);I++) {
540: if ( length(R[I])>=D ) {
541: D=length(R[I]); E=I;
542: }
543: }
544: S = listminus(VL,R[E]);
545: return S;
546: }
547:
548: def all_msis(Pr,VL,ORD)
549: {
550: dp_ord(ORD); N = length(Pr);
551: for (HD = [],I = N-1; I >= 0; I--)
552: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
553: R = dimrec(HD,0,[],[],VL);
554: return R;
555: }
556:
557: def dimrec(S,K,U,M,X)
558: {
559: MM = M;
560: for (I=K;I<length(X);I++) {
561: TX = append(U,[X[I]]);
562: for (J = 0; J < length(S); J++) {
563: if ( varinc(S[J],TX) )
564: break;
565: }
566: if ( J == length(S) )
567: MM = dimrec(S,I+1,TX,MM,X);
568: }
569: for (J = 0; J < length(MM); J++) {
570: if ( varinc(U,MM[J]) )
571: break;
572: }
573: if ( J == length(MM) )
574: MM = append(MM,[U]);
575: return MM;
576: }
577:
578: /* if V is min.alg.dep. return 1 */
579: def ismad(Pr,V,VL,ORD)
580: {
581: dp_ord(ORD); N = length(Pr);
582: for (HD = [],I = N-1; I >= 0; I--)
583: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
584: W = listminus(VL,V); L = append(W,V);
585: for (R = TX = [],I = 0; I < length(L); I++) {
586: TX = cons(L[I],TX);
587: for (J = 0; J < N; J++) {
588: if ( varinc(HD[J],TX) )
589: break;
590: }
591: if ((I<length(W)&&J!=N) || (I>=length(W)&&J==N))
592: return 0;
593:
594: }
595: return 1;
596: }
597:
598: /* output:list of [flag, generator, primes, sep.cond., intersection] */
599: def mksepext(RES,H,VL)
600: {
601: for (R = [],I=length(RES)-1;I>=0;I--) {
602: W = RES[I];
603: if ( COMPLETE || EXTENDED )
604: HI = idealcap(H,W[6],VL,PRIMAORD);
605: else
606: HI = H;
607: if (W[1] == "sep") {
608: E = mksep(W[2],W[3],W[4],VL);
609: F = append(E[0],W[0]);
610: if ( ( COMPLETE || EXTENDED ) && !REDUNDANT )
611: HI = idealsetcap(cons(HI,W[3]),VL,PRIMAORD);
612: R = cons(["sep",F,E[1],W[5],HI],R);
613: } else if (W[1] == "ext") {
614: E = mkext(W[2][0],W[3],W[4],VL);
615: F = cons(E[0],W[0]);
616: P = cons(E[1],W[2][1]); /* prime component of rem. comp. */
617: R = cons(["ext",F,P,W[5],HI],R);
618: } else debug;
619: }
620: return R;
621: }
622:
623: /* make extractor with exponents ; F must be G-Base */
624: /* output : extractor, max.sqfr.of extractor */
625: def mkext(Q,F,S,VL)
626: {
627: if ( PRINTLOG > 1 ) {
628: print("make extractor = ",0);
629: }
630: if ( S == 1 || S == -1 ) {
631: print(1); return [1,1];
632: }
633: DD=dp_mkdist([Q],F,VL,PRIMAORD);
634: DQ=DD[0][0]; DF=DD[1];
635: for (L=fctr(S),SL=[],J=length(L)-1;J>0;J--)
636: SL=cons(dp_ptod(L[J][0],VL),SL);
637: K = dp_fexponent(dp_ptod(S,VL),DF,DQ);
638: E = dp_vectexp(SL,EX=dp_minexp(K,SL,DF,DQ));
639: if ( E != 1 ) E = dp_dtop(E,VL);
640: for (U=1, I=0;I<size(EX)[0];I++)
641: if ( EX[I] )
642: U *= L[I+1][0];
643: /* for (L=sqfr(E),U=1,J=length(L)-1;J>0;J--)
644: U *= L[J][0]; */
645: if ( PRINTLOG > 1 ) {
646: print(E);
647: }
648: return [E,U];
649: }
650:
651: /* make separators with exponents ; F must be G-Base */
652: /* output [separator, list of components of radical] */
653: def mksep(PP,QQ,F,VL)
654: {
655: if ( PRINTLOG > 1 ) {
656: print("make separators, ",0); print(length(PP));
657: }
658: DD=dp_mkdist(QQ,F,VL,PRIMAORD);
659: DQ=DD[0]; DF=DD[1];
660: E = dp_mksep(PP,DQ,DF,VL,PRIMAORD);
661: for (RA=[],I=length(E)-1;I>=0;I--) {
662: for (L=fctr(E[I]),J=length(L)-1;J>0;J--) {
663: ES=cons(L[J][0],PP[I]);
664: RA = cons(ES,RA); /* radical of remaining comp. */
665: }
666: }
667: return [E,RA];
668: }
669:
670: def dp_mksep(PP,DQ,DF,VL,ORD)
671: {
672: dp_ord(ORD);
673: for (E=[],I=0;I<length(PP);I++) {
674: T0=time()[0];
675: if ( CONTINUE ) {
676: E = bload("sepa");
677: I = length(E);
678: CONTINUE = 0;
679: }
680: S = selects(PP[I],PP,VL,0)[0];
681: for (L=fctr(S),SL=[],J=length(L)-1;J>0;J--)
682: SL=cons(dp_ptod(L[J][0],VL),SL);
683: K = dp_fexponent(dp_ptod(S,VL),DF,DQ[I]);
684: M = dp_vectexp(SL,dp_minexp(K,SL,DF,DQ[I]));
685: if (M != 1) M = dp_dtop(M,VL);
686: E = append(E,[M]);
687: if ( PRINTLOG > 1 ) {
688: print("separator of ",0); print(I,0);
689: print(" in ",0); print(length(PP),0);
690: print(", time = ",0); print(time()[0]-T0,0); print(", ",0);
691: printsep(cdr(fctr(M))); print("");
692: }
693: if (BSAVE) {
694: bsave(E,"sepa");
695: print("bsave to sepa");
696: }
697: }
698: return E;
699: }
700:
701: extern AFO$
702:
703: /* F,Id,Q is given by distributed polynomials. Id is a G-Base.*/
704: /* get the exponent K s.t. (Id : F^K) = Q, Id is a G-Basis */
705: def dp_fexponent(F,Id,Q)
706: {
707: if (!AFO)
708: for (IN=[],I=size(Id)[0]-1;I>=0;I--) IN = cons(I,IN);
709: else
710: for (IN=[],I=0;I<size(Id)[0];I++) IN = cons(I,IN);
711: NF(G,IN,F,Id,0);
712: if ( F != G ) { F = G; }
713: N=size(Q)[0]; U = newvect(N);
714: for (I=0;I<N;I++)
715: U[I]=Q[I];
716: for (K = 1;; K++) { /* U[I] = Q[I]*F^K mod Id */
717: for (FLAG = I = 0; I < N; I++) {
718: NF(U[I],IN,U[I]*F,Id,1);
719: if ( U[I] ) FLAG = 1;
720: }
721: if ( !FLAG )
722: return K;
723: }
724: }
725:
726: /* F must be a G-Base. SL is a set of factors of S.*/
727: def dp_minexp(K,SL,F,Q)
728: {
729: if (!AFO)
730: for (IN=[],I=size(F)[0]-1;I>=0;I--) IN = cons(I,IN);
731: else
732: for (IN=[],I=0;I<size(F)[0];I++) IN = cons(I,IN);
733: N = length(SL); M = size(Q)[0];
734: EX = newvect(N); U = newvect(M); T = newvect(M);
735: for (I=0;I<N;I++)
736: EX[I]=K;
737: for (I=0;I<M;I++) {
738: NF(Q[I],IN,Q[I],F,1);
739: }
740: for (I=0;I<N;I++) {
741: EX[I] = 0;
742: for (FF=1,II=0;II<N;II++)
743: for (J = 0; J < EX[II]; J++) {
744: NF(FF,IN,FF*SL[II],F,1);
745: }
746: for (J = 0; J < M; J++)
747: U[J] = Q[J]*FF;
748: for (; EX[I] < K; EX[I]++) {
749: FLAG = 0;
750: for (J = 0; J < M; J++) {
751: NF(U[J],IN,U[J],F,1);
752: if ( U[J] ) FLAG = 1;
753: }
754: if ( !FLAG )
755: break;
756: if ( EX[I] < K-1 )
757: for (J = 0; J < M; J++) U[J] = U[J]*SL[I];
758: }
759: }
760: return EX;
761: }
762:
763: def dp_vectexp(SS,EE)
764: {
765: N = length(SS);
766: for (A=1,I=0;I<N;I++)
767: for (J = 0; J < EE[I]; J++)
768: A *= SS[I];
769: return A;
770: }
771:
772: def dp_mkdist(QQ,F,VL,ORD)
773: {
774: dp_ord(ORD);
775: for (DQ=[],I=length(QQ)-1;I>=0;I--) {
776: for (T=[],J=length(QQ[I])-1;J>=0;J--)
777: T = cons(dp_ptod(QQ[I][J],VL),T);
778: DQ=cons(newvect(length(T),T),DQ);
779: }
780: for (T=[],J=length(F)-1;J>=0;J--)
781: T = cons(dp_ptod(F[J],VL),T);
782: DF = newvect(length(T),T);
783: return [DQ,DF]$
784: }
785:
786: /* if FLAG == 0 return singltonset */
787: def selects(P,PD,VL,FLAG)
788: {
789: PP=dp_mkdist([],P,VL,PRIMEORD)[1];
790: for (M=[],I=length(P)-1;I>=0;I--)
791: M = cons(I,M);
792: for (SL = [],I = length(PD)-1; I >= 0; I--) {
793: B = PD[I];
794: if (B == P) continue;
795: for (L = [],J = 0; J < length(B); J++) {
796: A = B[J];
797: NF(E,M,dp_ptod(A,VL),PP,0);
798: if ( E ) {
799: L = append(L,[A]); /* A \notin Id(P) */
800: }
801: }
802: SL = cons(L,SL); /* candidate of separators */
803: }
804: for (S=[],I=length(SL)-1;I>=0;I--) {/* choose an element from L */
805: if ( FLAG >= 2 ) {
806: S = append(SL[I],S);
807: continue;
808: }
809: C = minselect(SL[I],VL,PRIMAORD);
810: S = cons(C,S);
811: }
812: S = setunion([S]);
813: if ( FLAG == 3 || length(S) == 1 )
814: return S;
815: if ( FLAG == 2 ) {
816: for (R=1,I=0;I<length(S);I++)
817: R *= S[I];
818: return [R];
819: }
820: /* FLAG == 0 || FLAG == 1 */
821: for (T=[]; S!=[];S = cdr(S)) { /* FLAG <= 1 and length(S) != 1 */
822: A = car(S);
823: U = append(T,cdr(S));
824: for (R=1,I=0;I<length(U);I++)
825: R *= U[I];
826: for (I=0;I<length(PD);I++) {
827: if ( PD[I] != P && !idealinc([R],PD[I],VL,PRIMAORD) )
828: break;
829: }
830: if ( I != length(PD) )
831: T = append(T,[A]);
832: }
833: if ( FLAG )
834: return T;
835: for (R=1,I=0;I<length(T);I++)
836: R *= T[I];
837: return [R]; /* if FLAG == 0 */
838: }
839:
840: def minselect(L,VL,ORD)
841: {
842: dp_ord(ORD);
843: F = dp_ptod(L[0],VL); MAG=dp_mag(F); DEG=dp_td(F);
844: for (J = 0, I = 1; I < length(L); I++) {
845: A = dp_ptod(L[I],VL); M=dp_mag(A); D=dp_td(A);
846: if ( dp_ht(A) > dp_ht(F) )
847: continue;
848: else if ( dp_ht(A) == dp_ht(F) && (D > DEG || (D == DEG && M > MAG)))
849: continue;
850: F = A; J = I; MAG = M; DEG=D;
851: }
852: return L[J];
853: }
854:
855: /* localization Id(F)_MI \cap Q[VL] */
856: /* output is the G-base w.r.t ordering O */
857: def localize(F,MI,VL,O)
858: {
859: if ( MI == [1] || MI == [-1] )
860: return F;
861: for (SVL = [],R = [],I = 0; I < length(MI); I++) {
862: V = strtov("zz"+rtostr(I));
863: R = append(R,[V*MI[I]-1]);
864: SVL = append(SVL,[V]);
865: }
866: if ( O == 0 ) {
867: dp_nelim(length(SVL)); ORD = 6;
868: } else if ( O == 1 || O == 2 ) {
869: ORD = [[0,length(SVL)],[O,length(VL)]];
870: } else if ( O == 3 || O == 4 || O == 5 )
871: ORD = [[0,length(SVL)],[O-3,length(VL)-1],[2,1]];
872: else
873: error("invarid ordering");
874: GR(G,append(F,R),append(SVL,VL),ORD);
875: S = varminus(G,SVL);
876: return S;
877: }
878:
879: def varminus(G,VL)
880: {
881: for (S = [],I = 0; I < length(G); I++) {
882: V = vars(G[I]);
883: if (listminus(V,VL) == V)
884: S = append(S,[G[I]]);
885: }
886: return S;
887: }
888:
889: def idealnormal(L)
890: {
891: R=[];
892: for (I=length(L)-1;I>=0;I--) {
893: A = ptozp(L[I]);
894: V = vars(A);
895: for (B = A,J = 0;J < length(V);J++)
896: B = coef(B,deg(B,V[J]),V[J]);
897: R = cons((B>0?A:-A),R);
898: }
899: return R;
900: }
901:
902: def idealsav(L) /* sub procedure of welldec and normposdec */
903: {
904: if ( PRINTLOG >= 4 )
905: print("multiple check ",2);
906: for (R = [],I = 0,N=length(L); I < N; I++) {
907: if ( PRINTLOG >= 4 && !irem(I,10) )
908: print(".",2);
909: if ( R == [] )
910: R = append(R,[L[I]]);
911: else {
912: for (J = 0; J < length(R); J++)
913: if ( idealeq(L[I],R[J]) )
914: break;
915: if ( J == length(R) )
916: R = append(R,[L[I]]);
917: }
918: }
919: if ( PRINTLOG >= 4 )
920: print("done.");
921: return R;
922: }
923:
924: /* intersection of ideals in PL, PL : list of ideals */
925: /* VL : variable list, O : output the G-base w.r.t order O */
926: def idealsetcap(PL,VL,O)
927: {
928: for (U = [], I = 0; I < length(PL); I++)
929: if ( PL[I] != [] )
930: U = append(U,[PL[I]]);
931: if ( U == [] )
932: return [];
933: if ( PRINTLOG >= 4 ) {print("intersection of ideal ",0); print(length(U),0);}
934: for (L = U[0],I=1;I<length(U);I++) {
935: if ( PRINTLOG >=4 ) print(".",2);
936: L = idealcap(L,U[I],VL,O);
937: }
938: if ( PRINTLOG >=4 ) print("");
939: return L;
940: }
941:
942: /* return intersection set of ideals P and Q */
943: def idealcap(P,Q,VL,O)
944: {
945: if ( P == [] )
946: if ( VL )
947: return Q;
948: else
949: return [];
950: if ( Q == [] )
951: return P;
952: T = tmp;
953: for (PP = [],I = 0; I < length(P); I++)
954: PP = append(PP,[P[I]*T]);
955: for (QQ = [],I = 0; I < length(Q); I++)
956: QQ = append(QQ,[Q[I]*(T-1)]);
957: if ( O == 0 ) {
958: dp_nelim(1); ORD = 6;
959: } else if ( O == 1 || O == 2 )
960: ORD = [[2,1],[O,length(VL)]];
961: else if ( O == 3 || O == 4 || O == 5 )
962: ORD = [[2,1],[O-3,length(VL)-1],[2,1]];
963: else
964: error("invarid ordering");
965: GR(G,append(PP,QQ),append([T],VL),ORD);
966: S = varminus(G,[T]);
967: return S;
968: }
969:
970: /* return ideal P : Q */
971: def idealquo(P,Q,VL,O)
972: {
973: for (R = [], I = 0; I < length(Q); I++) {
974: F = car(Q);
975: G = idealcap(P,[F],VL,O);
976: for (H = [],J = 0; J < length(G); J++) {
977: H = append(H,[sdiv(G[J],F)]);
978: }
979: R = idealcap(R,H,VL,O);
980: }
981: }
982:
983: /* if ideal Q includes P then return 1 else return 0 */
984: /* Q must be a Groebner Basis w.r.t ordering ORD */
985: def idealinc(P,Q,VL,ORD)
986: {
987: dp_ord(ORD);
988: for (T=IN=[],I=length(Q)-1;I>=0;I--) {
989: T=cons(dp_ptod(Q[I],VL),T);
990: IN=cons(I,IN);
991: }
992: DQ = newvect(length(Q),T);
993: for (I = 0; I < length(P); I++) {
994: NF(E,IN,dp_ptod(P[I],VL),DQ,0);
995: if ( E )
996: return 0;
997: }
998: return 1;
999: }
1000:
1001: /* FL : list of polynomial set */
1002: def primedec_main(FL,VL)
1003: {
1004: if ( PRINTLOG ) {
1005: print("prime decomposition of the radical");
1006: }
1007: if ( PRINTLOG >= 2 )
1008: print("G-base factorization");
1009: PP = gr_fctr(FL,VL);
1010: if ( PRINTLOG >= 2 )
1011: print("irreducible variety decomposition");
1012: RP = welldec(PP,VL);
1013: SP = normposdec(RP,VL);
1014: for (L=[],I=length(SP)-1;I>=0;I--) {
1015: L=cons(idealnormal(SP[I]),L); /* head coef. > 0 */
1016: }
1017: SP = L;
1018: return SP;
1019: }
1020:
1021: def gr_fctr(FL,VL)
1022: {
1023: for (TP = [],I = 0; I<length(FL); I++ ) {
1024: F = FL[I];
1025: SF = idealsqfr(F);
1026: if ( !idealeq(F,SF) ) { GR(F,SF,VL,PRIMEORD); }
1027: DIRECTRY = []; COMMONIDEAL=[1];
1028: SP = gr_fctr_sub(F,VL);
1029: TP = append(TP,SP);
1030: }
1031: TP = idealsav(TP);
1032: TP = prime_irred(TP,VL);
1033: return TP;
1034: }
1035:
1036: def gr_fctr_sub(G,VL)
1037: {
1038: CONTINUE;
1039: if ( length(G) == 1 && type(G[0]) == 1 )
1040: return [G];
1041: RL = [];
1042: for (I = 0; I < length(G); I++) {
1043: FL = fctr(G[I]); L = length(FL); N = FL[1][1];
1044: if (L > 2 || N > 1) {
1045: TLL = [];
1046: for (J = 1; J < L; J++) {
1047: if ( CONTINUE ) {
1048: JCOPP=bload("jcopp");
1049: J = JCOPP[0]+1;
1050: COMMONIDEAL = JCOPP[1];
1051: RL = JCOPP[2];
1052: CONTINUE = 0;
1053: }
1054: if ( PRINTLOG >= 5 ) {
1055: for (D = length(DIRECTRY)-1; D >= 0; D--)
1056: print(DIRECTRY[D],0);
1057: print([L-1,J,length(RL)]);
1058: /*
1059: L-1:a number of factors of polynomial
1060: J:a factor executed now [0,...L-1]
1061: length(RL):a number of prime components
1062: */
1063: }
1064: W = cons(FL[J][0],G);
1065: GR(NG,W,VL,PRIMEORD);
1066: TNG = idealsqfr(NG);
1067: if ( NG != TNG ) { GR(NG,TNG,VL,PRIMEORD); }
1068: if ( idealinc(COMMONIDEAL,NG,VL,PRIMEORD) )
1069: continue;
1070: else {
1071: DIRECTRY=cons(L-J-1,DIRECTRY);
1072: DG = gr_fctr_sub(NG,VL);
1073: DIRECTRY=cdr(DIRECTRY);
1074: RL = append(RL,DG);
1075: if ( J <= L-2 && !idealinc(COMMONIDEAL,NG,VL,PRIMEORD)
1076: && COMMONCHECK ) {
1077: if ( PRINTLOG >= 5 ) {
1078: for (D = 0; D < length(DIRECTRY); D++) print(" ",2);
1079: print("intersection ideal");
1080: }
1081: COMMONIDEAL=idealcap(COMMONIDEAL,NG,VL,PRIMEORD);
1082: }
1083: }
1084: if ( BSAVE && !length(DIRECTRY) ) {
1085: bsave([J,COMMONIDEAL,RL],"jcopp");
1086: if ( PRINTLOG >= 2 )
1087: print("bsave j, intersection ideal and primes to icopp ");
1088: }
1089: }
1090: break;
1091: }
1092: }
1093: if (I == length(G)) {
1094: RL = append([G],RL);
1095: if ( PRINTLOG >= 5 ) {
1096: for (D = 0; D < length(DIRECTRY)-1; D++) print(" ",0);
1097: print("prime G-base ",0);
1098: if ( PRINTLOG >= 6 )
1099: print(G);
1100: else
1101: print("");
1102: }
1103: }
1104: return RL;
1105: }
1106:
1107: def prime_irred(TP,VL)
1108: {
1109: if ( PRINTLOG >= 4 )
1110: print("irredundancy check for prime ideal");
1111: N = length(TP);
1112: for (P = [], I = 0; I < N; I++) {
1113: if ( PRINTLOG >= 4 )
1114: print(".",2);
1115: for (J = 0; J < N; J++)
1116: if ( I != J && idealinc(TP[J],TP[I],VL,PRIMEORD) )
1117: break;
1118: if (J == N)
1119: P = append(P,[TP[I]]);
1120: }
1121: if ( PRINTLOG >= 4 )
1122: print("done.");
1123: return P;
1124: }
1125:
1126: /* if V1 \subset V2 then return 1 else return 0 */
1127: def varinc(V1,V2)
1128: {
1129: N1=length(V1); N2=length(V2);
1130: for (I=N1-1;I>=0;I--) {
1131: for (J=N2-1;J>=0;J--)
1132: if (V1[I]==V2[J])
1133: break;
1134: if (J==-1)
1135: return 0;
1136: }
1137: return 1;
1138: }
1139:
1140: def setunion(PS)
1141: {
1142: for (R=C=[]; PS != [] && car(PS); PS=cdr(PS)) {
1143: for (L = car(PS); L != []; L = cdr(L)) {
1144: A = car(L);
1145: for (G = C; G != []; G = cdr(G)) {
1146: if ( A == car(G) || -A == car(G) )
1147: break;
1148: }
1149: if ( G == [] ) {
1150: R = append(R,[A]);
1151: C = cons(A,C);
1152: }
1153: }
1154: }
1155: return R;
1156: }
1157:
1158: def idealeq(L,M)
1159: {
1160: if ((N1=length(L)) != (N2=length(M)))
1161: return 0;
1162: for (I = 0; I < length(L); I++) {
1163: for (J = 0; J < length(M); J++)
1164: if (L[I] == M[J] || L[I] == - M[J])
1165: break;
1166: if (J == length(M))
1167: return 0;
1168: }
1169: return 1;
1170: }
1171:
1172: def listminus(L,M)
1173: {
1174: for (R = []; L != []; L = cdr(L) ) {
1175: A = car(L);
1176: for (G = M; G != []; G = cdr(G)) {
1177: if ( A == car(G) )
1178: break;
1179: }
1180: if ( G == [] ) {
1181: R = append(R,[A]);
1182: M = cons(A,M);
1183: }
1184: }
1185: return R;
1186: }
1187:
1188: def idealsqfr(G)
1189: {
1190: for(Flag=0,LL=[],I=length(G)-1;I>=0;I--) {
1191: for(A=1,L=sqfr(G[I]),J=1;J<length(L);J++)
1192: A*=L[J][0];
1193: LL=cons(A,LL);
1194: }
1195: return LL;
1196: }
1197:
1198: def welldec(PRIME,VL)
1199: {
1200: PP = PRIME; NP = [];
1201: while ( !Flag ) {
1202: LP = [];
1203: Flag = 1;
1204: for (I=0;I<length(PP);I++) {
1205: U = welldec_sub(PP[I],VL);
1206: if (length(U) >= 2) {
1207: Flag = 0;
1208: if ( PRINTLOG >= 3 ) {
1209: print("now decomposing to irreducible varieties ",0);
1210: print(I,0); print(" ",0); print(length(PP));
1211: }
1212: DIRECTRY = []; COMMONIDEAL=[1];
1213: PL1 = gr_fctr([U[0]],VL);
1214: DIRECTRY = []; COMMONIDEAL=[1];
1215: PL2 = gr_fctr([U[1]],VL);
1216: LP = append(LP,append(PL1,PL2));
1217: }
1218: else
1219: NP = append(NP,U);
1220: }
1221: PP = LP;
1222: if ( PRINTLOG > 3 ) {
1223: print("prime length and non-prime length = ",0);
1224: print(length(NP),0); print(" ",0); print(length(PP));
1225: }
1226: }
1227: if ( length(PRIME) != length(NP) ) {
1228: NP = idealsav(NP);
1229: NP = prime_irred(NP,VL);
1230: }
1231: return NP;
1232: for (I = 0; I<length(PP);I++) {
1233: }
1234: }
1235:
1236: def welldec_sub(PP,VL)
1237: {
1238: VV = minalgdep(PP,VL,PRIMEORD);
1239: S = wellsep(PP,VV,VL);
1240: if ( S == 1 )
1241: return [PP];
1242: P1 = localize(PP,[S],VL,PRIMEORD);
1243: if ( idealeq(P1,PP) )
1244: return([PP]);
1245: GR(P2,cons(S,PP),VL,PRIMEORD);
1246: return [P1,P2];
1247: }
1248:
1249: def wellsep(PP,VV,VL)
1250: {
1251: TMPORD = 0;
1252: V0 = listminus(VL,VV);
1253: V1 = append(VV,V0);
1254: /* ORD = [[TMPORD,length(VV)],[0,length(V0)]]; */
1255: dp_nelim(length(VV)); ORD = 6;
1256: GR(PP,PP,V1,ORD);
1257: dp_ord(TMPORD);
1258: for (E=1,I=0;I<length(PP);I++)
1259: E = lcm(E,dp_hc(dp_ptod(PP[I],VV)));
1260: for (F=1,L=sqfr(E),K=1;K<length(L);K++)
1261: F *= L[K][0];
1262: return F;
1263: }
1264:
1265: /* prime decomposition by using primitive element method */
1266: def normposdec(NP,VL)
1267: {
1268: if ( PRINTLOG >= 3 )
1269: print("radical decomposition by using normal position.");
1270: for (R=MP=[],I=0;I<length(NP);I++) {
1271: V=minalgdep(NP[I],VL,PRIMEORD);
1272: L=raddec(NP[I],V,VL,1);
1273: if ( PRINTLOG >= 3 )
1274: if ( length(L) == 1 ) {
1275: print(".",2);
1276: } else {
1277: print(" ",2); print(length(L),2); print(" ",2);
1278: }
1279: /* if (length(L)==1) {
1280: MP=append(MP,[NP[I]]);
1281: continue;
1282: } */
1283: R=append(R,L);
1284: }
1285: if ( PRINTLOG >= 3 )
1286: print("done");
1287: if ( length(R) )
1288: MP = idealsav(append(R,MP));
1289: LP = prime_irred(MP,VL);
1290: return LP;
1291: }
1292:
1293: /* radical decomposition program using by Zianni, Trager, Zacharias */
1294: /* R : factorized G-base in K[X], if FLAG then A is well comp. */
1295: def raddec(R,V,X,FLAG)
1296: {
1297: GR(A,R,V,irem(PRIMEORD,3));
1298: ZQ=zraddec(A,V); /* zero dimensional */
1299: /* if ( FLAG && length(ZQ) == 1 )
1300: return [R]; */
1301: if ( length(V) != length(X) )
1302: CQ=radcont(ZQ,V,X); /* contraction */
1303: else {
1304: for (CQ=[],I=length(ZQ)-1;I>=0;I--) {
1305: GR(R,ZQ[I],X,PRIMEORD);
1306: CQ=cons(R,CQ);
1307: }
1308: }
1309: return CQ;
1310: }
1311:
1312: /* radical decomposition for zero dimensional ideal */
1313: /* F is G-base w.r.t irem(PRIMEORD,3) */
1314: def zraddec(F,X)
1315: {
1316: if (length(F) == 1)
1317: return [F];
1318: Z=vvv;
1319: C=normposit(F,X,Z);
1320: /* C=[minimal polynomial, adding polynomial] */
1321: T=cdr(fctr(C[0]));
1322: if ( length(T) == 1 && T[0][1] == 1 )
1323: return [F];
1324: for (Q=[],I=0;I<length(T);I++) {
1325: if ( !C[1] ) {
1326: GR(P,cons(T[I][0],F),X,irem(PRIMEORD,3));
1327: } else {
1328: P=idealgrsub(append([T[I][0],C[1]],F),X,Z);
1329: }
1330: Q=cons(P,Q);
1331: }
1332: return Q;
1333: }
1334:
1335: /* contraction from V to X */
1336: def radcont(Q,V,X)
1337: {
1338: dp_ord(irem(PRIMEORD,3));
1339: for (R=[],I=length(Q)-1;I>=0;I--) {
1340: G=Q[I];
1341: for (E=1,J=0;J<length(G);J++)
1342: E = lcm(E,dp_hc(dp_ptod(G[J],V)));
1343: for (F=1,L=sqfr(E),K=1;K<length(L);K++)
1344: F *= L[K][0];
1345: H=localize(G,[F],X,PRIMEORD);
1346: R=cons(H,R);
1347: }
1348: return R;
1349: }
1350:
1351: /* A : polynomials (G-Base) */
1352: /* return [T,Y], T : minimal polynomial of Z, Y : append polynomial */
1353: def normposit(A,X,Z)
1354: {
1355: D = idim(A,X,irem(PRIMEORD,3)); /* dimension of the ideal Id(A) */
1356: for (I = length(A)-1;I>=0;I--) {
1357: for (J = length(X)-1; J>= 0; J--) {
1358: T=deg(A[I],X[J]);
1359: if ( T == D ) /* A[I] is a primitive polynomial of A */
1360: return [A[I],0];
1361: }
1362: }
1363: N=length(X);
1364: for (C = [],I = 0; I < N-1; I++)
1365: C=newvect(N-1);
1366: V=append(X,[Z]);
1367: ZD = (vars(A)==vars(X)); /* A is zero dim. over Q[X] or not */
1368: while( 1 ) {
1369: C = nextweight(C,cdr(X));
1370: for (Y=Z-X[0],I=0;I<N-1;I++) {
1371: Y-=C[I]*X[I+1]; /* new polynomial */
1372: }
1373: if ( !ZD ) {
1374: if ( version() == 940420 ) NOREDUCE = 1;
1375: else dp_gr_flags(["NoRA",1]);
1376: GR(G,cons(Y,A),V,3);
1377: if ( version() == 940420 ) NOREDUCE = 0;
1378: else dp_gr_flags(["NoRA",0]);
1379: for (T=G[length(G)-1],I = length(G)-2;I >= 0; I--)
1380: if ( deg(G[I],Z) > deg(T,Z) )
1381: T = G[I]; /* minimal polynomial */
1382: } else {
1383: T = minipoly(A,X,irem(PRIMEORD,3),Z-Y,Z);
1384: }
1385: if (deg(T,Z)==D)
1386: return [T,Y];
1387: }
1388: }
1389:
1390: def nextweight(C,V) /* degrevlex */
1391: {
1392: dp_ord(2);
1393: N = length(V);
1394: for (D=I=0;I<size(C)[0];I++)
1395: D += C[I];
1396: if ( C[N-1] == D ) {
1397: for (L=[],I=0;I<N-1;I++)
1398: L=cons(0,L);
1399: L = cons(D+1,L);
1400: return newvect(N,L);
1401: }
1402: CD=dp_vtoe(C);
1403: for (F=I=0;I<N;I++)
1404: F+=V[I];
1405: for (DF=dp_ptod(F^D,V);dp_ht(DF)!=CD;DF=dp_rest(DF));
1406: MD = dp_ht(dp_rest(DF));
1407: CC = dp_etov(MD);
1408: return CC;
1409: }
1410:
1411: def printsep(L)
1412: {
1413: for (I=0;I<length(L);I++) {
1414: if ( nmono(L[I][0]) == 1 )
1415: print(L[I][0],0);
1416: else {
1417: print("(",0); print(L[I][0],0); print(")",0);
1418: }
1419: if (L[I][1] > 1) {
1420: print("^",0); print(L[I][1],0);
1421: }
1422: if (I<length(L)-1)
1423: print("*",0);
1424: }
1425: }
1426:
1427: def idealgrsub(A,X,Z)
1428: {
1429: if ( PRIMEORD == 0 ) {
1430: dp_nelim(1); ORD = 8;
1431: } else
1432: ORD = [[2,1],[PRIMERORD,length(X)]];
1433: GR(G,A,cons(Z,X),ORD);
1434: for (R=[],I=length(G)-1;I>=0;I--) {
1435: V=vars(G[I]);
1436: for (J=0;J<length(V);J++)
1437: if (V[J]==Z)
1438: break;
1439: if (J==length(V))
1440: R=cons(G[I],R);
1441: }
1442: return R;
1443: }
1444:
1445: /* dimension of ideal */
1446: def idim(F,V,ORD)
1447: {
1448: return length(modbase(F,V,ORD));
1449: }
1450:
1451: def modbase(F,V,ORD)
1452: {
1453: dp_ord(ORD);
1454: for(G=[],I=length(F)-1;I>=0;I--)
1455: G = cons(dp_ptod(F[I],V),G);
1456: R = dp_modbase(G,length(V));
1457: for(S=[],I=length(R)-1;I>=0;I--)
1458: S = cons(dp_dtop(R[I],V),S);
1459: return S;
1460: }
1461:
1462: def dp_modbase(G,N)
1463: {
1464: E = newvect(N);
1465: R = []; I = 1;
1466: for (L = [];G != []; G = cdr(G))
1467: L = cons(dp_ht(car(G)),L);
1468: while ( I <= N ) {
1469: R = cons(dp_vtoe(E),R);
1470: E[0]++; I = 1;
1471: while( s_redble(dp_vtoe(E),L) ) {
1472: E[I-1] = 0;
1473: if (I < N)
1474: E[I]++;
1475: I++;
1476: }
1477: }
1478: return R;
1479: }
1480:
1481: def s_redble(M,L)
1482: {
1483: for (; L != []; L = cdr(L))
1484: if ( dp_redble(M,car(L)) )
1485: return 1;
1486: return 0;
1487: }
1488:
1489: /* FL : list of ideal Id(P,s) */
1490: def primedec_special(FL,VL)
1491: {
1492: if ( PRINTLOG ) {
1493: print("prime decomposition of the radical");
1494: }
1495: if ( PRINTLOG >= 2 )
1496: print("G-base factorization");
1497: PP = gr_fctr(FL,VL);
1498: if ( PRINTLOG >= 2 )
1499: print("special decomposition");
1500: SP = yokodec(PP,VL);
1501: for (L=[],I=length(SP)-1;I>=0;I--) {
1502: L=cons(idealnormal(SP[I]),L); /* head coef. > 0 */
1503: }
1504: SP = L;
1505: return SP;
1506: }
1507:
1508: /* PL : list of ideal Id(P,s) */
1509: def yokodec(PL,VL)
1510: {
1511: MSISTIME = 0; T = time()[0];
1512: for (R = [],I = 0; I<length(PL);I++) {
1513: PP = PL[I];
1514: if ( PRINTLOG >= 3 ) print(".",2);
1515: V = minalgdep(PP,VL,PRIMEORD);
1516: E = raddec(PP,V,VL,0);
1517: if ( length(E) >= 2 || !idealeq(E[0],PP) ) { /* prime check */
1518: T0 = time()[0];
1519: L=all_msis(PP,VL,PRIMEORD);
1520: MSISTIME += time()[0]-T0;
1521: E = yokodec_main(PP,L,VL);
1522: }
1523: R = append(R,E);
1524: }
1525: R = idealsav(R);
1526: R = prime_irred(R,VL);
1527: if ( PRINTLOG >= 3 ) {
1528: print("special dec time = ",0); print(time()[0]-T0);
1529: /* print(", ",0); print(MSISTIME); */
1530: }
1531: return R;
1532: }
1533:
1534: def yokodec_main(PP,L,VL)
1535: {
1536: AL = misset(L,VL); H = [1]; R = [];
1537: for (P=PP,I=0; I<length(AL); I++) {
1538: V = AL[I];
1539: if ( I && !ismad(P,AL[I],VL,PRIMEORD) )
1540: continue;
1541: U = raddec(PP,V,VL,0);
1542: R = append(R,U);
1543: for (I=0;I<length(U);I++)
1544: H = idealcap(H,U[I],VL,PRIMEORD);
1545: if ( idealeq(H,P) )
1546: break;
1547: F = wellsep(P,V,VL);
1548: GR(P,cons(F,P),VL,PRIMEORD);
1549: }
1550: return R;
1551: }
1552:
1553: /* output M.A.D. sets. AL : M.S.I.S sets */
1554: def misset(AL,VL)
1555: {
1556: for (L=[],D=0,I=length(AL)-1;I>=0;I--) {
1557: E = length(AL[I]);
1558: C = listminus(VL,AL[I]);
1559: if ( E == D )
1560: L = cons(C,L);
1561: else if ( E > D ) {
1562: L = [C]; D = E;
1563: }
1564: }
1565: return L;
1566: }
1567:
1568: end$
1569: