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