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