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