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