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