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