![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to primdec_mod 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_mod CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 noro 1: extern Hom,GBTime$
2: extern DIVLIST,INTIDEAL,ORIGINAL,ORIGINALDIMENSION,STOP,Trials,REM$
3: extern T_GRF,T_INT,T_PD,T_MP$
4: extern BuchbergerMinipoly,PartialDecompByLex,ParallelMinipoly$
5: extern B_Win,D_Win$
6: extern COMMONCHECK_SF,CID_SF$
1.2 ! noro 7: extern FFF_LOADED_BY_PRIMDEC_MOD$
! 8:
! 9: if(FFF_LOADED_BY_PRIMDEC_MOD) load("fff"); else ;
! 10: FFF_LOADED_BY_PRIMDEC_MOD = 1$
1.1 noro 11:
12: /*==============================================*/
13: /* prime decomposition of ideals over */
14: /* finite fields */
15: /* part 1 */
16: /* radical computation */
17: /* 28/12/2002 for Basic Small Finite Field */
18: /* 4/1/2003 for Arbitary Small Finite Field */
19: /*==============================================*/
20:
21: /*==============================================*/
22: /* radical computation of ideals over */
23: /* finite fields by Matsumoto's algorithm */
24: /* */
25: /* The radical of I is computed as the */
26: /* kernel of a power of Frobenius map. */
27: /* */
28: /* radical */
29: /* Input: a list P of polynomials */
30: /* Output: a list P1 of polynomials */
31: /* such that <P1>=rad(<P>) */
32: /* */
33: /* frobeniuskernel_main{2} */
34: /* Input: a list of polynomials P, */
35: /* a list of variables VSet, */
36: /* and a list of other variables WSet */
37: /* Output: a list of polynomials Q */
38: /* such that Q is the kernel of */
39: /* single Frobenius map: x--> x^q */
40: /* where q is the characteristic. */
41: /* version 2 uses successive kernel */
42: /* computation suggested by Matsumoto */
43: /* */
44: /* coefficientfrobeniuskernel */
45: /* Input: a list of polynomials P, */
46: /* Output: a list of polynomials Q */
47: /* such that Q is the kernel of */
48: /* single coefficient Frobenius */
49: /* map: a in GF(q) --> a^q */
50: /* where q is the characteristic. */
51: /* */
52: /* elimination */
53: /* Input: a list P of polynomials */
54: /* a list V of variables */
55: /* Output: a list P1 of polynomials */
56: /* such that P1={f\in P | */
57: /* vars(f)\cap V =\emptyset} */
58: /* */
59: /* checkequality */
60: /* Input: a list P of polynomials */
61: /* a list Q of polynomials */
62: /* Output: 1 if <P>=<Q> */
63: /* 0 otherwise */
64: /* */
65: /* */
66: /* */
67: /*==============================================*/
68:
69: def radicalideal(P,Choice,VSet)
70: {
71:
72: GBTime=0;
73: Ord=0;
74:
75: /* VSet=vars(P);*/
76: WSet=makecounterpart(VSet);
77:
78: T000=time()[0];
79: P1=dp_gr_f_main(P,VSet,Hom,Ord);
80: T001=time()[0];
81: GBTime +=T001-T000;
82:
83: if ( Choice == 3 )
84: {
85: P2=frobeniuskernel_main3(P1,VSet,WSet);
86: }
87: else if ( Choice == 2 )
88: {
89: P2=frobeniuskernel_main2(P1,VSet,WSet);
90: }
91: else if ( Choice == 4 )
92: {
93: P2=frobeniuskernel_main4(P1,VSet,WSet);
94: }
95: else
96: {
97: P2=frobeniuskernel_main(P1,VSet,WSet);
98: }
99:
100: M=checkequality(P1,P2,VSet);
101:
102: while ( M !=1 )
103: {
104: P1=P2;
105: P2=frobeniuskernel_main2(P2,VSet,WSet);
106: M=checkequality(P1,P2,VSet);
107: }
108:
109: return P1;
110: }
111:
112:
113: def frobeniuskernel_main(P,VSet,WSet)
114: {
115: NV=length(VSet);
116:
117: NewP=coefficientfrobeniuskernel(P);
118:
119: NW=length(NewP);
120:
121: TempP=NewP;
122:
123: for (K=NW-1;K>=0;K--)
124: {
125: Poly=TempP[K];
126: for (J=0;J<NV;J++)
127: {
128: Poly=subst(Poly,VSet[J],WSet[J]);
129: }
130: NewP=cons(Poly,NewP);
131: }
132:
133: XSet=append(VSet,WSet);
134: NewOrder=[[0,length(VSet)],[0,length(WSet)]];
135:
136: Char=setmod_ff()[0];
137:
138: for (I=0;I<NV;I++)
139: {
140: AddPoly=VSet[I]^Char-WSet[I];
141: NewP=cons(AddPoly,NewP);
142: }
143:
144: G=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
145: PW=elimination(G,WSet);
146:
147: NW=length(PW);
148: ANS=[];
149: for (I=NW-1;I>=0;I--)
150: {
151: Poly=PW[I];
152: for (J=0;J<NV;J++)
153: {
154: Poly=subst(Poly,WSet[J],VSet[J]);
155: }
156: ANS=cons(Poly,ANS);
157: }
158: return ANS;
159: }
160:
161: def frobeniuskernel_main2(P,VSet,WSet)
162: {
163: NV=length(VSet);
164:
165: NewP=coefficientfrobeniuskernel(P);
166:
167: XSet=append(VSet,WSet);
168: NewOrder=[[0,NV],[0,NV]];
169:
170: Char=setmod_ff()[0];
171:
172: for (I=0;I<NV;I++)
173: {
174: for (J=0;J<NV;J++)
175: {
176: if ( I == J )
177: {
178: AddPoly=VSet[J]^Char-WSet[J];
179: }
180: else
181: {
182: AddPoly=VSet[J]-WSet[J];
183: }
184: NewP=cons(AddPoly,NewP);
185: }
186:
187: T000=time()[0];
188: GP=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
189: T001=time()[0];
190: GBTime += T001-T000;
191:
192: PW=elimination(GP,WSet);
193: NW=length(PW);
194: NewP=[];
195: for (K=NW-1;K>=0;K--)
196: {
197: Poly=PW[K];
198: for (J=0;J<NV;J++)
199: {
200: Poly=subst(Poly,WSet[J],VSet[J]);
201: }
202: NewP=cons(Poly,NewP);
203: }
204: }
205: return NewP;
206: }
207:
208:
209: def frobeniuskernel_main4(P,VSet,WSet)
210: {
211: NV=length(VSet);
212:
213: NewP=coefficientfrobeniuskernel(P);
214:
215: XSet=append(VSet,WSet);
216: NewOrder=[[0,NV],[0,NV]];
217:
218: Char=setmod_ff()[0];
219:
220: for (I=0;I<NV;I++)
221: {
222: NW=length(NewP);
223: TempP=NewP;
224:
225: for (J=0;J<NV;J++)
226: {
227: if ( I == J )
228: {
229: AddPoly=VSet[J]^Char-WSet[J];
230: }
231: else
232: {
233: AddPoly=VSet[J]-WSet[J];
234: }
235: NewP=cons(AddPoly,NewP);
236: }
237:
238: /* for (K=NW-1;K>=0;K--)
239: {
240: Poly=TempP[K];
241: for (J=0;J<NV;J++)
242: {
243: if ( J == I )
244: {
245: Poly=subst(Poly,VSet[J],WSet[J]);
246: }
247: else
248: {
249: Poly=subst(Poly,VSet[J],WSet[J]^Char);
250: }
251: }
252: NewP=cons(Poly,NewP);
253: }*/
254: T000=time()[0];
255: GP=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
256: T001=time()[0];
257: GBTime += T001-T000;
258:
259: PW=elimination(GP,WSet);
260: NW=length(PW);
261: NewP=[];
262: for (K=NW-1;K>=0;K--)
263: {
264: Poly=PW[K];
265: for (J=0;J<NV;J++)
266: {
267: Poly=subst(Poly,WSet[J],VSet[J]);
268: }
269: NewP=cons(Poly,NewP);
270: }
271: }
272: return NewP;
273: }
274:
275: def frobeniuskernel_main3(P,VSet,WSet)
276: {
277: NV=length(VSet);
278:
279: NewP=coefficientfrobeniuskernel(P);
280:
281: Char=setmod_ff()[0];
282:
283: for (I=0;I<NV;I++)
284: {
285: XWSet=[];
286: for (J=0;J<NV;J++)
287: {
288: if ( I == J )
289: {
290: AddPoly=VSet[I]^Char-WSet[I];
291: NewP=cons(AddPoly,NewP);
292: XWSet=append(XWSet,[WSet[I]]);
293: }
294: else
295: {
296: XWSet =append(XWSet,[VSet[J]]);
297: }
298: }
299:
300: XSet=cons(VSet[I],XWSet);
301:
302: NewOrder=[[0,1],[0,NV]];
303: T000=time()[0];
304: GP=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
305: T001=time()[0];
306: GBTime += T001-T000;
307:
308: PW=elimination(GP,XWSet);
309: NW=length(PW);
310: NewP=[];
311: for (K=NW-1;K>=0;K--)
312: {
313: Poly=PW[K];
314: Poly=subst(Poly,WSet[I],VSet[I]);
315: NewP=cons(Poly,NewP);
316: }
317: }
318: return NewP;
319: }
320:
321: def coefficientfrobeniuskernel(P)
322: {
323: if ( setmod_ff()[1] == 0 )
324: {
325: return P;
326: }
327:
328: NP=length(P);
329: ANS=[];
330: for (I=0;I<NP;I++)
331: {
332: TEST=P[I];
333: Q=coefficientfrobeniuskernel_main(TEST);
334: ANS=append(ANS,[Q]);
335: }
336: return ANS;
337: }
338:
339: def coefficientfrobeniuskernel_main(Poly)
340: {
341: Vars=vars(Poly);
342: QP=dp_ptod(Poly,Vars);
343: ANS=0;
344: FOrd=deg(setmod_ff()[1],x);
345: Char=setmod_ff()[0];
346: Pow=Char^(FOrd-1);
347:
348: while(QP !=0 )
349: {
350: Coef=dp_hc(QP);
351: HT=dp_ht(QP);
352:
353: ANS=ANS+Coef^Pow*dp_dtop(HT,Vars);
354:
355: QP=dp_rest(QP);
356: }
357: return ANS;
358: }
359:
360:
361: def elimination(G,V)
362: {
363: ANS=[];
364: NG=length(G);
365:
366: for (I=NG-1;I>=0;I--)
367: {
368: VSet=vars(G[I]);
369: DIFF=setminus(VSet,V);
370:
371: if ( DIFF ==[] )
372: {
373: ANS=cons(G[I],ANS);
374: }
375: }
376: return ANS;
377: }
378:
379: def setminus(A,B)
380: {
381: NA=length(A);
382: NB=length(B);
383:
384: ANS=[];
385:
386: for (I=0;I<NA;I++)
387: {
388: Flag =0;
389: for (J=0;J<NB;J++)
390: {
391: if (A[I] == B[J])
392: {
393: Flag=1;
394: break;
395: }
396: }
397:
398: if ( Flag == 0 )
399: {
400: ANS=append(ANS,[A[I]]);
401: }
402: }
403: return ANS;
404: }
405:
406: def makecounterpart(V)
407: {
408: NV=length(V);
409:
410: A="tmp";
411:
412: ANS=[];
413: for (I=NV-1;I>=0;I--)
414: {
415: T=strtov(A+rtostr(I));
416: ANS=cons(T,ANS);
417: }
418: return ANS;
419: }
420:
421: def checkequality(P,Q,VSet)
422: {
423: QQ=dp_gr_f_main(Q,VSet,Hom,Ord);
424: PP=dp_gr_f_main(P,VSet,Hom,Ord);
425:
426: NP=length(PP);
427: NQ=length(QQ);
428:
429: VarsP=vars(P);
430: VarsQ=vars(Q);
431: if ( NP != NQ || VarsP !=VarsQ )
432: {
433: return 0;
434: }
435:
436: for (I=0;I<NP;I++)
437: {
438: HCP=dp_hc(dp_ptod(PP[I],VSet));
439: HCQ=dp_hc(dp_ptod(QQ[I],VSet));
440:
441: if ( HCP*QQ[I]-HCQ*PP[I] != 0 )
442: {
443: return 0;
444: }
445: }
446:
447: return 1;
448: }
449:
450: /*==============================================*/
451: /* prime decomposition of ideals over */
452: /* finite fields */
453: /* part 2 */
454: /* ideal dimension */
455: /* 3/1/2003 */
456: /*==============================================*/
457:
458: /*==============================================*/
459: /* ideal dimension computation */
460: /* */
461: /* The dimension of I is computed as the */
462: /* maximal size of MSI sets. */
463: /* */
464: /* idealdimension */
465: /* Input: a list P of polynomials */
466: /* a list V of variables */
467: /* Output: the dimension of I */
468: /* and the MSI with max size */
469: /* */
470: /* findmsi */
471: /* Input: a list of head monomials */
472: /* and a list of variables */
473: /* Output: a list of all MSI sets */
474: /* and SI sets */
475: /* */
476: /* findmsi_main */
477: /* Input: a list P of polynomials */
478: /* a list V of variables */
479: /* Output: a list P1 of polynomials */
480: /* such that P1={f\in P | */
481: /* vars(f)\cap V =\emptyset} */
482: /* */
483: /* headtermset */
484: /* Input: a list P of polynomials */
485: /* a list V of variables */
486: /* Output: a list of head monomials */
487: /* */
488: /* */
489: /* */
490: /*==============================================*/
491:
492: def idealdimension(P,V,Ord)
493: {
494: GP=dp_gr_f_main(P,V,Hom,Ord);
495: HeadTermset=headtermset(GP,V,Ord);
496: MSIset=findmsi(HeadTermset,V);
497: if ( MSIset == [] )
498: {
499: return [0,[]];
500: }
501: Index=findmaxindex(MSIset);
502: MSI=MSIset[Index];
503: Dimension=length(MSI);
504: return [Dimension, MSI];
505: }
506:
507: def findmaxindex(S)
508: {
509: NS=length(S);
510: Index=0;
511: Defaultsize=length(S[0]);
512: for (I=1;I<NS;I++)
513: {
514: Targetsize=length(S[I]);
515:
516: if ( Targetsize > Defaultsize )
517: {
518: Index=I;
519: }
520: Defaultsize= Targetsize;
521: }
522: return Index;
523: }
524:
525: def headtermset(P,V,Ord)
526: {
527: ANS=[];
528: NP=length(P);
529: for ( I=NP-1;I>=0;I--)
530: {
531: Head=headterm(P[I],V,Ord);
532: ANS=cons(Head,ANS);
533: }
534: return ANS;
535: }
536:
537: def headterm(P,V,Ord)
538: {
539: dp_ord(Ord);
540: Q=dp_ptod(P,V);
541: Headdp=dp_hm(Q);
542: Head=dp_dtop(Headdp,V);
543: return Head;
544: }
545:
546: def findmsi_main(TermsetIndex,MSIsetIndex,Int,N)
547: {
548: ANS=[];
549: BNS=[];
550: TempMSIsetIndex=MSIsetIndex;
551:
552: if (Int == 0)
553: {
554: for ( I=0;I<N;I++)
555: {
556: TEST1=[I];
557: Check=intersection(TermsetIndex,TEST1);
558: if ( Check == 0 )
559: {
560: ANS=cons(TEST1,ANS);
561: }
562: }
563: }
564:
565: while( Int !=0 && TempMSIsetIndex != [] )
566: {
567: TEST=TempMSIsetIndex[0];
568: Flag=0;
569: for ( I=TEST[0]+1;I<N;I++)
570: {
571: TEST1=cons(I,TEST);
572: Check=intersection(TermsetIndex,TEST1);
573: if ( Check == 0 )
574: {
575: Flag=1;
576: ANS=cons(TEST1,ANS);
577: }
578: }
579: if ( Flag == 0 )
580: {
581: BNS=cons(TEST,BNS);
582: }
583: TempMSIsetIndex=cdr(TempMSIsetIndex);
584: }
585: return [ANS,BNS];
586: }
587:
588: def findmsi(Termset,Vset)
589: {
590: ANS=[];
591: MSIsetIndex=[];
592: N=length(Vset);
593: TermsetIndex=makeindex(Termset,Vset);
594:
595: for (I=0;I<N;I++)
596: {
597: Temp=findmsi_main(TermsetIndex,ANS,I,N);
598: ANS=Temp[0];
599: BNS=Temp[1];
600: MSIsetIndex=append(BNS,MSIsetIndex);
601: }
602:
603: MSIsetIndex=append(ANS,MSIsetIndex);
604: MSIset=makevarset(MSIsetIndex,Vset);
605: return MSIset;
606: }
607:
608: def makeindex(Termset,Vset)
609: {
610: ANS=[];
611:
612: N=length(Vset);
613:
614: TempTermset=Termset;
615: while( TempTermset !=[] )
616: {
617: TEST=TempTermset[0];
618: Index=[];
619: for (I=0;I<N;I++)
620: {
621: Q=srem(TEST,Vset[I]);
622: if ( Q == 0 )
623: {
624: Index=cons(I,Index);
625: }
626: }
627: ANS=cons(Index,ANS);
628: TempTermset=cdr(TempTermset);
629: }
630: return ANS;
631: }
632:
633: def makevarset(IndexSet,Vset)
634: {
635: ANS=[];
636: NI=length(IndexSet);
637:
638: for (I=0;I<NI;I++)
639: {
640: TEST=IndexSet[I];
641:
642: TestV=makevar(TEST,Vset);
643: ANS=append(ANS,[TestV]);
644: }
645: return ANS;
646: }
647:
648: def makevar(Index,Vset)
649: {
650: ANS=[];
651: TempIndex=Index;
652:
653: while( TempIndex !=[] )
654: {
655: ANS=cons(Vset[TempIndex[0]],ANS);
656: TempIndex=cdr(TempIndex);
657: }
658:
659: return ANS;
660: }
661:
662: def intersection(IndexSet,TestIndex)
663: {
664: TempIndexSet=IndexSet;
665:
666: while(TempIndexSet !=[] )
667: {
668: Sample=TempIndexSet[0];
669:
670: Inclusion=testinclusion(Sample,TestIndex);
671:
672: if ( Inclusion == 1 )
673: {
674: return 1;
675: }
676:
677: TempIndexSet=cdr(TempIndexSet);
678: }
679: return 0;
680: }
681:
682: def testinclusion(Sample,Test)
683: {
684: NS=length(Sample);
685: NT=length(Test);
686:
687: for (I=0;I<NS;I++)
688: {
689: Flag=0;
690: for (J=0;J<NT;J++)
691: {
692: if (Sample[I]==Test[J])
693: {
694: Flag=1;
695: break;
696: }
697: }
698: if ( Flag == 0 )
699: {
700: return 0;
701: }
702: }
703: return 1;
704: }
705:
706:
707: /*======================================================*/
708: /* prime decomposition of ideals over */
709: /* small finite fields */
710: /* part 3 */
711: /* prime decomposition */
712: /* 9/1/2003 1st Implementation */
713: /* 11/1/2003 Divid 3 cases in zeroprime */
714: /* decomposition (3a) */
715: /* 12/1/2003 Redundant elimination (3e) */
716: /* 13/1/2003 gr_fctr_sf version (3f) */
717: /* 14/1/2003 Flag on Early Termination */
718: /* 15/1/2002 drl version */
719: /* 17/1/2002 depth oriented search */
720: /* 17/1/2002 dimension oriented search */
721: /* */
722: /* global variables */
723: /* Hom: used in dp_gr_f_main */
724: /* DIVLIST: the set of prime ideals */
725: /* ORIGINAL: the GB of the original input ideal */
726: /* ORIGINALDIMENSION: the dimension */
727: /* STOP: the flag on terminatation */
728: /* Trials: a bound on the number of trials */
729: /* REM: the vector list of remaining ideals */
730: /* */
731: /*======================================================*/
732:
733: /*======================================================*/
734: /* prime decomposition */
735: /* */
736: /* */
737: /* primedec(P,V,Ord,F) */
738: /* Input: a list P of polynomials */
739: /* a list V of variables */
740: /* a term order Ord */
741: /* a flag F on Early Termination */
742: /* Output: 0 */
743: /* DVILIST: the set of prime divisors of <P> */
744: /* [subprocedure] */
745: /* o --gr_fctr_sf (mplex2) */
746: /* o --primedecomposition */
747: /* Remark: if F = 1, we empoly Early Termination */
748: /* */
749: /* */
750: /* primedecomposition(P,V,Ord,C,F) */
751: /* Input: a list P of polynomials */
752: /* a list V of variables */
753: /* a term order Ord */
754: /* a level C of recursive call */
755: /* a flag F on Early Termination */
756: /* Output: 0 */
757: /* DVILIST: the set of prime divisors of <P>^ec */
758: /* [subprocedure] */
759: /* o --idealdimension (part2) */
760: /* o --setminus (part1) */
761: /* o --zeroprimedecomposition */
762: /* o --extcont */
763: /* o --checkadd2 */
764: /* o --ideal_intersection_sfrat */
765: /* o --radical_equality */
766: /* */
767: /* zeroprimedecomposition(P,V,W) */
768: /* Input: a list P of polynomials */
769: /* lists V,W of variables */
770: /* such that <P> is 0-dimensional */
771: /* in K[V] and W is the original set */
772: /* of variables */
773: /* Output: the set of prime/primary */
774: /* divisors of <P> */
775: /* (P is a GB w.r.t. DRL) */
776: /* [subprocedure] */
777: /* o --partical_decomp (mplex2) */
778: /* o --checkgeneric */
779: /* o --separableclosure */
780: /* o --zerogenericprimedecomposition */
781: /* o --zeroseparableprimedecomposition */
782: /* o --convertdivisor */
783: /* o --contraction */
784: /* o --radicalideal (part1) */
785: /* */
786: /* extcont(P,V,Ord) */
787: /* Input: a list P of polynomials */
788: /* a set V of variables */
789: /* Output: a polynomial f */
790: /* such that \sqrt<P>^c=\sqrt(P:f) */
791: /* Then \sqrt<P>=\sqrt(P:f)\cap \sqrt(P+<f>) */
792: /* */
793: /* separableclosure([P,MP],V,W) */
794: /* Input: a pair [P,MP] of a list P of */
795: /* polynomials and a list MP of */
796: /* minimal polynomials of variables */
797: /* list V,W of variables */
798: /* such that <P> is 0-dimensional */
799: /* in K[V] and W is the original set */
800: /* of variables */
801: /* Output: [Q,E], a list Q of polynomials */
802: /* an exponent vector E */
803: /* such that <Q> is the separable */
804: /* closure of <P> */
805: /* [subprocedure] */
806: /* o --makecounterpart(part1) */
807: /* o --elimination(part1) */
808: /* o --checkseparable */
809: /* */
810: /* zeroseparableprimedecomposition(P,V,W) */
811: /* Input: a list P of polynomials */
812: /* lists V,W of variables */
813: /* such that <P> is 0-dimensional */
814: /* and all minimal polynomial */
815: /* are irreducible (i.e. <P> is an */
816: /* intermediate ideal in the paper */
817: /* or the output of partial_decomp) */
818: /* in K[V] and W is the original set */
819: /* of variables */
820: /* Output: the set of prime/primary */
821: /* divisors of <P> */
822: /* [subprocedure] */
823: /* o --findgeneric */
824: /* o --elimination(part1) */
825: /* */
826: /* zerogenericprimedecomposition([P,MP],M,V,W) */
827: /* Input: a pair [P,MP] of a list P of */
828: /* polynomials and a list MP of */
829: /* minimal polynomials of variables */
830: /* such that <P> is 0-dimensional */
831: /* and M is a minimal poly of a */
832: /* variable in generic position */
833: /* in K[V] and W is the original */
834: /* set of variables */
835: /* Output: the set of prime/primary */
836: /* divisors of <P> */
837: /* [subprocedure] */
838: /* */
839: /* convertdivisor(P,V,W,E) */
840: /* Input: a list P of polynomials */
841: /* list V,W of variables */
842: /* such that <P> is a primary/prime */
843: /* 0-dimensional ideal in K[V] and */
844: /* W is the original set */
845: /* of variables */
846: /* the exponent vector E */
847: /* Output: the corresponding prime ideal */
848: /* in K[W] */
849: /* [subprocedure] */
850: /* o --elimination(part1) */
851: /* */
852: /* contraction(P,V,W) */
853: /* Input: a list P of polynomials */
854: /* list V,W of variables */
855: /* such that <P> is a primary/prime */
856: /* 0-dimensional ideal in K[V] and */
857: /* W is the original set */
858: /* of variables */
859: /* Output: the contraction <P>^c in K[W] */
860: /* */
861: /* findgeneric(P,V) */
862: /* Input: a list P of polynomials */
863: /* a list V of variables */
864: /* such that <P> is a separable */
865: /* 0-dimensional ideal in K[V] */
866: /* Output: [F,Min,X], a polynomial F in */
867: /* generic position, its minimal */
868: /* polynomial Min in a new variable */
869: /* X */
870: /* [subprocedure] */
871: /* o --lineardimension */
872: /* o --makemagic */
873: /* o --minipoly_sf (mplex2) */
874: /* */
875: /*======================================================*/
876:
877: def primedec_mod(P,VSet,Ord,Mod,Strategy)
878: {
879: for ( Q = Mod, E = 1; Q < 2^14; Q *= Mod, E++ );
880: Q /= Mod;
881: E--;
882: if ( !(E%2) ) {
883: Q /= Mod;
884: E--;
885: }
886: setmod_ff(Mod,E);
887: Pq = map(simp_ff,P);
888: primedec_sf(Pq,VSet,Ord,Strategy);
889: R = convsf(DIVLIST,VSet,0,1);
890: return R;
891: }
892:
893:
894: def primedec_sf(P,VSet,Ord,Strategy)
895: {
896: /* DIVLIST = the set of all computed candidates for iredundant divisors */
897: /* INTIDEAL = the intersection of all computed candidates */
898: /* ORIGINAL = a Groebner basis of the input ideal <P> w.r.t. Ord */
899: /* ORIGINALDIMENSION = the dimension of the input ideal <P> */
900: /* STOP = the flag on termination */
901: /* REM = the list of remaining ideals */
902:
903: T0 = time();
904: T_GRF=0;
905: T_INT=0;
906: T_PD=0;
907: T_MP=0;
908: DIVLIST=[];
909: INTIDEAL=[];
910: STOP=0;
911: B_Win=0;
912: D_Win=0;
913:
914: ORIGINAL=dp_gr_f_main(P,VSet,Hom,Ord);
915:
916: Dimeset=idealdimension(ORIGINAL,VSet,Ord);
917: ORIGINALDIMENSION=Dimeset[0];
918:
919: /* REM0=newvect(ORIGINALDIMENSION+1);*/
920: REM=newvect(ORIGINALDIMENSION+1);
921:
922: for (I=0;I<ORIGINALDIMENSION+1;I++)
923: {
924: /* REM0[I]=[];*/
925: REM[I]=[];
926: }
927:
928: print("The dimension of the ideal is ",2);print(ORIGINALDIMENSION,2);
929: print(". ");
930:
931: if ( ORIGINALDIMENSION == 0 )
932: {
933: Strategy = 0;
934: STOP = 0;
935: }
936:
937: ANS=gr_fctr_sf([ORIGINAL],VSet,Ord);
938: NANS=length(ANS);
939: print("There are ",2);print(NANS,2);print(" partial components. ");
940:
941: for (I=0;I<NANS;I++)
942: {
943: TempI=ANS[I];
944: DimI=idealdimension(TempI,VSet,Ord)[0];
945:
946: REM[ORIGINALDIMENSION-DimI]=cons(TempI,REM[ORIGINALDIMENSION-DimI]);
947: }
948:
949: REM[0]=ANS;
950:
951: INDEX=0;
952:
953: while ( INDEX != -1 )
954: {
955: /* print("We deal with the ",2);print(I,2);print("-th partial component. ");*/
956:
957: TESTIDEAL=car(REM[INDEX]);
958: REM[INDEX]=cdr(REM[INDEX]);
959: primedecomposition(TESTIDEAL,VSet,Ord,INDEX,Strategy);
960:
961: if (STOP == 1 )
962: {
963: DIVLIST = prime_irred_sf_by_first(DIVLIST,VSet,0);
964: DIVLIST = monic_sf_first(DIVLIST,VSet);
965: print("We finish the computation. ");
966: T_TOTAL = time()[3]-T0[3];
967: print(["T_TOTAL",T_TOTAL,"T_GRF",T_GRF,"T_PD",T_PD,"T_MP",T_MP,"T_INT",T_INT,"B_Win",B_Win,"D_Win",D_Win]);
968: return 0;
969: }
970:
971: INDEX=findindex(REM);
972: }
973:
974: DIVLIST = prime_irred_sf_by_first(DIVLIST,VSet,0);
975: DIVLIST = monic_sf_first(DIVLIST,VSet);
976: T_TOTAL = time()[3]-T0[3];
977: print(["T_TOTAL",T_TOTAL,"T_GRF",T_GRF,"T_PD",T_PD,"T_MP",T_MP,"T_INT",T_INT,"B_Win",B_Win,"D_Win",D_Win]);
978: return 0;
979: }
980:
981:
982:
983:
984:
985: def findindex(VECTORLIST)
986: {
987: NL=size(VECTORLIST)[0];
988:
989: for (I=0;I<NL;I++)
990: {
991: if ( VECTORLIST[I] == [] )
992: {
993: continue;
994: }
995: return I;
996: }
997: return -1;
998: }
999:
1000: #if 0
1001: def checkadd2(ADD,TESTADDLIST,VSet)
1002: {
1003: /* This function will be used for eliminating redundant divisors */
1004:
1005: NTESTADDLIST=length(TESTADDLIST);
1006: for (I=0;I<NTESTADDLIST;I++)
1007: {
1008: TESTLIST=TESTADDLIST[I][0];
1009:
1010: if ( setminus(TESTLIST,ADD) == [] )
1011: {
1012: return 1;
1013: }
1014: if ( inclusion_test(TESTLIST,ADD,VSet,Ord) == 1 )
1015: {
1016: return 1;
1017: }
1018: }
1019: return 0;
1020: }
1021: #endif
1022:
1023: def checkadd2a(ADD,DIM,TESTADDLIST,VSet)
1024: {
1025: /* This function will be used for eliminating redundant divisors */
1026:
1027: NTESTADDLIST=length(TESTADDLIST);
1028: for (I=0;I<NTESTADDLIST;I++)
1029: {
1030: TESTLIST=TESTADDLIST[I][0];
1031: TESTDIM=TESTADDLIST[I][1];
1032:
1033: if (DIM > TESTDIM )
1034: {
1035: continue;
1036: }
1037:
1038: if ( setminus(TESTLIST,ADD) == [] )
1039: {
1040: return 1;
1041: }
1042: if ( inclusion_test(TESTLIST,ADD,VSet,Ord) == 1 )
1043: {
1044: return 1;
1045: }
1046: }
1047: return 0;
1048: }
1049:
1050: def checkadd3(ADD,TESTADDLIST,VSet)
1051: {
1052: /* This function will be used for eliminating redundant divisors */
1053:
1054: NTESTADDLIST=length(TESTADDLIST);
1055: for (I=0;I<NTESTADDLIST;I++)
1056: {
1057: TESTLIST=TESTADDLIST[I];
1058:
1059: if ( setminus(TESTLIST,ADD) == [] )
1060: {
1061: return 1;
1062: }
1063:
1064:
1065: /* if ( inclusion_test(TESTLIST,ADD,VSet,0) == 1 )
1066: {
1067: return 1;
1068: }*/
1069: }
1070: return 0;
1071: }
1072:
1073: def radical_equality(A,B,VSet,Ord)
1074: {
1075: /* This function will be used for checking the termination. */
1076:
1077: NA=length(A);
1078:
1079: Ord1=[[0,1],[0,length(VSet)]];
1080: Vt=newt;
1081:
1082: for (I=0;I<NA;I++)
1083: {
1084: NewB=cons(1-newt*A[I],B);
1085: GB = dp_gr_f_main(NewB,cons(Vt,VSet),0,Ord1);
1086: K=elimination(GB,VSet);
1087:
1088: if ( vars(K) !=[] )
1089: {
1090: return 0;
1091: }
1092:
1093: }
1094:
1095: return 1;
1096: }
1097:
1098: def primedecomposition(P,VSet,Ord,COUNTER,Strategy)
1099: {
1100: /* DIVLIST = the set of all computed candidates for iredundant divisors */
1101: /* INTIDEAL = the intersection of all computed candidates */
1102: /* ORIGINAL = a Groebner basis of the input ideal <P> w.r.t. Ord */
1103: /* ORIGINALDIMENSION = the dimension of the input ideal <P> */
1104:
1105: RP=P;
1106: GP=dp_gr_f_main(RP,VSet,Hom,Ord);
1107:
1108: /* First we compute the MSI and the dimension of the */
1109: /* ideal <P> in K[Vset], where K is the ground field. */
1110:
1111: Dimeset=idealdimension(GP,VSet,Ord);
1112: Dimension=Dimeset[0];
1113: MSI=Dimeset[1];
1114:
1115: print("The dimension of the ideal is ",2); print(Dimension,2);
1116: print(".");
1117:
1118: TargetVSet=setminus(VSet,MSI);
1119: NewGP=dp_gr_f_main(GP,TargetVSet,Hom,Ord);
1120:
1121: if ( vars(NewGP) == [] )
1122: {
1123: return [];
1124: }
1125:
1126: /* Then the ideal is 0-dimension in K[TargetVSet]. */
1127:
1128: print("We enter Zero-dimension Prime Decomposition. ",2);
1129:
1130: QP=zeroprimedecomposition(NewGP,TargetVSet,VSet);
1131:
1132: ANS=[];
1133: NQP=length(QP);
1134:
1135: print("The number of the newly found component is ",2);
1136: print(NQP,2);print(". ",2);
1137:
1138: for (I=0;I<NQP;I++)
1139: {
1140: ZPrimeideal=QP[I];
1141: Primedivisor=ZPrimeideal;
1142: CHECKADD=checkadd2a(Primedivisor,Dimension,DIVLIST,VSet);
1143: if (CHECKADD == 0 )
1144: {
1145: DIVLIST=append(DIVLIST,[[Primedivisor,Dimension]]);
1146:
1147: if (Strategy != 1 )
1148: {
1149: /* NO-OPERATION */
1150: }
1151: else if ( COUNTER == 0 && Dimension == 0 && ORIGINALDIMENSION == 0 )
1152: {
1153: /* NO-OPERATION */
1154: }
1155: else if ( INTIDEAL == [] )
1156: {
1157: INTIDEAL=Primedivisor;
1158: }
1159: else
1160: {
1161: INTIDEAL=ideal_intersection_sfrat(Primedivisor,INTIDEAL,VSet);
1162: }
1163: }
1164: }
1165:
1166: /* We compute prime decomposition for remaining ideal */
1167: /* I+<ExtPoly> by recursive call. */
1168:
1169: if ( Strategy == 1 )
1170: {
1171: CHECK=radical_equality(INTIDEAL,ORIGINAL,VSet,Ord);
1172:
1173: if (CHECK==1)
1174: {
1175: print("We already obtain all divisor. ");
1176: STOP = 1;
1177: return 0;
1178: }
1179: }
1180:
1181:
1182: if ( Dimension == 0 )
1183: {
1184: return 0;
1185: }
1186:
1187: Ord1 = [[0,length(TargetVSet)],[0,length(MSI)]];
1188: GP1 = dp_gr_f_main(GP,append(TargetVSet,MSI),Hom,Ord1);
1189: ExtpolyFactor=extcont_factor(GP1,TargetVSet,0);
1190:
1191: COUNTER=COUNTER+1;
1192:
1193: for ( Tmp = ExtpolyFactor; Tmp != []; Tmp = cdr(Tmp) )
1194: {
1195: AddPoly=car(Tmp);
1196:
1197: NewGP=cons(AddPoly,GP);
1198:
1199: GP1 = dp_gr_f_main(cons(AddPoly,GP),VSet,Hom,Ord);
1200:
1201: if ( Strategy == 1 )
1202: {
1203: CHECKADD=radical_equality(INTIDEAL,NewGP,VSet,Ord);
1204:
1205: if ( CHECKADD != 0 )
1206: {
1207: print("Avoid unnecessary computation. ",2);
1208: continue;
1209: }
1210: }
1211:
1212:
1213: DimI=idealdimension(GP1,VSet,Ord)[0];
1214:
1215: if ( REM[ORIGINALDIMENSION-DimI] !=[] )
1216: {
1217: REM[ORIGINALDIMENSION-DimI]=cons(GP1,REM[ORIGINALDIMENSION-DimI]);
1218: }
1219: else
1220: {
1221: REM[ORIGINALDIMENSION-DimI]=[GP1];
1222: }
1223:
1224: /* primedecomposition(GP1,VSet,Ord,COUNTER,Strategy);
1225: if ( STOP == 1)
1226: {
1227: return 0;
1228: }*/
1229: }
1230: return 0;
1231: }
1232:
1233: /* returns 1 if Poly is a subset of Id(GB) */
1234: /* we assume GB is a gb wrt (V,Ord) */
1235:
1236: def inclusion_test(Poly,GB,V,Ord)
1237: {
1238: Len = length(GB);
1239: PS = newvect(Len);
1240: dp_ord(Ord);
1241: for ( J = 0, T = GB; T != []; T = cdr(T), J++ )
1242: PS[J] = dp_ptod(car(T),V);
1243: for ( J = Len-1, Ind = []; J >= 0; J-- )
1244: Ind = cons(J,Ind);
1245:
1246: for ( T = Poly; T != []; T = cdr(T) )
1247: if ( nf_sfrat(Ind,dp_ptod(car(T),V),1,PS)[0] )
1248: return 0;
1249: return 1;
1250: }
1251:
1252:
1253: def zeroprimedecomposition(P,TargetVSet,VSet)
1254: {
1255: /* We compute prime decomposition for an ideal <P> */
1256: /* such that <P> is 0-dimensional in K[TargetVSet]. */
1257:
1258: PD=partial_decomp(P,TargetVSet);
1259:
1260: /* each ideal in PD is an intermediate ideal, */
1261: /* i.e., all minimal polynomials are irreducible. */
1262: /* Thus, each ideal is very likely a prime ideal. */
1263: /* PD=[[P,MP],...], where P is a GB w.r.t. DRL and */
1264: /* MP is the set of all minimal polynomials. */
1265:
1266: NPD=length(PD);
1267: ANS=[];
1268:
1269: for (I=0;I<NPD;I++)
1270: {
1271: COMP=PD[I];
1272:
1273: /* check if the ideal has a variable in generic position*/
1274:
1275: CHECK=checkgeneric(COMP,TargetVSet,VSet);
1276:
1277: /* If there is a variable in generic position, then we apply a special */
1278: /* procedure (zerogenericprimedecomposition). Otherwise we apply a */
1279: /* general procedure (zeroseparableprimedecomposition). */
1280: /* */
1281: /* CHECK=[1,M], where M is the minimal polynomail of a variable x */
1282: /* such that x is in generic position, if there is such a variable x. */
1283: /* Otherwise, CHECK=0. */
1284:
1285: if (CHECK != 0 )
1286: {
1287: if ( TargetVSet != VSet )
1288: {
1289: PDiv=contraction(COMP[0],TargetVSet,VSet)[0];
1290: }
1291: else
1292: {
1293: PDiv=COMP[0];
1294: }
1295:
1296: ZDecomp=[PDiv];
1297:
1298: print("An intermediate ideal is of generic type. ");
1299: }
1300: else
1301: {
1302: print("An intermediate ideal is not of generic type. ",2);
1303:
1304: /* We compute the separable closure of <P> by using minimal polynomails.*/
1305: /* separableclosure outputs */
1306: /* [a GB Q of separable closure, the exponent vector]. */
1307: /* If <P> is already separable, then the exponent vector is 0. */
1308:
1309: Sep=separableclosure(COMP,TargetVSet,VSet);
1310:
1311: /* Then, we compute prime divisors of the separable closure by using */
1312: /* generic position. */
1313: /* If Sep[1]=0, COMP[0] is already separable ideal, and hence, */
1314: /* we do not apply radicalideal to its divisors. */
1315:
1316: if ( Sep[1] != 0 )
1317: {
1318: print("The ideal is inseparable. ",2);
1319: CHECK2=checkgeneric2(Sep[2]);
1320: }
1321: else
1322: {
1323: print("The ideal is already separable. ",2);
1324: }
1325:
1326: if ( Sep[1] !=0 && CHECK2 == 1 )
1327: {
1328: print("The separable closure is of generic type. ",2);
1329: print("So, the intermediate ideal is prime or primary. ",2);
1330:
1331: PDiv=convertdivisor(Sep[0],TargetVSet,VSet,Sep[1]);
1332: if ( TargetVSet != VSet )
1333: {
1334: PDiv=contraction(PDiv,TargetVSet,VSet)[0];
1335: }
1336: Divisor=radicalideal(PDiv,2,VSet);
1337: ZDecomp=[Divisor];
1338: }
1339: else
1340: {
1341: #if 0
1342: ZSDecomp=zeroseparableprimedecomposition(Sep[0],TargetVSet,VSet);
1343: #else
1344: ZSDecomp=zerosepdec(Sep[0],TargetVSet,VSet,Ord);
1345: #endif
1346: /* We convert a prime separable ideal in K[TargetVSet] to its */
1347: /* corresponding prime ideal in K[VSet]. */
1348:
1349: NComp=length(ZSDecomp);
1350: ZDecomp=[];
1351:
1352: for (J=0;J<NComp;J++)
1353: {
1354: SDiv=ZSDecomp[J];
1355:
1356: /* First we convert a prime separable ideal in K[TargetVSet] to its */
1357: /* corresponding prime/primary ideal in K[TargetVSet] by computing the */
1358: /* image of Frobenius map. */
1359:
1360: PDiv=convertdivisor(SDiv,TargetVSet,VSet,Sep[1]);
1361:
1362: /* Then we compute the true prime divisor by contraction and radical */
1363: /* computation. */
1364:
1365: if ( TargetVSet != VSet )
1366: {
1367: PDiv=contraction(PDiv,TargetVSet,VSet)[0];
1368: }
1369:
1370: if (Sep[1] != 0 )
1371: {
1372: Divisor=radicalideal(PDiv,2,VSet);
1373: }
1374: else
1375: {
1376: Divisor=PDiv;
1377: }
1378:
1379:
1380: ZDecomp=append(ZDecomp,[Divisor]);
1381: }
1382: }
1383: }
1384:
1385: ANS=append(ANS, ZDecomp);
1386: }
1387: return map(monic_hc,ANS,VSet);
1388: }
1389:
1390: def monic_hc(Poly,V)
1391: {
1392: if ( type(Poly) == 4 )
1393: map(monic_hc,Poly,V);
1394: else {
1395: T = dp_ptod(Poly,V);
1396: return dp_dtop(T/dp_hc(T),V);
1397: }
1398: }
1399:
1400: def zeroseparableprimedecomposition(P,TargetVSet,VSet)
1401: {
1402: /* We compute prime decomposition for a separable */
1403: /* 0-dimensional ideal <P> in K[TargetVSet] by using */
1404: /* generic position. */
1405:
1406: ANS=[];
1407: Ord=0;
1408:
1409: NVSet=length(TargetVSet);
1410:
1411: /* The P is already a GB w.r.t. DRL. So, the following */
1412: /* must be replaced with NewP=P. */
1413:
1414: /* NewGP=dp_gr_f_main(P,TargetVSet,Hom,Ord); */
1415: NewGP = P;
1416:
1417: /* We search a polynomial f in generic position. */
1418: /* Generic=[f, minimal polynomial of f in newt, newt], */
1419: /* where newt (X) is a newly introduced variable. */
1420:
1421: print("We search for a linear sum of variables in generic position. ",2);
1422:
1423: Generic=findgeneric(NewGP,TargetVSet,VSet);
1424:
1425: X=Generic[2]; /* newly introduced variable */
1426: Factors=sffctr(Generic[1]);
1427: NFactors=length(Factors);
1428: /* irreducible case */
1429: if ( NFactors == 2 && Factors[1][1] == 1 )
1430: return [NewGP];
1431: #if 0
1432: for (J=1;J<NFactors;J++)
1433: {
1434: AddPoly=Factors[J][0];
1435: AddPoly2=X-Generic[0];
1436:
1437: TestGP=append(NewGP,[AddPoly, AddPoly2]);
1438: VS=cons(X,TargetVSet);
1439: Q=dp_gr_f_main(TestGP,VS,Hom,[[0,1],[Ord,NVSet]]);
1440: QR=elimination(Q,VSet);
1441: ANS=append(ANS,[QR]);
1442: }
1443: #else
1444: /* noro */
1445: dp_ord([[0,1],[Ord,NVSet]]);
1446: XVSet = cons(X,TargetVSet);
1447: PS = newvect(length(NewGP)+1,map(dp_ptod,cons(X-Generic[0],NewGP),XVSet));
1448: for ( I = length(NewGP), Ind = [];I >= 0; I-- )
1449: Ind = cons(I,Ind);
1450: Ind = reverse(Ind);
1451: YSet = setminus(VSet,TargetVSet);
1452: NYSet = length(YSet);
1453: ElimVSet = append(TargetVSet,YSet);
1454: ElimOrd = [[Ord,NVSet],[0,NYSet]];
1455: for ( J = 1; J < NFactors; J++ ) {
1456: dp_ord([[0,1],[Ord,NVSet]]);
1457: Factor = dp_dtop(nf_sfrat(Ind,dp_ptod(Factors[J][0],XVSet),1,PS)[0],XVSet);
1458: #if 0
1459: Q = dp_gr_f_main(cons(Factor,NewGP),ElimVSet,Hom,ElimOrd);
1460: #else
1461: Q0 = dp_gr_f_main(cons(Factor,NewGP),ElimVSet,Hom,0);
1462: Q = dp_gr_f_main(Q0,ElimVSet,Hom,ElimOrd);
1463: #endif
1464: Q = dp_gr_f_main(Q,TargetVSet,Hom,Ord);
1465: ANS = cons(Q,ANS);
1466: }
1467: #endif
1468: return ANS;
1469: }
1470:
1471: /* partial decomposition by sums of variables
1472: -> zeroseparableprimedecomposition */
1473:
1474: #if 0
1475: def zerosepdec(P,W,V,Ord)
1476: {
1477: /* P is GB wrt (W,Ord) */
1478: dp_ord(Ord);
1479: DIM = length(dp_mbase(map(dp_ptod,P,W)));
1480: /* preprocessing */
1481: N = length(W);
1482: C = newvect(N);
1483: C[0] = 1;
1484: Vt = ttttt;
1485: do {
1486: for ( I = 0, LF = 0; I < N; I++ )
1487: if ( C[I] ) LF += W[I];
1488: LF = simp_ff(LF);
1489: MP = minipoly_sf(P,W,Ord,LF,Vt);
1490: Factors = map(first,cdr(sffctr(MP)));
1491: if ( deg(MP,Vt) == DIM ) {
1492: if ( length(Factors) == 1 )
1493: return [P];
1494: else
1495: return zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));
1496: } else if ( length(Factors) > 1 ) {
1497: R = zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));
1498: S = [];
1499: for ( TR = R; TR != []; TR = cdr(TR) )
1500: S = append(zerosepdec(car(TR),W,V,Ord),S);
1501: return S;
1502: }
1503: } while ( !nextchoice(C) );
1504: /* we could not find any useful combination */
1505: R = zeroseparableprimedecomposition(P,W,V);
1506: return R;
1507: }
1508: #else
1509:
1510: /* we only search for useful poly among v[i]+v[j] */
1511: def zerosepdec(P,W,V,Ord)
1512: {
1513: /* P is GB wrt (W,Ord) */
1514: dp_ord(Ord);
1515: DIM = length(dp_mbase(map(dp_ptod,P,W)));
1516: /* preprocessing */
1517: N = length(W);
1518: C = newvect(N);
1519: C[0] = 1;
1520: Vt = ttttt;
1521: for ( I1 = 0; I1 < N; I1++ )
1522: for ( I2 = I1+1; I2 < N; I2++ ) {
1523: LF = simp_ff(W[I1]+W[I2]);
1524: MP = minipoly_sf(P,W,Ord,LF,Vt);
1525: Factors = map(first,cdr(sffctr(MP)));
1526: if ( deg(MP,Vt) == DIM ) {
1527: if ( length(Factors) == 1 )
1528: return [P];
1529: else
1530: return zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));
1531: } else if ( length(Factors) > 1 ) {
1532: R = zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));
1533: S = [];
1534: for ( TR = R; TR != []; TR = cdr(TR) )
1535: S = append(zerosepdec(car(TR),W,V,Ord),S);
1536: return S;
1537: }
1538: }
1539: /* we could not find any useful combination */
1540: R = zeroseparableprimedecomposition(P,W,V);
1541: return R;
1542: }
1543: #endif
1544:
1545: def zerosepdec_main(P,W,V,Ord,Factors)
1546: {
1547: dp_ord(Ord);
1548: N = length(P);
1549: NFactors = length(Factors);
1550: PS = newvect(length(P),map(dp_ptod,P,W));
1551: for ( I = length(P)-1, Ind = [];I >= 0; I-- )
1552: Ind = cons(I,Ind);
1553: Y= setminus(V,W);
1554: NW = length(W);
1555: NY = length(Y);
1556: ElimV = append(W,Y);
1557: ElimOrd = [[Ord,NW],[0,NY]];
1558: ANS = [];
1559: for ( J = 0; J < NFactors; J++ ) {
1560: Factor = dp_dtop(nf_sfrat(Ind,dp_ptod(Factors[J],W),1,PS)[0],W);
1561: Q = dp_gr_f_main(cons(Factor,P),ElimV,Hom,ElimOrd);
1562: Q = dp_gr_f_main(Q,W,Hom,Ord);
1563: ANS = cons(Q,ANS);
1564: }
1565: return ANS;
1566: }
1567:
1568: def nextchoice(C)
1569: {
1570: N = size(C)[0];
1571: for ( I = N-1, Carry = 1; I >= 0; I-- ) {
1572: if ( C[I] ) {
1573: C[I] = !Carry;
1574: Carry = !C[I];
1575: } else {
1576: C[I] = Carry;
1577: Carry = 0;
1578: }
1579: }
1580: return Carry;
1581: }
1582:
1583: def separableclosure(CP,TargetVSet,VSet)
1584: {
1585: /* We compute a separable ideal <Q> from <P> by */
1586: /* computing inverse Frobenius map. */
1587:
1588: Ord=0;
1589:
1590: NVSet=length(TargetVSet);
1591: IndexVector=newvect(NVSet);
1592:
1593: /* First we check if <P> is already separable. */
1594: /* CHECK == 1 if <P> is already separable, */
1595: /* Otherwise, CHECK is the list of pairs of */
1596: /* the separable closure sc(m) of the minimal polynomial*/
1597: /* m and the exponent e, i.e. m=sc(m)(t^(q^e)). */
1598:
1599: CHECK=checkseparable(CP,TargetVSet,VSet);
1600:
1601: if ( CHECK == 1 )
1602: {
1603: print("This is already a separable ideal.", 2);
1604: return [CP[0],0];
1605: }
1606:
1607: print("This is not a separable ideal, so we make its separable closure.", 2);
1608:
1609: WSet=makecounterpart(TargetVSet);
1610: Char=setmod_ff()[0];
1611:
1612: NewP=CP[0];
1613: EXPVECTOR=newvect(NVSet);
1614: CHECKEXPVECTOR=newvect(NVSet);
1615:
1616: for (I=0;I<NVSet;I++)
1617: {
1618: EXPVECTOR[I]=CHECK[NVSet-I-1][0];
1619: CHECKEXPVECTOR[I]=deg(CHECK[NVSet-I-1][1],TargetVSet[I]);
1620:
1621: POW=Char^CHECK[NVSet-I-1][0];
1622:
1623: AddPoly=TargetVSet[I]^POW-WSet[I];
1624: /*SepClosure=subst(CHECK[I][0],TargetVSet[I],WSet[I]);*/
1625:
1626: NewP=cons(AddPoly,NewP);
1627: }
1628:
1629: NewOrder=[[Ord,NVSet],[Ord,NVSet]];
1630: XSet=append(TargetVSet,WSet);
1631: YSet=append(VSet,WSet);
1632:
1633: NewG=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
1634:
1635: NWSet=setminus(YSet,TargetVSet);
1636:
1637: Q=elimination(NewG,NWSet);
1638:
1639: NQ=length(Q);
1640:
1641: ANS=[];
1642: for (I=NQ-1;I>=0;I--)
1643: {
1644: Poly=Q[I];
1645: for (J=0;J<NVSet;J++)
1646: {
1647: Poly=subst(Poly,WSet[J],TargetVSet[J]);
1648: }
1649: ANS=cons(Poly,ANS);
1650: }
1651: return [ANS,EXPVECTOR,CHECKEXPVECTOR];
1652: }
1653:
1654: def convertdivisor(P,TargetVSet,VSet,ExVector)
1655: {
1656: /* We compute a corresponding ideal <Q> from <P> by */
1657: /* computing Frobenius map. */
1658:
1659: if (ExVector == 0 )
1660: {
1661: return P;
1662: }
1663:
1664: NVSet=length(TargetVSet);
1665: WSet=makecounterpart(TargetVSet);
1666: Char=setmod_ff()[0];
1667: Ord=0;
1668:
1669: NewP=P;
1670:
1671: for (I=0;I<NVSet;I++)
1672: {
1673: POW=Char^ExVector[I];
1674:
1675: AddPoly=TargetVSet[I]-WSet[I]^POW;
1676: /*SepClosure=subst(CHECK[I][0],TargetVSet[I],WSet[I]);*/
1677:
1678: NewP=cons(AddPoly,NewP);
1679: }
1680:
1681: NewOrder=[[Ord,NVSet],[Ord,NVSet]];
1682:
1683: XSet=append(TargetVSet,WSet);
1684: YSet=append(VSet,WSet);
1685:
1686: NewG=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
1687:
1688: NWSet=setminus(YSet,TargetVSet);
1689:
1690: Q=elimination(NewG,NWSet);
1691:
1692: NQ=length(Q);
1693:
1694: ANS=[];
1695: for (I=NQ-1;I>=0;I--)
1696: {
1697: Poly=Q[I];
1698: for (J=0;J<NVSet;J++)
1699: {
1700: Poly=subst(Poly,WSet[J],TargetVSet[J]);
1701: }
1702: ANS=cons(Poly,ANS);
1703: }
1704: return ANS;
1705: }
1706:
1707:
1708: def findgeneric(P,TargetVSet,VSet)
1709: {
1710: /* The number of trials is set as Trails. */
1711: NTargetVSet=length(TargetVSet);
1712: MNumber=0;
1713: /* Trials=100;*/
1714:
1715: if ( Trials == 0 )
1716: Trials = 2;
1717:
1718: Lineardimension=lineardimension(P,TargetVSet);
1719: NewVar=newt;
1720: Count=0;
1721: Ord=0;
1722:
1723: #if 0
1724: while(Count < Trials )
1725: {
1726: Candidate=0;
1727: MNumber=MNumber+1;
1728:
1729: MAGIC=makemagic(TargetVSet,MNumber);
1730:
1731: for (I=0;I<NTargetVSet;I++)
1732: {
1733: Candidate=Candidate+MAGIC[I]*TargetVSet[I];
1734: }
1735: MinPoly=minipoly_sf(P,TargetVSet,Ord,Candidate,NewVar);
1736: Deg=deg(MinPoly,NewVar);
1737: if ( Deg == Lineardimension )
1738: {
1739: return [Candidate,MinPoly,NewVar];
1740: }
1741: Count=Count+1;
1742: }
1743: #else
1744: YSet = setminus(VSet,TargetVSet);
1745: NYSet = length(YSet);
1746: Eval = newvect(NYSet);
1747: Q1 = field_order_ff()-1;
1748: HM = hmlist(P,TargetVSet,Ord);
1749: for ( Count = 0; Count < Trials; Count++ ) {
1750: Candidate = ptosfp(2)*TargetVSet[0];
1751: Candidate = 0;
1752: for ( I = 0; I < NTargetVSet; I++ )
1753: Candidate += random_ff()*TargetVSet[I];
1754: do {
1755: for ( I = 0; I < NYSet; I++ )
1756: Eval[I] = random()%Q1;
1757: } while ( !valid_modulus_sfrat(HM,YSet,Eval) );
1758: P0 = map(eval_sfrat,P,YSet,Eval);
1759: MinPoly0 = minipoly_sf(P0,TargetVSet,Ord,Candidate,NewVar);
1760: Deg = deg(MinPoly0,NewVar);
1761: if ( Deg == Lineardimension ) {
1762: MinPoly=minipoly_sf(P,TargetVSet,Ord,Candidate,NewVar);
1763: if ( deg(MinPoly,NewVar) != Deg )
1764: error("findgeneric : cannot happen");
1765: return [Candidate,MinPoly,NewVar];
1766: }
1767: }
1768: #endif
1769: print("Extend the ground field. ",2);
1770: error();
1771: }
1772:
1773: def makemagic(VSet,MNumber)
1774: {
1775: NVSet=length(VSet);
1776: MAGIC=[1];
1777: for (I=1;I<NVSet;I++)
1778: {
1779: U=ptosfp(MNumber^I);
1780: MAGIC=append(MAGIC,[U]);
1781: }
1782: return MAGIC;
1783: }
1784:
1785: #if 0
1786: def zerogenericprimedecomposition(CP,MinP,TargetVSet,VSet)
1787: {
1788: P=CP[0];
1789:
1790: FACTORS=sffctr(MinP);
1791: NFACTORS=length(FACTORS);
1792:
1793: if ( NFACTORS == 2 )
1794: {
1795: return [P];
1796: }
1797:
1798: ANS=[];
1799:
1800: for (I=1;I<NFACTORS;I++)
1801: {
1802: AddPoly=NFACTORS[I][0];
1803: NewP=cos(AddPoly,P);
1804: NewG=dp_gr_f_main(NewP,TargetVSet,Hom,Ord);
1805: ANS=cons(NewG,ANS);
1806: }
1807: return ANS;
1808: }
1809: #endif
1810:
1811:
1812: def lineardimension(P,VSet)
1813: {
1814: Ord=0;
1815:
1816: PP=dp_gr_f_main(P,VSet,Hom,Ord);
1817: dp_ord(Ord);
1818: Dimension=length(dp_mbase(map(dp_ptod,PP,VSet)));
1819: return Dimension;
1820: }
1821:
1822: def checkgeneric(CP,TargetVSet,VSet)
1823: {
1824: Dimension=lineardimension(CP[0],TargetVSet);
1825: NVSet=length(TargetVSet);
1826: Flag=0;
1827:
1828: for (I=0;I<NVSet;I++)
1829: {
1830: MinP=CP[1][I];
1831:
1832: if ( deg(MinP,TargetVSet[NVSet-I-1]) == Dimension )
1833: {
1834: return [1,MinP];
1835: }
1836: }
1837: return 0;
1838: }
1839:
1840: def checkseparable(CP,TargetVSet,VSet)
1841: {
1842: NVSet=length(TargetVSet);
1843: EXPVector=newvect(NVSet);
1844: Flag=1;
1845:
1846: for (I=0;I<NVSet;I++)
1847: {
1848: MinP=CP[1][I];
1849: CHECK=checkseparablepoly(MinP,TargetVSet[NVSet-I-1]);
1850: if ( CHECK[0] != 0 )
1851: {
1852: Flag = 0;
1853: }
1854: EXPVector[I]=CHECK;
1855: }
1856: if ( Flag == 0 )
1857: {
1858: return EXPVector;
1859: }
1860: return 1;
1861: }
1862:
1863: def checkseparablepoly(F,V)
1864: {
1865: P = characteristic_ff();
1866: E = 0;
1867: Vt = newt;
1868: while ( 1 ) {
1869: if ( diff(F,V) != 0 )
1870: return [E,F];
1871: T = srem(F,V^P-Vt,V);
1872: F = subst(T,Vt,V);
1873: E++;
1874: }
1875: }
1876:
1877: def extcont(P,V,Ord)
1878: {
1879: /* We assume P is a Groebner Basis w.r.t. Ord. */
1880: /* We compute a polynomial f such that */
1881: /* \sqrt{<P>}=\sqrt{<P>+<f>}\cap \sqrt{(<P>)^{ec}} */
1882: /* by B.W. Proposition 8.96. */
1883: /* Remark: The square free part of f is also OK. */
1884:
1885: dp_ord(Ord);
1886: HC = map(dp_hc,map(dp_ptod,P,V));
1887: LCM = car(HC);
1888: for ( T = cdr(HC); T != []; T = cdr(T) )
1889: LCM = lcm_sfrat(LCM,car(T));
1890: return LCM;
1891: }
1892:
1893: def first(A)
1894: {
1895: return A[0];
1896: }
1897:
1898: def set_union(A,B)
1899: {
1900: for ( T = A; T != []; T = cdr(T) )
1901: B = insert_element(car(T),B);
1902: return B;
1903: }
1904:
1905: def insert_element(E,A)
1906: {
1907: for ( T = A; T != []; T = cdr(T) )
1908: if ( E == car(T) )
1909: break;
1910: if ( T == [] )
1911: return cons(E,A);
1912: else
1913: return A;
1914: }
1915:
1916: def extcont_factor(P,V,Ord)
1917: {
1918: /* We assume P is a Groebner Basis w.r.t. Ord. */
1919: /* We compute a polynomial f such that */
1920: /* \sqrt{<P>}=\sqrt{<P>+<f>}\cap \sqrt{(<P>)^{ec}} */
1921: /* by B.W. Proposition 8.96. */
1922: /* Remark: The square free part of f is also OK. */
1923:
1924: dp_ord(Ord);
1925: HC = map(dp_hc,map(dp_ptod,P,V));
1926: F = [];
1927: for ( T = HC; T != []; T = cdr(T) ) {
1928: F1 = map(first,cdr(sffctr(car(T))));
1929: F = set_union(F1,F);
1930: }
1931: return F;
1932: }
1933:
1934: def contraction(P,V,W)
1935: {
1936: /* We assume P is a Groebner Basis w.r.t. Ord. */
1937: /* We compute the contraction <P>^{c} by */
1938: /* <P>^c=(<G>:f^\infty) */
1939: /* by B.W. Proposition 8.91. */
1940: /* This procedure is called by zeroprimedecomposition. */
1941: /* So, P is supposed to be a GB w.r.t. DRL. */
1942:
1943: Ord=0;
1944: YSet=setminus(W,V);
1945:
1946: Ord1 = [[Ord,length(V)],[0,length(YSet)]];
1947: GP1 = dp_gr_f_main(P,W,Hom,Ord1);
1948:
1949: Factor = extcont_factor(GP1,V,Ord);
1950: for ( F = 1, T = Factor; T != []; T = cdr(T) )
1951: F *= car(T);
1952: Vt = newt;
1953:
1954: G = dp_gr_f_main(cons(1-Vt*F,P),cons(Vt,W),0,[[0,1],[Ord,length(W)]]);
1955: R = [];
1956: for ( T = G; T != []; T = cdr(T) )
1957: if ( !member(Vt,vars(car(T))) )
1958: R = cons(car(T),R);
1959: return [R,F];
1960: }
1961:
1962: def checkgeneric2(LIST)
1963: {
1964: NList=size(LIST)[0];
1965:
1966: FLAG=0;
1967:
1968: for (I=0;I<NList;I++)
1969: {
1970: if ( LIST[I] > 1 )
1971: {
1972: FLAG=FLAG+1;
1973: }
1974: }
1975:
1976: if (FLAG < 2 )
1977: {
1978: return 1;
1979: }
1980: return 0;
1981: }
1982:
1983: #if 0
1984: def checkseparablepoly(P,V)
1985: {
1986: TestP=P;
1987: CHECK=diff(TestP,V);
1988: Count=0;
1989:
1990: while ( CHECK !=0 )
1991: {
1992: if ( deg(TestP,V) != 0 )
1993: {
1994: break;
1995: }
1996: TestP=pdivide(TestP);
1997: CHECK=diff(TestP,V);
1998: Count=Count+1;
1999: }
2000:
2001: return [TestP,Count];
2002: }
2003:
2004: def pdivide(F,V)
2005: {
2006: Char=setmod_ff()[0];
2007: TestP=P;
2008:
2009: Deg=ideg(TestP,V);
2010: DegP=idiv(Deg,Char);
2011:
2012: if ( irem(Deg,Char) != 0 )
2013: {
2014: error;
2015: }
2016: ANS=0;
2017:
2018: for (I=O;I<DegP;I++)
2019: {
2020: TempDeg=I*Char;
2021: ANS=ANS+coeff(TestP,V,TempDeg)*X^I;
2022: }
2023: return ANS;
2024: }
2025: #endif
2026:
2027:
2028: def convsf(PP,VSet,Ord,Flag)
2029: {
2030: /* Flag = 0 or 1 */
2031: /* If Flag = 1, we compute the intersection */
2032: /* of the Galois orbit. */
2033:
2034: CHECK=checkgaloisorbit(PP,VSet,Ord,Flag);
2035:
2036: NewPP=CHECK[0];
2037:
2038: ANS=[];
2039:
2040: NPP=length(NewPP);
2041:
2042: for (I=0;I<NPP;I++)
2043: {
2044: NewComp=convertsmallfield(CHECK[Flag][I],VSet,Ord);
2045: ANS=cons(NewComp,ANS);
2046: }
2047: return ANS;
2048: }
2049:
2050: def convertsmallfield(PP,VSet,Ord)
2051: {
2052: dp_ord(Ord);
2053: NVSet=length(VSet);
2054: Char=setmod_ff()[0];
2055: ExtDeg=deg(setmod_ff()[1],x);
2056:
2057: NewV=pg;
2058: MPP=map(monic_hc,PP,VSet);
2059: MPP=map(sfptopsfp,MPP,NewV);
2060:
2061: MinPoly=subst(setmod_ff()[1],x,NewV);
2062: XSet=cons(NewV,VSet);
2063:
2064: Ord1=[[0,1],[Ord,NVSet]];
2065:
2066: /* setmod_ff(Char,1);*/
2067:
2068: NewP=dp_gr_mod_main(cons(MinPoly,MPP),XSet,0,Char,Ord1);
2069:
2070: ANS=elimination(NewP,VSet);
2071:
2072: /* setmod_ff(Char,ExtDeg);*/
2073:
2074: return ANS;
2075: }
2076:
2077: def checkgaloisorbit(PP,VSet,Ord,Flag)
2078: {
2079: NPP=length(PP);
2080: TmpPP=PP;
2081: ExtDeg=deg(setmod_ff()[1],x);
2082:
2083: ANS=[];
2084: BNS=[];
2085:
2086: while (TmpPP != [] )
2087: {
2088: TestP=car(TmpPP);
2089: Dim=TestP[1];
2090: TmpPP=cdr(TmpPP);
2091: NewP=TestP[0];
2092:
2093: ANS=cons(TestP[0],ANS);
2094:
2095: for (J=1;J<ExtDeg;J++)
2096: {
2097: G1=map(sf_galois_action,TestP[0],J);
2098:
2099: if ( setminus(G1,TestP[0]) == [] )
2100: {
2101: break;
2102: }
2103: if (Flag == 1 )
2104: {
2105: NewP=ideal_intersection_sfrat(NewP,G1,VSet);
2106: }
2107: TmpPP=deletecomponent(TmpPP,[G1,Dim]);
2108: }
2109: BNS=cons(NewP,BNS);
2110:
2111: }
2112: return [ANS,BNS];
2113: }
2114:
2115: def deletecomponent(PP,G)
2116: {
2117: TmpPP=PP;
2118:
2119: while( TmpPP !=[] )
2120: {
2121: Test=car(TmpPP)[0];
2122:
2123: if ( setminus(Test,G[0]) == [] )
2124: {
2125: ANS=setminus(PP,[G]);
2126: return ANS;
2127: }
2128:
2129: TmpPP=cdr(TmpPP);
2130: }
2131: error();
2132: }
2133:
2134:
2135: def pfctr(F)
2136: {
2137: if ( type(F) == 1 )
2138: return [[F,1]];
2139: P = characteristic_ff();
2140: E = extdeg_ff();
2141: F = sfptop(F);
2142: check_coef(F,P);
2143: D = modfctr(F,P);
2144: setmod_ff(P,E);
2145: return map(mapsf_first,D);
2146: }
2147:
2148: def check_coef(F,P)
2149: {
2150: V = vars(F);
2151: D = dp_ptod(F,V);
2152: for ( T = D; T; T = dp_rest(T) )
2153: if ( dp_hc(T) >= P )
2154: error("invalid coef");
2155: }
2156:
2157: def mapsf_first(D)
2158: {
2159: return [ptosfp(D[0]),D[1]];
2160: }
2161:
2162: def partial_decomp(B,V)
2163: {
2164: T0 = time();
2165: if ( ParallelMinipoly ) {
2166: map(ox_cmo_rpc,ParallelMinipoly,"setmod_ff",characteristic_ff(),extdeg_ff());
2167: map(ox_pop_cmo,ParallelMinipoly);
2168: }
2169: B = map(ptosfp,B);
2170: B = dp_gr_f_main(B,V,0,0);
2171: R = partial_decomp0(B,V,length(V)-1);
2172: if ( PartialDecompByLex ) {
2173: R0 = [];
2174: for ( TR = R; TR != []; TR = cdr(TR) ) {
2175: T = car(TR);
2176: S = dp_gr_f_main(T[0],V,0,0);
2177: R0 = cons([S,T[1]],R0);
2178: }
2179: R = reverse(R0);
2180: }
2181: T_PD += time()[3]-T0[3];
2182: return R;
2183: }
2184:
2185: def setintersection(A,B)
2186: {
2187: for ( L = []; A != []; A = cdr(A) )
2188: if ( member(car(A),B) )
2189: L = cons(car(A),L);
2190: return L;
2191: }
2192:
2193: /* returns [[Plist,Mlist],...] */
2194:
2195: def partial_decomp0(B,V,I)
2196: {
2197: N = length(V);
2198: if ( I < 0 )
2199: return [[B,[]]];
2200: Ord = PartialDecompByLex ? [[0,I],[2,N-I]] : 0;
2201: #if 0
2202: if ( setminus(vars(B),V) == [] )
2203: B = fglm_sf_0dim(B,V,Ord0,Ord);
2204: else
2205: #endif
2206: B = dp_gr_f_main(B,V,0,Ord);
2207: if ( type(B[0]) == 1 )
2208: return [];
2209: if ( !zero_dim(B,V,Ord) )
2210: error("non zero-dimensional ideal");
2211: /* XXX */
2212: Vt = ttttt;
2213: VI = V[I];
2214: if ( PartialDecompByLex ) {
2215: for ( J = 0, W = V; J < I; J++ )
2216: W = cdr(W);
2217: VW = setminus(V,W);
2218: for ( Bw = [], T = B; T != []; T = cdr(T) )
2219: if ( setintersection(vars(car(T)),VW) == [] )
2220: Bw = cons(car(T),Bw);
2221: MI = minipoly_sf(Bw,W,2,VI,Vt);
2222: } else
2223: MI = minipoly_sf(B,V,0,VI,Vt);
2224: MIF = sffctr(MI);
2225: /* if MI is irreducible, then process the next variable */
2226: if ( length(MIF) == 2 && MIF[1][1] == 1 ) {
2227: L = partial_decomp0(B,V,I-1);
2228: R = [];
2229: for ( S = L; S != []; S = cdr(S) ) {
2230: L = car(S);
2231: R = cons([L[0],cons(subst(MI,Vt,VI),L[1])],R);
2232: }
2233: return R;
2234: }
2235:
2236: /* for nf computation */
2237: Ord1 = PartialDecompByLex ? 2 : [[0,1],[0,N]];
2238: Len = length(B);
2239: PS = newvect(Len+1);
2240: dp_ord(Ord1);
2241: XV = cons(Vt,V);
2242: for ( J = 0, T = cons(Vt-VI,B); T != []; T = cdr(T), J++ )
2243: PS[J] = dp_ptod(car(T),XV);
2244: for ( J = 0, GI = []; J <= Len; J++ )
2245: GI = cons(J,GI);
2246: if ( PartialDecompByLex )
2247: GI = reverse(GI);
2248:
2249: R = [];
2250: for ( T = MIF; T != []; T = cdr(T) ) {
2251: Mt = car(car(T));
2252: if ( type(Mt) == 1 )
2253: continue;
2254: dp_ord(Ord1);
2255: NfMt = dp_dtop(nf_sfrat(GI,dp_ptod(Mt,XV),1,PS)[0],XV);
2256: if ( NfMt )
2257: B1 = dp_gr_f_main(cons(NfMt,B),V,0,Ord);
2258: else
2259: B1 = B;
2260: Rt = partial_decomp0(B1,V,I-1);
2261: for ( S = Rt; S != []; S = cdr(S) ) {
2262: L = car(S);
2263: R = cons([L[0],cons(subst(Mt,Vt,VI),L[1])],R);
2264: }
2265: }
2266: return R;
2267: }
2268:
2269:
2270: /* G is a gb wrt (V,O) */
2271:
2272: def minipoly_sf(G,V,O,F,V0)
2273: {
2274: T0 = time();
2275: Vc = cons(V0,setminus(vars(G),V));
2276: if ( ParallelMinipoly ) {
2277: Proc0 = ParallelMinipoly[0];
2278: Proc1 = ParallelMinipoly[1];
2279: if ( length(Vc) == 1 ) {
2280: ox_rpc(Proc0,"minipoly_sf_by_buchberger",G,V,O,F,V0,1);
2281: ox_rpc(Proc1,"minipoly_sf_0dim",G,V,O,F,V0,1);
2282: map(ox_get,ParallelMinipoly);
2283: /* 258=SM_popSerializedLocalObject */
2284: map(ox_push_cmd,ParallelMinipoly,258);
2285: F = ox_select(ParallelMinipoly);
2286: MP = ox_get(F[0]);
2287: if ( F[0] == Proc0 ) {
2288: if ( length(F) == 1 )
2289: B_Win++;
2290: else
2291: ox_get(Proc1);
2292: } else {
2293: if ( length(F) == 1 )
2294: D_Win++;
2295: else
2296: ox_get(Proc0);
2297: }
2298: ox_reset(Proc0);
2299: ox_reset(Proc1);
2300: } else if ( length(Vc) == 2 ) {
2301: ox_rpc(Proc0,"minipoly_sf_by_buchberger",G,V,O,F,V0,1);
2302: ox_rpc(Proc1,"minipoly_sfrat",G,V,O,F,V0,1);
2303: map(ox_get,ParallelMinipoly);
2304: /* 258=SM_popSerializedLocalObject */
2305: map(ox_push_cmd,ParallelMinipoly,258);
2306: F = ox_select(ParallelMinipoly);
2307: MP = ox_get(F[0]);
2308: if ( F[0] == Proc0 ) {
2309: if ( length(F) == 1 )
2310: B_Win++;
2311: else
2312: ox_get(Proc1);
2313: } else {
2314: if ( length(F) == 1 )
2315: D_Win++;
2316: else
2317: ox_get(Proc0);
2318: }
2319: ox_reset(Proc0);
2320: ox_reset(Proc1);
2321: } else
2322: MP = minipoly_sf_by_buchberger(G,V,O,F,V0,0);
2323: } else if ( BuchbergerMinipoly )
2324: MP = minipoly_sf_by_buchberger(G,V,O,F,V0,0);
2325: else {
2326: if ( length(Vc) == 1 )
2327: MP = minipoly_sf_0dim(G,V,O,F,V0,0);
2328: else if ( length(Vc) == 2 )
2329: MP = minipoly_sfrat(G,V,O,F,V0,0);
2330: else
2331: MP = minipoly_sf_by_buchberger(G,V,O,F,V0,0);
2332: }
2333: T_MP += time()[3]-T0[3];
2334: return MP;
2335: }
2336:
2337: def minipoly_sf_by_buchberger(G,V,O,F,V0,Server)
2338: {
2339: if ( Server )
2340: ox_sync(0);
2341: Vc = cons(V0,setminus(vars(G),V));
2342: Gf = cons(ptosfp(V0-F),G);
2343: Vf = append(V,Vc);
2344: Gelim = dp_gr_f_main(Gf,Vf,1,[[0,length(V)],[0,length(Vc)]]);
2345: for ( Gc = [], T = Gelim; T != []; T = cdr(T) ) {
2346: Vt = setminus(vars(car(T)),Vc);
2347: if ( Vt == [] )
2348: Gc = cons(car(T),Gc);
2349: }
2350: Gt = dp_gr_f_main(Gc,Vc,1,[[0,1],[0,length(Vc)-1]]);
2351: Pmin = car(Gt); Dmin = deg(Pmin,V0);
2352: for ( T = cdr(Gt); T != []; T = cdr(T) ) {
2353: Dt = deg(car(T),V0);
2354: if ( Dt < Dmin ) {
2355: Pmin = car(T); Dmin = Dt;
2356: }
2357: }
2358: Cont = sfcont(Pmin,V0);
2359: return sdiv(Pmin,Cont);
2360: }
2361:
2362: def minipoly_sf_0dim(G,V,O,F,V0,Server)
2363: {
2364: if ( Server )
2365: ox_sync(0);
2366: N = length(V);
2367: Len = length(G);
2368: dp_ord(O);
2369: PS = newvect(Len);
2370: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
2371: PS[I] = dp_ptod(car(T),V);
2372: for ( I = Len-1, HL = []; I >= 0; I-- )
2373: HL = cons(dp_ht(PS[I]),HL);
2374: for ( I = Len - 1, GI = []; I >= 0; I-- )
2375: GI = cons(I,GI);
2376: MB = dp_mbase(HL); DIM = length(MB); UT = newvect(DIM);
2377: U = dp_ptod(ptosfp(F),V);
2378: U = dp_nf_f(GI,U,PS,1);
2379: for ( I = 0; I < DIM; I++ )
2380: UT[I] = [MB[I],dp_nf_f(GI,U*MB[I],PS,1)];
2381:
2382: T = dp_ptod(ptosfp(1),[V0]);
2383: TT = dp_ptod(ptosfp(1),V);
2384: G = H = [[TT,T]];
2385:
2386: for ( I = 1; ; I++ ) {
2387: if ( dp_gr_print() )
2388: print(".",2);
2389: T = dp_ptod(ptosfp(V0^I),[V0]);
2390: TT = dp_nf_tab_f(H[0][0],UT);
2391: H = cons([TT,T],H);
2392: L = dp_lnf_f([TT,T],G);
2393: if ( !L[0] ) {
2394: if ( dp_gr_print() )
2395: print("");
2396: return dp_dtop(L[1],[V0]); /* XXX */
2397: } else
2398: G = insert(G,L);
2399: }
2400: }
2401:
2402: #if 1
2403: def minipoly_sf_rat(G,V,F,V0)
2404: {
2405: Vc = setminus(vars(G),V);
2406: Gf = cons(V0-F,G);
2407: Vf = append(V,[V0]);
2408: G3 = dp_gr_f_main(map(ptosfp,Gf),Vf,0,3);
2409: for ( T = G3; T != []; T = cdr(T) ) {
2410: Vt = setminus(vars(car(T)),Vc);
2411: if ( Vt == [V0] )
2412: return car(T);
2413: }
2414: }
2415:
2416: def minipoly_mod(G,V,Mod,F,V0)
2417: {
2418: Vc = setminus(vars(G),V);
2419: Gf = cons(V0-F,G);
2420: Vf = append(V,[V0]);
2421: G3 = dp_gr_mod_main(Gf,Vf,1,Mod,3);
2422: for ( T = G3; T != []; T = cdr(T) ) {
2423: Vt = setminus(vars(car(T)),Vc);
2424: if ( Vt == [V0] )
2425: return car(T);
2426: }
2427: }
2428: #endif
2429:
2430: /* find N/D s.t. F = N/D mod V^(2M), deg(N),deg(D)<M */
2431:
2432: def pade(F,M)
2433: {
2434: V = var(F);
2435: A = newvect(M);
2436: D = 0;
2437: for ( I = 0; I < M; I++ ) {
2438: A[I] = strtov("b"+rtostr(I));
2439: D += A[I]*V^I;
2440: }
2441: T = F*D;
2442: M2 = M*2;
2443: R = [];
2444: for ( I = M; I < M2; I++ )
2445: R = cons(coef(T,I,V),R);
2446: }
2447:
2448: def minipoly_sfrat(G0,V,O,P,V0,Server)
2449: {
2450: if ( Server )
2451: ox_sync(0);
2452: if ( !zero_dim(hmlist(G0,V,O),V,O) )
2453: error("minipoly_sfrat : ideal is not zero-dimensional!");
2454:
2455: G1 = cons(V0-P,G0);
2456: if ( type(O) <= 1 )
2457: O1 = [[0,1],[O,length(V)]];
2458: else
2459: O1 = cons([0,1],O);
2460: V1 = cons(V0,V);
2461: W = append(V,[V0]);
2462: Vc = setminus(vars(G0),V);
2463:
2464: N = length(V1);
2465: dp_ord(O1);
2466: HM = hmlist(G1,V1,O1);
2467: MB = dp_mbase(map(dp_ptod,HM,V1));
2468: dp_ord(O);
2469:
2470: Nc = length(Vc);
2471: Eval = newvect(Nc);
2472:
2473: /* heristic lower bound of deg(MP) */
2474: Q1 = field_order_ff()-1;
2475: do {
2476: for ( I = 0; I < Nc; I++ )
2477: Eval[I] = random()%Q1;
2478: } while ( !valid_modulus_sfrat(HM,Vc,Eval) );
2479: G0E = map(eval_sfrat,G0,Vc,Eval);
2480: P0 = eval_sfrat(P,Vc,Eval);
2481: MP = minipoly_sf(G0E,V,O,P0,V0);
2482: DMP = deg(MP,V0);
2483:
2484: for ( I = 0; I < Nc; I++ )
2485: Eval[I] = 0;
2486: while ( 1 ) {
2487: if ( valid_modulus_sfrat(HM,Vc,Eval) ) {
2488: G0E = map(eval_sfrat,G0,Vc,Eval);
2489: P0 = eval_sfrat(P,Vc,Eval);
2490: MP = minipoly_sf(G0E,V,O,P0,V0);
2491: D = deg(MP,V0);
2492: if ( D >= DMP ) {
2493: DMP = D;
2494: for ( TL = [], MPT = 0, J = 0; J <= D; J++ ) {
2495: TL = cons(V0^J,TL);
2496: MPT += car(TL);
2497: }
2498: NF = gennf_sfrat(G1,TL,V1,O1,V0,1)[0];
2499: R = tolex_sfrat_main(V1,O1,NF,[MPT],Vc,Eval,MB);
2500: if ( R )
2501: return sdiv(R[0],sfcont(R[0],V0));
2502: }
2503: }
2504: next_eval_sfrat(Eval);
2505: }
2506: }
2507:
2508: def next_eval_sfrat(Eval)
2509: {
2510: N = size(Eval)[0];
2511: for ( I = N-1; I >= 0; I-- )
2512: if ( Eval[I] ) break;
2513: if ( I < 0 ) Eval[N-1] = 1;
2514: else if ( I == 0 ) {
2515: T = Eval[0]; Eval[0] = 0; Eval[N-1] = T+1;
2516: } else {
2517: Eval[I-1]++; T = Eval[I];
2518: for ( J = I; J < N-1; J++ )
2519: Eval[J] = 0;
2520: Eval[N-1] = T-1;
2521: }
2522: }
2523:
2524: def eval_sfrat(F,V,Eval)
2525: {
2526: for ( I = 0; V != []; V = cdr(V), I++ )
2527: F = subst(F,car(V),ptosfp(Eval[I]));
2528: return F;
2529: }
2530:
2531: def valid_modulus_sfrat(HL,Vc,Eval) {
2532: for ( T = HL; T != []; T = cdr(T) ) {
2533: C = car(T);
2534: for ( S = Vc, I = 0; S != []; S = cdr(S), I++ )
2535: C = subst(C,car(S),ptosfp(Eval[I]));
2536: if ( !C )
2537: break;
2538: }
2539: return T == [] ? 1 : 0;
2540: }
2541:
2542: def gennf_sfrat(G,TL,V,O,V0,FLAG)
2543: {
2544: N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
2545: for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
2546: PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
2547: }
2548: for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
2549: DTL = cons(dp_ptod(car(TL),V),DTL);
2550: for ( I = Len - 1, GI = []; I >= 0; I-- )
2551: GI = cons(I,GI);
2552: T = car(DTL); DTL = cdr(DTL);
2553: H = [nf_sfrat(GI,T,T,PS)];
2554:
2555: USE_TAB = (FLAG != 0);
2556: if ( USE_TAB ) {
2557: T0 = time()[0];
2558: MB = dp_mbase(HL); DIM = length(MB);
2559: U = dp_ptod(V0,V);
2560: UTAB = newvect(DIM);
2561: for ( I = 0; I < DIM; I++ ) {
2562: UTAB[I] = [MB[I],remove_cont_sfrat(nf_sfrat(GI,U*MB[I],1,PS))];
2563: if ( dp_gr_print() )
2564: print(".",2);
2565: }
2566: if ( dp_gr_print() )
2567: print("");
2568: TTAB = time()[0]-T0;
2569: }
2570:
2571: T0 = time()[0];
2572: for ( LCM = 1; DTL != []; ) {
2573: if ( dp_gr_print() )
2574: print(".",2);
2575: T = car(DTL); DTL = cdr(DTL);
2576: if ( L = search_redble(T,H) ) {
2577: DD = dp_subd(T,L[1]);
2578: if ( USE_TAB && (DD == U) ) {
2579: NF = nf_tab_sfrat(L[0],UTAB);
2580: NF = [NF[0],dp_hc(L[1])*NF[1]*T];
2581: } else
2582: NF = nf_sfrat(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
2583: } else
2584: NF = nf_sfrat(GI,T,T,PS);
2585: NF = remove_cont_sfrat(NF);
2586: H = cons(NF,H);
2587: LCM = lcm_sfrat(LCM,dp_hc(NF[1]));
2588: }
2589: TNF = time()[0]-T0;
2590: if ( dp_gr_print() )
2591: print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
2592: return [[map(adj_dn_sfrat,H,LCM),LCM],PS,GI];
2593: }
2594:
2595: def lcm_sfrat(A,B)
2596: {
2597: G = sfgcd(A,B);
2598: return sdiv(A,G)*B;
2599: }
2600:
2601: #define REDCONT(f) ((f)/sfcont(f))
2602:
2603: def remove_cont_sfrat(L)
2604: {
2605: if ( type(L[1]) <= 2 ) {
2606: T = remove_cont_sfrat([L[0],L[1]*<<0>>]);
2607: return [T[0],dp_hc(T[1])];
2608: } else if ( !L[0] )
2609: return [0,REDCONT(L[1])];
2610: else if ( !L[1] )
2611: return [REDCONT(L[0]),0];
2612: else {
2613: A0 = REDCONT(L[0]); A1 = REDCONT(L[1]);
2614: C0 = sdiv(dp_hc(L[0]),dp_hc(A0)); C1 = sdiv(dp_hc(L[1]),dp_hc(A1));
2615: GCD = sfgcd(C0,C1); M0 = sdiv(C0,GCD); M1 = sdiv(C1,GCD);
2616: return [M0*A0,M1*A1];
2617: }
2618: }
2619:
2620: def nf_sfrat(B,G,M,PS)
2621: {
2622: for ( D = 0; G; ) {
2623: for ( U = 0, L = B; L != []; L = cdr(L) ) {
2624: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
2625: GCD = sfgcd(dp_hc(G),dp_hc(R));
2626: CG = sdiv(dp_hc(R),GCD); CR = sdiv(dp_hc(G),GCD);
2627: U = CG*G-dp_subd(G,R)*CR*R;
2628: if ( !U )
2629: return [D,M];
2630: D *= CG; M *= CG;
2631: break;
2632: }
2633: }
2634: if ( U )
2635: G = U;
2636: else {
2637: D += dp_hm(G); G = dp_rest(G);
2638: }
2639: }
2640: return [D,M];
2641: }
2642:
2643: def nf_tab_sfrat(F,TAB)
2644: {
2645: for ( NM = 0, DN = 1, I = 0; F; F = dp_rest(F) ) {
2646: T = dp_ht(F);
2647: for ( ; TAB[I][0] != T; I++);
2648: NF = TAB[I][1]; N = NF[0]; D = NF[1];
2649: G = sfgcd(DN,D); DN1 = sdiv(DN,G); D1 = sdiv(D,G);
2650: NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
2651: }
2652: return [NM,DN];
2653: }
2654:
2655: def adj_dn_sfrat(P,D)
2656: {
2657: return [(sdiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
2658: }
2659:
2660: def tolex_sfrat_main(V,O,NF,GM,Vc,Eval,MB)
2661: {
2662: if ( MB ) {
2663: PosDim = 0;
2664: DIM = length(MB);
2665: DV = newvect(DIM);
2666: } else
2667: PosDim = 1;
2668: for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
2669: S = p_terms(car(T),V,2);
2670: if ( PosDim ) {
2671: MB = gather_nf_terms(S,NF,V,O);
2672: DV = newvect(length(MB));
2673: }
2674: dp_ord(O); RHS = termstomat_sfrat(NF,map(dp_ptod,cdr(S),V),MB,Vc,Eval);
2675: dp_ord(O); NHT = nf_tab_gsl_sfrat(dp_ptod(LCM*car(S),V),NF);
2676: dptov(NHT[0],DV,MB);
2677: dp_ord(O); B = hen_ttob_gsl_sfrat([DV,NHT[1]],RHS,cdr(S),Vc,Eval);
2678: if ( !B )
2679: return 0;
2680: Len = length(S);
2681: LCM *= B[1];
2682: for ( U = LCM*car(S), I = 1; I < Len; I++ )
2683: U += B[0][I-1]*S[I];
2684: SL = cons(U,SL);
2685: if ( dp_gr_print() )
2686: print(["DN",B[1]]);
2687: }
2688: return SL;
2689: }
2690:
2691: def termstomat_sfrat(NF,TERMS,MB,Vc,Eval)
2692: {
2693: DN = NF[1];
2694: NF = NF[0];
2695: N = length(MB);
2696: M = length(TERMS);
2697: MAT = newmat(N,M);
2698: W = newvect(N);
2699: Len = length(NF);
2700: for ( I = 0; I < M; I++ ) {
2701: T = TERMS[I];
2702: for ( K = 0; K < Len; K++ )
2703: if ( T == NF[K][1] )
2704: break;
2705: dptov(NF[K][0],W,MB);
2706: for ( J = 0; J < N; J++ )
2707: MAT[J][I] = W[J];
2708: }
2709: return [henleq_prep_sfrat(MAT,Vc,Eval),DN];
2710: }
2711:
2712: def henleq_prep_sfrat(A,Vc,Eval)
2713: {
2714: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
2715: L = geninv_sf_swap(map(eval_sfrat,A,Vc,Eval)); INV = L[0]; INDEX = L[1];
2716: AA = newmat(COL,COL);
2717: for ( I = 0; I < COL; I++ )
2718: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL; J++ )
2719: T[J] = S[J];
2720: if ( COL != ROW ) {
2721: RESTA = newmat(ROW-COL,COL);
2722: for ( ; I < ROW; I++ )
2723: for ( J = 0, T = RESTA[I-COL], S = A[INDEX[I]]; J < COL; J++ )
2724: T[J] = S[J];
2725: } else
2726: RESTA = 0;
2727: return [[A,AA,RESTA],L];
2728: }
2729:
2730: def hen_ttob_gsl_sfrat(LHS,RHS,TERMS,Vc,Eval)
2731: {
2732: LDN = LHS[1]; RDN = RHS[1]; LCM = lcm_sfrat(LDN,RDN);
2733: L1 = sdiv(LCM,LDN); R1 = sdiv(LCM,RDN);
2734: T0 = time()[0];
2735: S = henleq_gsl_sfrat(RHS[0],LHS[0]*L1,Vc,Eval);
2736: if ( dp_gr_print() )
2737: print(["henleq_gsl_sfrat",time()[0]-T0]);
2738: if ( !S )
2739: return 0;
2740: N = length(TERMS);
2741: return [S[0],S[1]*R1];
2742: }
2743:
2744: def nf_tab_gsl_sfrat(A,NF)
2745: {
2746: DN = NF[1];
2747: NF = NF[0];
2748: TLen = length(NF);
2749: for ( R = 0; A; A = dp_rest(A) ) {
2750: HM = dp_hm(A); C = dp_hc(HM); T = dp_ht(HM);
2751: for ( I = 0; I < TLen; I++ )
2752: if ( NF[I][1] == T )
2753: break;
2754: R += C*NF[I][0];
2755: }
2756: return remove_cont_sfrat([R,DN]);
2757: }
2758:
2759: def henleq_gsl_sfrat(L,B,Vc,Eval)
2760: {
2761: Nc = length(Vc);
2762: if ( Nc > 1 )
2763: return henleq_gsl_sfrat_higher(L,B,Vc,Eval);
2764:
2765: V0 = Vc[0]; E0 = ptosfp(Eval[0]);
2766:
2767: AL = L[0]; INVL = L[1];
2768: A = AL[0]; AA = AL[1]; RESTA = AL[2];
2769: INV = INVL[0]; INDEX = INVL[1];
2770:
2771: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
2772: BB = newvect(COL);
2773: for ( I = 0; I < COL; I++ )
2774: BB[I] = B[INDEX[I]];
2775: if ( COL != ROW ) {
2776: RESTB = newvect(ROW-COL);
2777: for ( ; I < ROW; I++ )
2778: RESTB[I-COL] = B[INDEX[I]];
2779: } else
2780: RESTB = 0;
2781:
2782: COUNT = 2;
2783: if ( !COUNT )
2784: COUNT = 2;
2785: X = newvect(size(AA)[0]);
2786: for ( I = 0, C = BB, CCC = 0, ITOR_FAIL = -1; ; I++ ) {
2787: if ( zerovector(C) ) {
2788: X = map(subst,X,V0,V0-E0);
2789: if ( zerovector(RESTA*X+RESTB) ) {
2790: if ( dp_gr_print() ) print("end",0);
2791: return [X,ptosfp(1)];
2792: } else
2793: return 0;
2794: } else if ( COUNT == CCC ) {
2795: CCC = 0;
2796: ND = polyvtoratv(X,V0,I);
2797: if ( ND ) {
2798: F = map(subst,ND[0],V0,V0-E0);
2799: LCM = subst(ND[1],V0,V0-E0);
2800: T = AA*F+LCM*BB;
2801: if ( zerovector(T) ) {
2802: T = RESTA*F+LCM*RESTB;
2803: if ( zerovector(T) ) {
2804: if ( dp_gr_print() ) print("end",0);
2805: return [F,LCM];
2806: } else
2807: return 0;
2808: }
2809: } else {
2810: }
2811: } else {
2812: if ( dp_gr_print() ) print(".",2);
2813: CCC++;
2814: }
2815: XT = -INV*map(subst,C,V0,E0);
2816: X += XT*V0^I;
2817: C += AA*XT;
2818: C = map(sdiv,C,V0-E0);
2819: }
2820: }
2821:
2822: def henleq_gsl_sfrat_higher(L,B,Vc,Eval)
2823: {
2824: Nc = length(Vc);
2825: E = map(ptosfp,Eval);
2826:
2827: AL = L[0]; INVL = L[1];
2828: A = AL[0]; AA0 = AL[1]; RESTA = AL[2];
2829: INV = INVL[0]; INDEX = INVL[1];
2830:
2831: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
2832:
2833: AA = map(mshift,AA0,Vc,E,1);
2834: BB = newvect(COL);
2835: for ( I = 0; I < COL; I++ )
2836: BB[I] = mshift(B[INDEX[I]],Vc,E,1);
2837:
2838: if ( COL != ROW ) {
2839: RESTB = newvect(ROW-COL);
2840: for ( ; I < ROW; I++ )
2841: RESTB[I-COL] = B[INDEX[I]];
2842: } else
2843: RESTB = 0;
2844:
2845: COUNT = 2;
2846: if ( !COUNT )
2847: COUNT = 2;
2848: X = newvect(size(AA)[0]);
2849: for ( I = 0, R = BB, CCC = 0, ITOR_FAIL = -1; ; I++ ) {
2850: if ( zerovector(R) ) {
2851: X = map(mshift,X,Vc,E,-1);
2852: if ( zerovector(RESTA*X+RESTB) ) {
2853: if ( dp_gr_print() ) print("end",0);
2854: return [X,ptosfp(1)];
2855: } else
2856: return 0;
2857: } else if ( COUNT == CCC ) {
2858: CCC = 0;
2859: ND = polyvtoratv_higher(X,Vc,I);
2860: if ( ND ) {
2861: F = map(mshift,ND[0],Vc,E,-1);
2862: LCM = mshift(ND[1],Vc,E,-1);
2863: T = AA*F+LCM*BB;
2864: if ( zerovector(T) ) {
2865: T = RESTA*F+LCM*RESTB;
2866: if ( zerovector(T) ) {
2867: if ( dp_gr_print() ) print("end",0);
2868: return [F,LCM];
2869: } else
2870: return 0;
2871: }
2872: } else {
2873: }
2874: } else {
2875: if ( dp_gr_print() ) print(".",2);
2876: CCC++;
2877: }
2878: RK = map(mtrunc,R,I+1);
2879: XT = -INV*RK;
2880: X += XT;
2881: R += AA*XT;
2882: }
2883: }
2884:
2885: /* V -> V+Sgn*E */
2886:
2887: def mshift(F,V,E,Sgn)
2888: {
2889: N = length(V);
2890: for ( I = 0; I < N; I++ )
2891: F = subst(F,V[I],V[I]+Sgn*E[I]);
2892: return F;
2893: }
2894:
2895: /* truncate terms whose degree are higher than or equal to D */
2896:
2897: def mtrunc(F,D)
2898: {
2899: if ( type(F) <= 1 )
2900: return F;
2901: else {
2902: R = 0;
2903: V = var(F);
2904: while ( F ) {
2905: K = deg(F,V);
2906: C = coef(F,K,V);
2907: if ( K < D )
2908: R += mtrunc(C,D-K)*V^K;
2909: F -= C*V^K;
2910: }
2911: return R;
2912: }
2913: }
2914:
2915: /* for 1-dim case */
2916:
2917: def polyvtoratv(Vect,V,K)
2918: {
2919: N = size(Vect)[0];
2920: R = newvect(N);
2921: K2 = idiv(K,2);
2922: Mat = newmat(K2,K2);
2923: for ( I = 0; I < N; I++ ) {
2924: T = polytorat(Vect[I],V,Mat,K2);
2925: if ( T )
2926: R[I] = T;
2927: else
2928: return 0;
2929: }
2930: LCM = R[0][1];
2931: for ( I = 1; I < N; I++ )
2932: LCM = lcm_sfrat(LCM,R[I][1]);
2933: for ( I = 0; I < N; I++ )
2934: R[I] = R[I][0]*sdiv(LCM,R[I][1]);
2935: return [R,LCM];
2936: }
2937:
2938: /* for higher dim case */
2939:
2940: def polyvtoratv_higher(Vect,Vc,K)
2941: {
2942: N = size(Vect)[0];
2943: R = newvect(N);
2944: K2 = idiv(K,2);
2945: for ( I = 0; I < N; I++ ) {
2946: T = polytorat_higher(Vect[I],Vc,K2);
2947: if ( T )
2948: R[I] = T;
2949: else
2950: return 0;
2951: }
2952: LCM = R[0][1];
2953: for ( I = 1; I < N; I++ )
2954: LCM = lcm_sfrat(LCM,R[I][1]);
2955: for ( I = 0; I < N; I++ )
2956: R[I] = R[I][0]*sdiv(LCM,R[I][1]);
2957: return [R,LCM];
2958: }
2959:
2960: /* find F = N/D mod V^(2K), deg(N), deg(D) < K */
2961: def polytorat_gcd(F,V,K)
2962: {
2963: if ( deg(F,V) < K )
2964: return [F,ptosfp(1)];
2965: F1 = Mod^(K*2); F2 = F;
2966: B1 = 0; B2 = 1;
2967: while ( 1 ) {
2968: Q = sdiv(F1,F2);
2969: F3 = F1-F2*Q;
2970: B3 = B1-B2*Q;
2971: if ( deg(F3,V) < K ) {
2972: if ( deg(B3,V) < K )
2973: return [F3,B3];
2974: else
2975: return 0;
2976: }
2977: F1 = F2; F2 = F3;
2978: B1 = B2; B2 = B3;
2979: }
2980: }
2981:
2982: /*
2983: * for 1-dim case
2984: *
2985: * solve
2986: *
2987: * [fk ... f1][ a0 ]
2988: * [f(k+1) ... f2][ a1 ] = 0
2989: * [ ... ][ ... ]
2990: * [f(2k-1) ... fk][a(k-1)]
2991: */
2992:
2993: def polytorat(F,V,Mat,K)
2994: {
2995: if ( deg(F,V) < K )
2996: return [F,ptosfp(1)];
2997: for ( I = 0; I < K; I++ )
2998: for ( J = 0; J < K; J++ )
2999: Mat[I][J] = coef(F,I+K-J);
3000: S = nullspace_ff(Mat);
3001: MT = S[0]; IND = S[1];
3002: for ( I = 0; I < K; I++ )
3003: if ( IND[I] )
3004: break;
3005: if ( I == K )
3006: return 0;
3007: D = null_to_poly_ff(MT,IND,V)[0];
3008: N = trunc(F*D,K-1);
3009: return [N,D];
3010: }
3011:
3012: /* find N,D s.t. tdeg(N), tdeg(D) < K, F = N/D mod <V0,...,VM>^(2K) */
3013:
3014: def polytorat_higher(F,V,K)
3015: {
3016: if ( K < 2 ) return 0;
3017: if ( homogeneous_deg(F) < K )
3018: return [F,ptosfp(1)];
3019: D = create_icpoly(V,K);
3020: C = extract_coef(D*F,V,K,2*K);
3021: Vc = vars(C);
3022: G = dp_gr_f_main(C,Vc,0,2);
3023: PS = newvect(length(G),map(dp_ptod,G,Vc));
3024: for ( I = length(G)-1, Ind = [];I >= 0; I-- )
3025: Ind = cons(I,Ind);
3026: D = dp_dtop(nf_sfrat(Ind,dp_ptod(D,Vc),1,PS)[0],Vc);
3027: if ( !D )
3028: return 0;
3029: Vp = setminus(vars(D),V);
3030: D = subst(D,car(Vp),1);
3031: for ( T = cdr(Vp); T != []; T = cdr(T) )
3032: D = subst(D,car(T),0);
3033: N = mtrunc(F*D,K);
3034: return [N,D];
3035: }
3036:
3037: def create_icpoly(V,K)
3038: {
3039: if ( V == [] )
3040: return uc();
3041: R = 0;
3042: for ( I = 0; I < K; I++ )
3043: R += create_icpoly(cdr(V),K-I)*V[0]^I;
3044: return R;
3045: }
3046:
3047: /* extract terms whose degrees are in [From,To) */
3048:
3049: def extract_coef(F,V,From,To)
3050: {
3051: if ( V == [] ) {
3052: if ( From == 0 )
3053: return [F];
3054: else
3055: return [];
3056: }
3057: R = [];
3058: V0 = V[0];
3059: for ( I = 0; I < To; I++ ) {
3060: C = coef(F,I,V0);
3061: if ( C ) {
3062: C = extract_coef(C,cdr(V),From>=I?From-I:0,To-I);
3063: R = append(C,R);
3064: }
3065: }
3066: return R;
3067: }
3068:
3069: def saturation_sfrat(F,G,V,Ord)
3070: {
3071: Vt = ttttt;
3072: G0 = dp_gr_f_main(cons(Vt*F-1,G),cons(Vt,V),0,[[0,1],[Ord,length(V)]]);
3073: return elimination(G0,V);
3074: }
3075:
3076: def ideal_list_intersection_sfrat(L,V)
3077: {
3078: R = car(L);
3079: for ( TL = cdr(L); TL != []; TL = cdr(TL) )
3080: R = ideal_intersection_sfrat(R,car(TL),V);
3081: return R;
3082: }
3083:
3084: def ideal_intersection_sfrat(A,B,V)
3085: {
3086: T0 = time();
3087: Vt = ttttt;
3088: C = [];
3089: for ( T = A; T != []; T = cdr(T) )
3090: C = cons(car(T)*Vt,C);
3091: for ( T = B; T != []; T = cdr(T) )
3092: C = cons(car(T)*(1-Vt),C);
3093: Ord = [[0,1],[0,length(V)]];
3094: #if 0
3095: G0 = dp_gr_f_main(C,cons(Vt,V),0,0);
3096: G = dp_gr_f_main(G0,cons(Vt,V),0,Ord);
3097: #else
3098: G = dp_gr_f_main(C,cons(Vt,V),0,Ord);
3099: #endif
3100: T_INT += time()[3]-T0[3];
3101: return elimination(G,V);
3102: }
3103:
3104: /* Shimoyama's gr_fctr */
3105:
3106: def idealsqfr_sf(G)
3107: {
3108: for(LL=[],I=length(G)-1;I>=0;I--) {
3109: for(A=1,L=sfsqfr(G[I]),J=1;J<length(L);J++)
3110: A*=L[J][0];
3111: LL=cons(A,LL);
3112: }
3113: return LL;
3114: }
3115:
3116: def ideal_uniq(L) /* sub procedure of welldec and normposdec */
3117: {
3118: for (R = [],I = 0,N=length(L); I < N; I++) {
3119: if ( R == [] )
3120: R = append(R,[L[I]]);
3121: else {
3122: for (J = 0; J < length(R); J++)
3123: if ( gb_comp(L[I],R[J]) )
3124: break;
3125: if ( J == length(R) )
3126: R = append(R,[L[I]]);
3127: }
3128: }
3129: return R;
3130: }
3131:
3132: def ideal_uniq_by_first(L) /* sub procedure of welldec and normposdec */
3133: {
3134: for (R = [],I = 0,N=length(L); I < N; I++) {
3135: if ( R == [] )
3136: R = append(R,[L[I]]);
3137: else {
3138: for (J = 0; J < length(R); J++)
3139: if ( gb_comp(L[I][0],R[J][0]) )
3140: break;
3141: if ( J == length(R) )
3142: R = append(R,[L[I]]);
3143: }
3144: }
3145: return R;
3146: }
3147:
3148: def prime_irred_sf(TP,VL,Ord)
3149: {
3150: TP = ideal_uniq(TP);
3151: N = length(TP);
3152: for (P = [], I = 0; I < N; I++) {
3153: for (J = 0; J < N; J++)
3154: if ( I != J && inclusion_test(TP[J],TP[I],VL,Ord) )
3155: break;
3156: if (J == N)
3157: P = append(P,[TP[I]]);
3158: }
3159: return P;
3160: }
3161:
3162: def prime_irred_sf_by_first(TP,VL,Ord)
3163: {
3164: TP = ideal_uniq_by_first(TP);
3165: N = length(TP);
3166: for (P = [], I = 0; I < N; I++) {
3167: for (J = 0; J < N; J++)
3168: if ( I != J && inclusion_test(car(TP[J]),car(TP[I]),VL,Ord) )
3169: break;
3170: if (J == N)
3171: P = append(P,[TP[I]]);
3172: }
3173: return P;
3174: }
3175:
3176: def monic_sf_first(L,V)
3177: {
3178: for ( R = [], T = L; T != []; T = cdr(T) )
3179: R = cons([monic_sf(car(T)[0],V),car(T)[1]],R);
3180: return reverse(R);
3181: }
3182:
3183: def monic_sf(P,V)
3184: {
3185: if ( type(P) == 4 )
3186: return map(monic_sf,P,V);
3187: else {
3188: D = dp_ptod(P,V);
3189: return dp_dtop(D/dp_hc(D),V);
3190: }
3191: }
3192:
3193: def gr_fctr_sf(FL,VL,Ord)
3194: {
3195: T0 = time();
3196: COMMONCHECK_SF = 1;
3197: for (TP = [],I = 0; I<length(FL); I++ ) {
3198: F = FL[I];
3199: SF = idealsqfr_sf(F);
3200: if ( !gb_comp(F,SF) )
3201: F = dp_gr_f_main(SF,VL,0,Ord);
3202: CID_SF=[1];
3203: SP = gr_fctr_sub_sf(F,VL,Ord);
3204: TP = append(TP,SP);
3205: }
3206: TP = prime_irred_sf(TP,VL,Ord);
3207: T_GRF += time()[3]-T0[3];
3208: return TP;
3209: }
3210:
3211: def gr_fctr_sub_sf(G,VL,Ord)
3212: {
3213: if ( length(G) == 1 && type(G[0]) == 1 )
3214: return [G];
3215: RL = [];
3216: for (I = 0; I < length(G); I++) {
3217: FL = sffctr(G[I]); L = length(FL); N = FL[1][1];
3218: if (L > 2 || N > 1) {
3219: TLL = [];
3220: for (J = 1; J < L; J++) {
3221: W = cons(FL[J][0],G);
3222: NG = dp_gr_f_main(W,VL,0,Ord);
3223: TNG = idealsqfr_sf(NG);
3224: if ( !gb_comp(NG,TNG) )
3225: NG = dp_gr_f_main(TNG,VL,0,Ord);
3226: if ( !inclusion_test(CID_SF,NG,VL,Ord) ) {
3227: DG = gr_fctr_sub_sf(NG,VL,Ord);
3228: RL = append(RL,DG);
3229: if ( J <= L-2
3230: && !inclusion_test(CID_SF,NG,VL,Ord)
3231: && COMMONCHECK_SF ) {
3232: CID_SF=ideal_intersection_sfrat(CID_SF,NG,VL);
3233: }
3234: }
3235: }
3236: break;
3237: }
3238: }
3239: if (I == length(G))
3240: RL = append([G],RL);
3241: return RL;
3242: }
3243: end$