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