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