![[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.8 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/primdec,v 1.7 2006/02/24 01:15:56 noro 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);
1.8 ! noro 460: #if 0
! 461: /* This may be a bug. GF is not a GB w.r.t the elimination order */
1.1 noro 462: if ( V0 != VL ) {
463: ORD = [[TMPORD1,length(V)],[TMPORD2,length(U)]];
464: GR(G,GF,V0,ORD);
465: } else
466: G = GF;
1.8 ! noro 467: #else
! 468: ORD = [[TMPORD1,length(V)],[TMPORD2,length(U)]];
! 469: GR(G,GF,V0,ORD);
! 470: #endif
1.1 noro 471: dp_ord(TMPORD1);
472: for (LL = [],HC = 1,I = 0; I < length(G); I++)
473: LL = append(LL,cdr(fctr(dp_hc(dp_ptod(G[I],V)))));
474: for (L=[],I=0;I<length(LL);I++)
475: L = cons(LL[I][0],L);
476: L = setunion([L]);
477: for (S=1,SL=[],I=0;I<length(L);I++) {
478: S *= L[I];
479: if ( SELECTFLAG )
480: SL = cons(L[I],SL);
481: }
482: if ( !SELECTFLAG )
483: SL= [S];
484: if ( PRINTLOG > 1 ) {
485: print("extractor = ",0);
486: print(S);
487: }
488: T0 = time()[0];
489: Q = localize(GF,SL,VL,PRIMAORD);
490: Q = idealnormal(Q);
491: DEC = [Q,Pr];
492: if ( PRINTLOG ) {
493: print("find a primary component! time = ",0);
494: print(time()[0]-T0);
495: if (PRINTLOG >= 3){
496: print(" associated prime of primary component = ",0);
497: print(DEC[1]);
498: }
499: if (PRINTLOG >= 4){print(" primary component = ",0); print(Q);}
500: }
501: if ( !ISOLATED && !NOEMBEDDED && SL != [1]
502: && length(V) != length(VL) /* nonzerodim */
503: && (REDUNDANT || !idealinc(Q,GF,VL,PRIMAORD)) ) {
504: REM = ["ext",[Q,Pr],GF,S];
505: if ( PRINTLOG ) {
506: print("find the remaining component of the extraction");
507: }
508: } else {
509: REM = [];
510: if ( PRINTLOG ) {
511: print("eliminate the remaining component of the extraction");
512: }
513: }
514: return [DEC,REM];
515: }
516:
517: /* sep. cond. for pseudo-primary components */
518: def sepcond_ps(P,SC,VL)
519: {
520: for (J = 0; J < length(SC); J++) {
521: if ( idealinc([SC[J]],P,VL,PRIMEORD) )
522: break; /* separating condition */
523: }
524: if ( J != length(SC) ) {
525: if ( PRINTLOG ) {
526: print("");
527: print("elim. the pseudo primary comp. by separating cond.");
528: }
529: return 0;
530: }
531: return 1;
532: }
533:
534: /* sep. cond. for rem. components. */
535: /* REM = ["ext",[Q,Pr],GF,S] or ["sep",PD,QD,GF,SEL], SC : sep.cond. */
536: def sepcond_re(REM,SC,VL)
537: {
538: for (S=1,I=0;I<length(SC);I++)
539: S *= SC[I];
540: if (REM[0] == "ext") {
541: F = cons(REM[3],REM[1][1]);
542: L = localize(F,[S],VL,PRIMAORD);
543: if ( L != [1] )
544: return 1;
545: else
546: return 0;
547: } else if (REM[0] == "sep") {
548: PL = REM[1]; SEL = REM[4];
549: for (I=0;I<length(PL);I++) {
550: F = append(SEL[I],PL[I]);
551: L = localize(F,[S],VL,PRIMAORD);
552: if ( L != [1] )
553: return 1;
554: }
555: return 0;
556: }
557: }
558:
559: def minalgdep(Pr,VL,ORD)
560: {
561: T0=newvect(4,time());
562: if (MAXALGIND)
563: R = minalgdep1(Pr,VL,ORD); /* M.I.S. */
564: else
565: R = minalgdep0(Pr,VL,ORD); /* M.S.I.S. */
566: T1=newvect(4,time())-T0; MISTIME += T1;
567: if ( PRINTLOG >= 6 ) {
568: if (MAXALGIND) print("M.I.S. time = ",0);
569: else print("M.S.I.S. time = ",0);
570: print(T1[0]);
571: }
572: return R;
573: }
574:
575: def minalgdep0(Pr,VL,ORD)
576: {
577: dp_ord(ORD); N = length(Pr);
578: for (HD = [],I = N-1; I >= 0; I--)
579: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
580: for (R = X = [], I = length(VL)-1; I >= 0; I--) {
581: V = VL[I];
582: TX = cons(V,X);
583: for (J = 0; J < N; J++) {
584: if ( varinc(HD[J],TX) )
585: break;
586: }
587: if ( J == N )
588: X = TX;
589: else
590: R = cons(V,R);
591: }
592: return R;
593: }
594:
595: def minalgdep1(Pr,VL,ORD)
596: {
597: R = all_msis(Pr,VL,ORD);
598: for (E=I=0,D=length(R[0]);I<length(R);I++) {
599: if ( length(R[I])>=D ) {
600: D=length(R[I]); E=I;
601: }
602: }
603: S = listminus(VL,R[E]);
604: return S;
605: }
606:
607: def all_msis(Pr,VL,ORD)
608: {
609: dp_ord(ORD); N = length(Pr);
610: for (HD = [],I = N-1; I >= 0; I--)
611: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
612: R = dimrec(HD,0,[],[],VL);
613: return R;
614: }
615:
616: def dimrec(S,K,U,M,X)
617: {
618: MM = M;
619: for (I=K;I<length(X);I++) {
620: TX = append(U,[X[I]]);
621: for (J = 0; J < length(S); J++) {
622: if ( varinc(S[J],TX) )
623: break;
624: }
625: if ( J == length(S) )
626: MM = dimrec(S,I+1,TX,MM,X);
627: }
628: for (J = 0; J < length(MM); J++) {
629: if ( varinc(U,MM[J]) )
630: break;
631: }
632: if ( J == length(MM) )
633: MM = append(MM,[U]);
634: return MM;
635: }
636:
637: /* if V is min.alg.dep. return 1 */
638: def ismad(Pr,V,VL,ORD)
639: {
640: dp_ord(ORD); N = length(Pr);
641: for (HD = [],I = N-1; I >= 0; I--)
642: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
643: W = listminus(VL,V); L = append(W,V);
644: for (R = TX = [],I = 0; I < length(L); I++) {
645: TX = cons(L[I],TX);
646: for (J = 0; J < N; J++) {
647: if ( varinc(HD[J],TX) )
648: break;
649: }
650: if ((I<length(W)&&J!=N) || (I>=length(W)&&J==N))
651: return 0;
652:
653: }
654: return 1;
655: }
656:
657: /* output:list of [flag, generator, primes, sep.cond., intersection] */
658: def mksepext(RES,H,VL)
659: {
660: for (R = [],I=length(RES)-1;I>=0;I--) {
661: W = RES[I];
662: if ( COMPLETE || EXTENDED )
663: HI = idealcap(H,W[6],VL,PRIMAORD);
664: else
665: HI = H;
666: if (W[1] == "sep") {
667: E = mksep(W[2],W[3],W[4],VL);
668: F = append(E[0],W[0]);
669: if ( ( COMPLETE || EXTENDED ) && !REDUNDANT )
670: HI = idealsetcap(cons(HI,W[3]),VL,PRIMAORD);
671: R = cons(["sep",F,E[1],W[5],HI],R);
672: } else if (W[1] == "ext") {
673: E = mkext(W[2][0],W[3],W[4],VL);
674: F = cons(E[0],W[0]);
675: P = cons(E[1],W[2][1]); /* prime component of rem. comp. */
676: R = cons(["ext",F,P,W[5],HI],R);
677: } else debug;
678: }
679: return R;
680: }
681:
682: /* make extractor with exponents ; F must be G-Base */
683: /* output : extractor, max.sqfr.of extractor */
684: def mkext(Q,F,S,VL)
685: {
686: if ( PRINTLOG > 1 ) {
687: print("make extractor = ",0);
688: }
689: if ( S == 1 || S == -1 ) {
690: print(1); return [1,1];
691: }
692: DD=dp_mkdist([Q],F,VL,PRIMAORD);
693: DQ=DD[0][0]; DF=DD[1];
694: for (L=fctr(S),SL=[],J=length(L)-1;J>0;J--)
695: SL=cons(dp_ptod(L[J][0],VL),SL);
696: K = dp_fexponent(dp_ptod(S,VL),DF,DQ);
697: E = dp_vectexp(SL,EX=dp_minexp(K,SL,DF,DQ));
698: if ( E != 1 ) E = dp_dtop(E,VL);
699: for (U=1, I=0;I<size(EX)[0];I++)
700: if ( EX[I] )
701: U *= L[I+1][0];
702: /* for (L=sqfr(E),U=1,J=length(L)-1;J>0;J--)
703: U *= L[J][0]; */
704: if ( PRINTLOG > 1 ) {
705: print(E);
706: }
707: return [E,U];
708: }
709:
710: /* make separators with exponents ; F must be G-Base */
711: /* output [separator, list of components of radical] */
712: def mksep(PP,QQ,F,VL)
713: {
714: if ( PRINTLOG > 1 ) {
715: print("make separators, ",0); print(length(PP));
716: }
717: DD=dp_mkdist(QQ,F,VL,PRIMAORD);
718: DQ=DD[0]; DF=DD[1];
719: E = dp_mksep(PP,DQ,DF,VL,PRIMAORD);
720: for (RA=[],I=length(E)-1;I>=0;I--) {
721: for (L=fctr(E[I]),J=length(L)-1;J>0;J--) {
722: ES=cons(L[J][0],PP[I]);
723: RA = cons(ES,RA); /* radical of remaining comp. */
724: }
725: }
726: return [E,RA];
727: }
728:
729: def dp_mksep(PP,DQ,DF,VL,ORD)
730: {
731: dp_ord(ORD);
732: for (E=[],I=0;I<length(PP);I++) {
733: T0=time()[0];
734: if ( CONTINUE ) {
735: E = bload("sepa");
736: I = length(E);
737: CONTINUE = 0;
738: }
739: S = selects(PP[I],PP,VL,0)[0];
740: for (L=fctr(S),SL=[],J=length(L)-1;J>0;J--)
741: SL=cons(dp_ptod(L[J][0],VL),SL);
742: K = dp_fexponent(dp_ptod(S,VL),DF,DQ[I]);
743: M = dp_vectexp(SL,dp_minexp(K,SL,DF,DQ[I]));
744: if (M != 1) M = dp_dtop(M,VL);
745: E = append(E,[M]);
746: if ( PRINTLOG > 1 ) {
747: print("separator of ",0); print(I,0);
748: print(" in ",0); print(length(PP),0);
749: print(", time = ",0); print(time()[0]-T0,0); print(", ",0);
750: printsep(cdr(fctr(M))); print("");
751: }
752: if (BSAVE) {
753: bsave(E,"sepa");
754: print("bsave to sepa");
755: }
756: }
757: return E;
758: }
759:
760: extern AFO$
761:
762: /* F,Id,Q is given by distributed polynomials. Id is a G-Base.*/
763: /* get the exponent K s.t. (Id : F^K) = Q, Id is a G-Basis */
764: def dp_fexponent(F,Id,Q)
765: {
766: if (!AFO)
767: for (IN=[],I=size(Id)[0]-1;I>=0;I--) IN = cons(I,IN);
768: else
769: for (IN=[],I=0;I<size(Id)[0];I++) IN = cons(I,IN);
770: NF(G,IN,F,Id,0);
771: if ( F != G ) { F = G; }
772: N=size(Q)[0]; U = newvect(N);
773: for (I=0;I<N;I++)
774: U[I]=Q[I];
775: for (K = 1;; K++) { /* U[I] = Q[I]*F^K mod Id */
776: for (FLAG = I = 0; I < N; I++) {
777: NF(U[I],IN,U[I]*F,Id,1);
778: if ( U[I] ) FLAG = 1;
779: }
780: if ( !FLAG )
781: return K;
782: }
783: }
784:
785: /* F must be a G-Base. SL is a set of factors of S.*/
786: def dp_minexp(K,SL,F,Q)
787: {
788: if (!AFO)
789: for (IN=[],I=size(F)[0]-1;I>=0;I--) IN = cons(I,IN);
790: else
791: for (IN=[],I=0;I<size(F)[0];I++) IN = cons(I,IN);
792: N = length(SL); M = size(Q)[0];
793: EX = newvect(N); U = newvect(M); T = newvect(M);
794: for (I=0;I<N;I++)
795: EX[I]=K;
796: for (I=0;I<M;I++) {
797: NF(Q[I],IN,Q[I],F,1);
798: }
799: for (I=0;I<N;I++) {
800: EX[I] = 0;
801: for (FF=1,II=0;II<N;II++)
802: for (J = 0; J < EX[II]; J++) {
803: NF(FF,IN,FF*SL[II],F,1);
804: }
805: for (J = 0; J < M; J++)
806: U[J] = Q[J]*FF;
807: for (; EX[I] < K; EX[I]++) {
808: FLAG = 0;
809: for (J = 0; J < M; J++) {
810: NF(U[J],IN,U[J],F,1);
811: if ( U[J] ) FLAG = 1;
812: }
813: if ( !FLAG )
814: break;
815: if ( EX[I] < K-1 )
816: for (J = 0; J < M; J++) U[J] = U[J]*SL[I];
817: }
818: }
819: return EX;
820: }
821:
822: def dp_vectexp(SS,EE)
823: {
824: N = length(SS);
825: for (A=1,I=0;I<N;I++)
826: for (J = 0; J < EE[I]; J++)
827: A *= SS[I];
828: return A;
829: }
830:
831: def dp_mkdist(QQ,F,VL,ORD)
832: {
833: dp_ord(ORD);
834: for (DQ=[],I=length(QQ)-1;I>=0;I--) {
835: for (T=[],J=length(QQ[I])-1;J>=0;J--)
836: T = cons(dp_ptod(QQ[I][J],VL),T);
837: DQ=cons(newvect(length(T),T),DQ);
838: }
839: for (T=[],J=length(F)-1;J>=0;J--)
840: T = cons(dp_ptod(F[J],VL),T);
841: DF = newvect(length(T),T);
842: return [DQ,DF]$
843: }
844:
845: /* if FLAG == 0 return singltonset */
846: def selects(P,PD,VL,FLAG)
847: {
848: PP=dp_mkdist([],P,VL,PRIMEORD)[1];
849: for (M=[],I=length(P)-1;I>=0;I--)
850: M = cons(I,M);
851: for (SL = [],I = length(PD)-1; I >= 0; I--) {
852: B = PD[I];
853: if (B == P) continue;
854: for (L = [],J = 0; J < length(B); J++) {
855: A = B[J];
856: NF(E,M,dp_ptod(A,VL),PP,0);
857: if ( E ) {
858: L = append(L,[A]); /* A \notin Id(P) */
859: }
860: }
861: SL = cons(L,SL); /* candidate of separators */
862: }
863: for (S=[],I=length(SL)-1;I>=0;I--) {/* choose an element from L */
864: if ( FLAG >= 2 ) {
865: S = append(SL[I],S);
866: continue;
867: }
868: C = minselect(SL[I],VL,PRIMAORD);
869: S = cons(C,S);
870: }
871: S = setunion([S]);
872: if ( FLAG == 3 || length(S) == 1 )
873: return S;
874: if ( FLAG == 2 ) {
875: for (R=1,I=0;I<length(S);I++)
876: R *= S[I];
877: return [R];
878: }
879: /* FLAG == 0 || FLAG == 1 */
880: for (T=[]; S!=[];S = cdr(S)) { /* FLAG <= 1 and length(S) != 1 */
881: A = car(S);
882: U = append(T,cdr(S));
883: for (R=1,I=0;I<length(U);I++)
884: R *= U[I];
885: for (I=0;I<length(PD);I++) {
886: if ( PD[I] != P && !idealinc([R],PD[I],VL,PRIMAORD) )
887: break;
888: }
889: if ( I != length(PD) )
890: T = append(T,[A]);
891: }
892: if ( FLAG )
893: return T;
894: for (R=1,I=0;I<length(T);I++)
895: R *= T[I];
896: return [R]; /* if FLAG == 0 */
897: }
898:
899: def minselect(L,VL,ORD)
900: {
901: dp_ord(ORD);
902: F = dp_ptod(L[0],VL); MAG=dp_mag(F); DEG=dp_td(F);
903: for (J = 0, I = 1; I < length(L); I++) {
904: A = dp_ptod(L[I],VL); M=dp_mag(A); D=dp_td(A);
905: if ( dp_ht(A) > dp_ht(F) )
906: continue;
907: else if ( dp_ht(A) == dp_ht(F) && (D > DEG || (D == DEG && M > MAG)))
908: continue;
909: F = A; J = I; MAG = M; DEG=D;
910: }
911: return L[J];
912: }
913:
914: /* localization Id(F)_MI \cap Q[VL] */
915: /* output is the G-base w.r.t ordering O */
916: def localize(F,MI,VL,O)
917: {
918: if ( MI == [1] || MI == [-1] )
919: return F;
920: for (SVL = [],R = [],I = 0; I < length(MI); I++) {
921: V = strtov("zz"+rtostr(I));
922: R = append(R,[V*MI[I]-1]);
923: SVL = append(SVL,[V]);
924: }
925: if ( O == 0 ) {
926: dp_nelim(length(SVL)); ORD = 6;
927: } else if ( O == 1 || O == 2 ) {
928: ORD = [[0,length(SVL)],[O,length(VL)]];
929: } else if ( O == 3 || O == 4 || O == 5 )
930: ORD = [[0,length(SVL)],[O-3,length(VL)-1],[2,1]];
931: else
932: error("invarid ordering");
933: GR(G,append(F,R),append(SVL,VL),ORD);
934: S = varminus(G,SVL);
935: return S;
936: }
937:
938: def varminus(G,VL)
939: {
940: for (S = [],I = 0; I < length(G); I++) {
941: V = vars(G[I]);
942: if (listminus(V,VL) == V)
943: S = append(S,[G[I]]);
944: }
945: return S;
946: }
947:
948: def idealnormal(L)
949: {
950: R=[];
951: for (I=length(L)-1;I>=0;I--) {
952: A = ptozp(L[I]);
953: V = vars(A);
954: for (B = A,J = 0;J < length(V);J++)
955: B = coef(B,deg(B,V[J]),V[J]);
956: R = cons((B>0?A:-A),R);
957: }
958: return R;
959: }
960:
961: def idealsav(L) /* sub procedure of welldec and normposdec */
962: {
963: if ( PRINTLOG >= 4 )
964: print("multiple check ",2);
965: for (R = [],I = 0,N=length(L); I < N; I++) {
966: if ( PRINTLOG >= 4 && !irem(I,10) )
967: print(".",2);
968: if ( R == [] )
969: R = append(R,[L[I]]);
970: else {
971: for (J = 0; J < length(R); J++)
972: if ( idealeq(L[I],R[J]) )
973: break;
974: if ( J == length(R) )
975: R = append(R,[L[I]]);
976: }
977: }
978: if ( PRINTLOG >= 4 )
979: print("done.");
980: return R;
981: }
982:
983: /* intersection of ideals in PL, PL : list of ideals */
984: /* VL : variable list, O : output the G-base w.r.t order O */
985: def idealsetcap(PL,VL,O)
986: {
987: for (U = [], I = 0; I < length(PL); I++)
988: if ( PL[I] != [] )
989: U = append(U,[PL[I]]);
990: if ( U == [] )
991: return [];
992: if ( PRINTLOG >= 4 ) {print("intersection of ideal ",0); print(length(U),0);}
993: for (L = U[0],I=1;I<length(U);I++) {
994: if ( PRINTLOG >=4 ) print(".",2);
995: L = idealcap(L,U[I],VL,O);
996: }
997: if ( PRINTLOG >=4 ) print("");
998: return L;
999: }
1000:
1001: /* return intersection set of ideals P and Q */
1002: def idealcap(P,Q,VL,O)
1003: {
1004: if ( P == [] )
1005: if ( VL )
1006: return Q;
1007: else
1008: return [];
1009: if ( Q == [] )
1010: return P;
1011: T = tmp;
1012: for (PP = [],I = 0; I < length(P); I++)
1013: PP = append(PP,[P[I]*T]);
1014: for (QQ = [],I = 0; I < length(Q); I++)
1015: QQ = append(QQ,[Q[I]*(T-1)]);
1016: if ( O == 0 ) {
1017: dp_nelim(1); ORD = 6;
1018: } else if ( O == 1 || O == 2 )
1019: ORD = [[2,1],[O,length(VL)]];
1020: else if ( O == 3 || O == 4 || O == 5 )
1021: ORD = [[2,1],[O-3,length(VL)-1],[2,1]];
1022: else
1023: error("invarid ordering");
1024: GR(G,append(PP,QQ),append([T],VL),ORD);
1025: S = varminus(G,[T]);
1026: return S;
1027: }
1028:
1029: /* return ideal P : Q */
1030: def idealquo(P,Q,VL,O)
1031: {
1032: for (R = [], I = 0; I < length(Q); I++) {
1033: F = car(Q);
1034: G = idealcap(P,[F],VL,O);
1035: for (H = [],J = 0; J < length(G); J++) {
1036: H = append(H,[sdiv(G[J],F)]);
1037: }
1038: R = idealcap(R,H,VL,O);
1039: }
1040: }
1041:
1042: /* if ideal Q includes P then return 1 else return 0 */
1043: /* Q must be a Groebner Basis w.r.t ordering ORD */
1044: def idealinc(P,Q,VL,ORD)
1045: {
1046: dp_ord(ORD);
1047: for (T=IN=[],I=length(Q)-1;I>=0;I--) {
1048: T=cons(dp_ptod(Q[I],VL),T);
1049: IN=cons(I,IN);
1050: }
1051: DQ = newvect(length(Q),T);
1052: for (I = 0; I < length(P); I++) {
1053: NF(E,IN,dp_ptod(P[I],VL),DQ,0);
1054: if ( E )
1055: return 0;
1056: }
1057: return 1;
1058: }
1059:
1060: /* FL : list of polynomial set */
1061: def primedec_main(FL,VL)
1062: {
1063: if ( PRINTLOG ) {
1064: print("prime decomposition of the radical");
1065: }
1066: if ( PRINTLOG >= 2 )
1067: print("G-base factorization");
1068: PP = gr_fctr(FL,VL);
1069: if ( PRINTLOG >= 2 )
1070: print("irreducible variety decomposition");
1071: RP = welldec(PP,VL);
1072: SP = normposdec(RP,VL);
1073: for (L=[],I=length(SP)-1;I>=0;I--) {
1074: L=cons(idealnormal(SP[I]),L); /* head coef. > 0 */
1075: }
1076: SP = L;
1077: return SP;
1078: }
1079:
1080: def gr_fctr(FL,VL)
1081: {
1082: for (TP = [],I = 0; I<length(FL); I++ ) {
1083: F = FL[I];
1084: SF = idealsqfr(F);
1085: if ( !idealeq(F,SF) ) { GR(F,SF,VL,PRIMEORD); }
1086: DIRECTRY = []; COMMONIDEAL=[1];
1087: SP = gr_fctr_sub(F,VL);
1088: TP = append(TP,SP);
1089: }
1090: TP = idealsav(TP);
1091: TP = prime_irred(TP,VL);
1092: return TP;
1093: }
1094:
1095: def gr_fctr_sub(G,VL)
1096: {
1097: CONTINUE;
1098: if ( length(G) == 1 && type(G[0]) == 1 )
1099: return [G];
1100: RL = [];
1101: for (I = 0; I < length(G); I++) {
1102: FL = fctr(G[I]); L = length(FL); N = FL[1][1];
1103: if (L > 2 || N > 1) {
1104: TLL = [];
1105: for (J = 1; J < L; J++) {
1106: if ( CONTINUE ) {
1107: JCOPP=bload("jcopp");
1108: J = JCOPP[0]+1;
1109: COMMONIDEAL = JCOPP[1];
1110: RL = JCOPP[2];
1111: CONTINUE = 0;
1112: }
1113: if ( PRINTLOG >= 5 ) {
1114: for (D = length(DIRECTRY)-1; D >= 0; D--)
1115: print(DIRECTRY[D],0);
1116: print([L-1,J,length(RL)]);
1117: /*
1118: L-1:a number of factors of polynomial
1119: J:a factor executed now [0,...L-1]
1120: length(RL):a number of prime components
1121: */
1122: }
1123: W = cons(FL[J][0],G);
1124: GR(NG,W,VL,PRIMEORD);
1125: TNG = idealsqfr(NG);
1126: if ( NG != TNG ) { GR(NG,TNG,VL,PRIMEORD); }
1127: if ( idealinc(COMMONIDEAL,NG,VL,PRIMEORD) )
1128: continue;
1129: else {
1130: DIRECTRY=cons(L-J-1,DIRECTRY);
1131: DG = gr_fctr_sub(NG,VL);
1132: DIRECTRY=cdr(DIRECTRY);
1133: RL = append(RL,DG);
1134: if ( J <= L-2 && !idealinc(COMMONIDEAL,NG,VL,PRIMEORD)
1135: && COMMONCHECK ) {
1136: if ( PRINTLOG >= 5 ) {
1137: for (D = 0; D < length(DIRECTRY); D++) print(" ",2);
1138: print("intersection ideal");
1139: }
1140: COMMONIDEAL=idealcap(COMMONIDEAL,NG,VL,PRIMEORD);
1141: }
1142: }
1143: if ( BSAVE && !length(DIRECTRY) ) {
1144: bsave([J,COMMONIDEAL,RL],"jcopp");
1145: if ( PRINTLOG >= 2 )
1146: print("bsave j, intersection ideal and primes to icopp ");
1147: }
1148: }
1149: break;
1150: }
1151: }
1152: if (I == length(G)) {
1153: RL = append([G],RL);
1154: if ( PRINTLOG >= 5 ) {
1155: for (D = 0; D < length(DIRECTRY)-1; D++) print(" ",0);
1156: print("prime G-base ",0);
1157: if ( PRINTLOG >= 6 )
1158: print(G);
1159: else
1160: print("");
1161: }
1162: }
1163: return RL;
1164: }
1165:
1166: def prime_irred(TP,VL)
1167: {
1168: if ( PRINTLOG >= 4 )
1169: print("irredundancy check for prime ideal");
1170: N = length(TP);
1171: for (P = [], I = 0; I < N; I++) {
1172: if ( PRINTLOG >= 4 )
1173: print(".",2);
1174: for (J = 0; J < N; J++)
1175: if ( I != J && idealinc(TP[J],TP[I],VL,PRIMEORD) )
1176: break;
1177: if (J == N)
1178: P = append(P,[TP[I]]);
1179: }
1180: if ( PRINTLOG >= 4 )
1181: print("done.");
1182: return P;
1183: }
1184:
1185: /* if V1 \subset V2 then return 1 else return 0 */
1186: def varinc(V1,V2)
1187: {
1188: N1=length(V1); N2=length(V2);
1189: for (I=N1-1;I>=0;I--) {
1190: for (J=N2-1;J>=0;J--)
1191: if (V1[I]==V2[J])
1192: break;
1193: if (J==-1)
1194: return 0;
1195: }
1196: return 1;
1197: }
1198:
1199: def setunion(PS)
1200: {
1201: for (R=C=[]; PS != [] && car(PS); PS=cdr(PS)) {
1202: for (L = car(PS); L != []; L = cdr(L)) {
1203: A = car(L);
1204: for (G = C; G != []; G = cdr(G)) {
1205: if ( A == car(G) || -A == car(G) )
1206: break;
1207: }
1208: if ( G == [] ) {
1209: R = append(R,[A]);
1210: C = cons(A,C);
1211: }
1212: }
1213: }
1214: return R;
1215: }
1216:
1217: def idealeq(L,M)
1218: {
1219: if ((N1=length(L)) != (N2=length(M)))
1220: return 0;
1221: for (I = 0; I < length(L); I++) {
1222: for (J = 0; J < length(M); J++)
1223: if (L[I] == M[J] || L[I] == - M[J])
1224: break;
1225: if (J == length(M))
1226: return 0;
1227: }
1228: return 1;
1229: }
1230:
1231: def listminus(L,M)
1232: {
1233: for (R = []; L != []; L = cdr(L) ) {
1234: A = car(L);
1235: for (G = M; G != []; G = cdr(G)) {
1236: if ( A == car(G) )
1237: break;
1238: }
1239: if ( G == [] ) {
1240: R = append(R,[A]);
1241: M = cons(A,M);
1242: }
1243: }
1244: return R;
1245: }
1246:
1247: def idealsqfr(G)
1248: {
1249: for(Flag=0,LL=[],I=length(G)-1;I>=0;I--) {
1250: for(A=1,L=sqfr(G[I]),J=1;J<length(L);J++)
1251: A*=L[J][0];
1252: LL=cons(A,LL);
1253: }
1254: return LL;
1255: }
1256:
1257: def welldec(PRIME,VL)
1258: {
1259: PP = PRIME; NP = [];
1260: while ( !Flag ) {
1261: LP = [];
1262: Flag = 1;
1263: for (I=0;I<length(PP);I++) {
1264: U = welldec_sub(PP[I],VL);
1265: if (length(U) >= 2) {
1266: Flag = 0;
1267: if ( PRINTLOG >= 3 ) {
1268: print("now decomposing to irreducible varieties ",0);
1269: print(I,0); print(" ",0); print(length(PP));
1270: }
1271: DIRECTRY = []; COMMONIDEAL=[1];
1272: PL1 = gr_fctr([U[0]],VL);
1273: DIRECTRY = []; COMMONIDEAL=[1];
1274: PL2 = gr_fctr([U[1]],VL);
1275: LP = append(LP,append(PL1,PL2));
1276: }
1277: else
1278: NP = append(NP,U);
1279: }
1280: PP = LP;
1281: if ( PRINTLOG > 3 ) {
1282: print("prime length and non-prime length = ",0);
1283: print(length(NP),0); print(" ",0); print(length(PP));
1284: }
1285: }
1286: if ( length(PRIME) != length(NP) ) {
1287: NP = idealsav(NP);
1288: NP = prime_irred(NP,VL);
1289: }
1290: return NP;
1291: for (I = 0; I<length(PP);I++) {
1292: }
1293: }
1294:
1295: def welldec_sub(PP,VL)
1296: {
1297: VV = minalgdep(PP,VL,PRIMEORD);
1298: S = wellsep(PP,VV,VL);
1299: if ( S == 1 )
1300: return [PP];
1301: P1 = localize(PP,[S],VL,PRIMEORD);
1302: if ( idealeq(P1,PP) )
1303: return([PP]);
1304: GR(P2,cons(S,PP),VL,PRIMEORD);
1305: return [P1,P2];
1306: }
1307:
1308: def wellsep(PP,VV,VL)
1309: {
1310: TMPORD = 0;
1311: V0 = listminus(VL,VV);
1312: V1 = append(VV,V0);
1313: /* ORD = [[TMPORD,length(VV)],[0,length(V0)]]; */
1314: dp_nelim(length(VV)); ORD = 6;
1315: GR(PP,PP,V1,ORD);
1316: dp_ord(TMPORD);
1317: for (E=1,I=0;I<length(PP);I++)
1318: E = lcm(E,dp_hc(dp_ptod(PP[I],VV)));
1319: for (F=1,L=sqfr(E),K=1;K<length(L);K++)
1320: F *= L[K][0];
1321: return F;
1322: }
1323:
1324: /* prime decomposition by using primitive element method */
1325: def normposdec(NP,VL)
1326: {
1327: if ( PRINTLOG >= 3 )
1328: print("radical decomposition by using normal position.");
1329: for (R=MP=[],I=0;I<length(NP);I++) {
1330: V=minalgdep(NP[I],VL,PRIMEORD);
1331: L=raddec(NP[I],V,VL,1);
1332: if ( PRINTLOG >= 3 )
1333: if ( length(L) == 1 ) {
1334: print(".",2);
1335: } else {
1336: print(" ",2); print(length(L),2); print(" ",2);
1337: }
1338: /* if (length(L)==1) {
1339: MP=append(MP,[NP[I]]);
1340: continue;
1341: } */
1342: R=append(R,L);
1343: }
1344: if ( PRINTLOG >= 3 )
1345: print("done");
1346: if ( length(R) )
1347: MP = idealsav(append(R,MP));
1348: LP = prime_irred(MP,VL);
1349: return LP;
1350: }
1351:
1352: /* radical decomposition program using by Zianni, Trager, Zacharias */
1353: /* R : factorized G-base in K[X], if FLAG then A is well comp. */
1354: def raddec(R,V,X,FLAG)
1355: {
1356: GR(A,R,V,irem(PRIMEORD,3));
1357: ZQ=zraddec(A,V); /* zero dimensional */
1358: /* if ( FLAG && length(ZQ) == 1 )
1359: return [R]; */
1360: if ( length(V) != length(X) )
1361: CQ=radcont(ZQ,V,X); /* contraction */
1362: else {
1363: for (CQ=[],I=length(ZQ)-1;I>=0;I--) {
1364: GR(R,ZQ[I],X,PRIMEORD);
1365: CQ=cons(R,CQ);
1366: }
1367: }
1368: return CQ;
1369: }
1370:
1371: /* radical decomposition for zero dimensional ideal */
1372: /* F is G-base w.r.t irem(PRIMEORD,3) */
1373: def zraddec(F,X)
1374: {
1375: if (length(F) == 1)
1376: return [F];
1377: Z=vvv;
1378: C=normposit(F,X,Z);
1379: /* C=[minimal polynomial, adding polynomial] */
1380: T=cdr(fctr(C[0]));
1381: if ( length(T) == 1 && T[0][1] == 1 )
1382: return [F];
1383: for (Q=[],I=0;I<length(T);I++) {
1384: if ( !C[1] ) {
1385: GR(P,cons(T[I][0],F),X,irem(PRIMEORD,3));
1386: } else {
1387: P=idealgrsub(append([T[I][0],C[1]],F),X,Z);
1388: }
1389: Q=cons(P,Q);
1390: }
1391: return Q;
1392: }
1393:
1394: /* contraction from V to X */
1395: def radcont(Q,V,X)
1396: {
1397: for (R=[],I=length(Q)-1;I>=0;I--) {
1.7 noro 1398: dp_ord(irem(PRIMEORD,3));
1.1 noro 1399: G=Q[I];
1400: for (E=1,J=0;J<length(G);J++)
1401: E = lcm(E,dp_hc(dp_ptod(G[J],V)));
1402: for (F=1,L=sqfr(E),K=1;K<length(L);K++)
1403: F *= L[K][0];
1404: H=localize(G,[F],X,PRIMEORD);
1405: R=cons(H,R);
1406: }
1407: return R;
1408: }
1409:
1410: /* A : polynomials (G-Base) */
1411: /* return [T,Y], T : minimal polynomial of Z, Y : append polynomial */
1412: def normposit(A,X,Z)
1413: {
1414: D = idim(A,X,irem(PRIMEORD,3)); /* dimension of the ideal Id(A) */
1415: for (I = length(A)-1;I>=0;I--) {
1416: for (J = length(X)-1; J>= 0; J--) {
1417: T=deg(A[I],X[J]);
1418: if ( T == D ) /* A[I] is a primitive polynomial of A */
1419: return [A[I],0];
1420: }
1421: }
1422: N=length(X);
1423: for (C = [],I = 0; I < N-1; I++)
1424: C=newvect(N-1);
1425: V=append(X,[Z]);
1426: ZD = (vars(A)==vars(X)); /* A is zero dim. over Q[X] or not */
1427: while( 1 ) {
1428: C = nextweight(C,cdr(X));
1429: for (Y=Z-X[0],I=0;I<N-1;I++) {
1430: Y-=C[I]*X[I+1]; /* new polynomial */
1431: }
1432: if ( !ZD ) {
1433: if ( version() == 940420 ) NOREDUCE = 1;
1434: else dp_gr_flags(["NoRA",1]);
1435: GR(G,cons(Y,A),V,3);
1436: if ( version() == 940420 ) NOREDUCE = 0;
1437: else dp_gr_flags(["NoRA",0]);
1438: for (T=G[length(G)-1],I = length(G)-2;I >= 0; I--)
1439: if ( deg(G[I],Z) > deg(T,Z) )
1440: T = G[I]; /* minimal polynomial */
1441: } else {
1442: T = minipoly(A,X,irem(PRIMEORD,3),Z-Y,Z);
1443: }
1444: if (deg(T,Z)==D)
1445: return [T,Y];
1446: }
1447: }
1448:
1449: def nextweight(C,V) /* degrevlex */
1450: {
1451: dp_ord(2);
1452: N = length(V);
1453: for (D=I=0;I<size(C)[0];I++)
1454: D += C[I];
1455: if ( C[N-1] == D ) {
1456: for (L=[],I=0;I<N-1;I++)
1457: L=cons(0,L);
1458: L = cons(D+1,L);
1459: return newvect(N,L);
1460: }
1461: CD=dp_vtoe(C);
1462: for (F=I=0;I<N;I++)
1463: F+=V[I];
1464: for (DF=dp_ptod(F^D,V);dp_ht(DF)!=CD;DF=dp_rest(DF));
1465: MD = dp_ht(dp_rest(DF));
1466: CC = dp_etov(MD);
1467: return CC;
1468: }
1469:
1470: def printsep(L)
1471: {
1472: for (I=0;I<length(L);I++) {
1473: if ( nmono(L[I][0]) == 1 )
1474: print(L[I][0],0);
1475: else {
1476: print("(",0); print(L[I][0],0); print(")",0);
1477: }
1478: if (L[I][1] > 1) {
1479: print("^",0); print(L[I][1],0);
1480: }
1481: if (I<length(L)-1)
1482: print("*",0);
1483: }
1484: }
1485:
1486: def idealgrsub(A,X,Z)
1487: {
1488: if ( PRIMEORD == 0 ) {
1489: dp_nelim(1); ORD = 8;
1490: } else
1491: ORD = [[2,1],[PRIMERORD,length(X)]];
1492: GR(G,A,cons(Z,X),ORD);
1493: for (R=[],I=length(G)-1;I>=0;I--) {
1494: V=vars(G[I]);
1495: for (J=0;J<length(V);J++)
1496: if (V[J]==Z)
1497: break;
1498: if (J==length(V))
1499: R=cons(G[I],R);
1500: }
1501: return R;
1502: }
1503:
1504: /* dimension of ideal */
1505: def idim(F,V,ORD)
1506: {
1507: return length(modbase(F,V,ORD));
1508: }
1509:
1510: def modbase(F,V,ORD)
1511: {
1512: dp_ord(ORD);
1513: for(G=[],I=length(F)-1;I>=0;I--)
1514: G = cons(dp_ptod(F[I],V),G);
1515: R = dp_modbase(G,length(V));
1516: for(S=[],I=length(R)-1;I>=0;I--)
1517: S = cons(dp_dtop(R[I],V),S);
1518: return S;
1519: }
1520:
1521: def dp_modbase(G,N)
1522: {
1523: E = newvect(N);
1524: R = []; I = 1;
1525: for (L = [];G != []; G = cdr(G))
1526: L = cons(dp_ht(car(G)),L);
1527: while ( I <= N ) {
1528: R = cons(dp_vtoe(E),R);
1529: E[0]++; I = 1;
1530: while( s_redble(dp_vtoe(E),L) ) {
1531: E[I-1] = 0;
1532: if (I < N)
1533: E[I]++;
1534: I++;
1535: }
1536: }
1537: return R;
1538: }
1539:
1540: def s_redble(M,L)
1541: {
1542: for (; L != []; L = cdr(L))
1543: if ( dp_redble(M,car(L)) )
1544: return 1;
1545: return 0;
1546: }
1547:
1548: /* FL : list of ideal Id(P,s) */
1549: def primedec_special(FL,VL)
1550: {
1551: if ( PRINTLOG ) {
1552: print("prime decomposition of the radical");
1553: }
1554: if ( PRINTLOG >= 2 )
1555: print("G-base factorization");
1556: PP = gr_fctr(FL,VL);
1557: if ( PRINTLOG >= 2 )
1558: print("special decomposition");
1559: SP = yokodec(PP,VL);
1560: for (L=[],I=length(SP)-1;I>=0;I--) {
1561: L=cons(idealnormal(SP[I]),L); /* head coef. > 0 */
1562: }
1563: SP = L;
1564: return SP;
1565: }
1566:
1567: /* PL : list of ideal Id(P,s) */
1568: def yokodec(PL,VL)
1569: {
1570: MSISTIME = 0; T = time()[0];
1571: for (R = [],I = 0; I<length(PL);I++) {
1572: PP = PL[I];
1573: if ( PRINTLOG >= 3 ) print(".",2);
1574: V = minalgdep(PP,VL,PRIMEORD);
1575: E = raddec(PP,V,VL,0);
1576: if ( length(E) >= 2 || !idealeq(E[0],PP) ) { /* prime check */
1577: T0 = time()[0];
1578: L=all_msis(PP,VL,PRIMEORD);
1579: MSISTIME += time()[0]-T0;
1580: E = yokodec_main(PP,L,VL);
1581: }
1582: R = append(R,E);
1583: }
1584: R = idealsav(R);
1585: R = prime_irred(R,VL);
1586: if ( PRINTLOG >= 3 ) {
1587: print("special dec time = ",0); print(time()[0]-T0);
1588: /* print(", ",0); print(MSISTIME); */
1589: }
1590: return R;
1591: }
1592:
1593: def yokodec_main(PP,L,VL)
1594: {
1595: AL = misset(L,VL); H = [1]; R = [];
1596: for (P=PP,I=0; I<length(AL); I++) {
1597: V = AL[I];
1598: if ( I && !ismad(P,AL[I],VL,PRIMEORD) )
1599: continue;
1600: U = raddec(PP,V,VL,0);
1601: R = append(R,U);
1602: for (I=0;I<length(U);I++)
1603: H = idealcap(H,U[I],VL,PRIMEORD);
1604: if ( idealeq(H,P) )
1605: break;
1606: F = wellsep(P,V,VL);
1607: GR(P,cons(F,P),VL,PRIMEORD);
1608: }
1609: return R;
1610: }
1611:
1612: /* output M.A.D. sets. AL : M.S.I.S sets */
1613: def misset(AL,VL)
1614: {
1615: for (L=[],D=0,I=length(AL)-1;I>=0;I--) {
1616: E = length(AL[I]);
1617: C = listminus(VL,AL[I]);
1618: if ( E == D )
1619: L = cons(C,L);
1620: else if ( E > D ) {
1621: L = [C]; D = E;
1622: }
1623: }
1624: return L;
1625: }
1626:
1627: end$
1628: