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