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