Changeset 597783 in git
- Timestamp:
- May 2, 2012, 1:41:01 PM (11 years ago)
- Branches:
- (u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
- Children:
- e2666785b249262bbfdcdc9e4ca1a39a4054e994
- Parents:
- cb7827a842fd2e9878498f940e9f7735fdbe488e
- git-author:
- Martin Lee <martinlee84@web.de>2012-05-02 13:41:01+02:00
- git-committer:
- Martin Lee <martinlee84@web.de>2012-05-07 14:16:48+02:00
- Location:
- factory
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
factory/cf_gcd.cc
rcb7827 r597783 1238 1238 CanonicalForm Dn, test= 0; 1239 1239 cl = gcd (f.lc(),g.lc()); 1240 CanonicalForm gcdcfcg= gcd (cf, cg); 1240 1241 CanonicalForm b= 1; 1241 1242 int minCommonDeg= 0; 1242 CanonicalForm gcdcfcg= gcd (cf, cg);1243 1243 for (i= tmax (f.level(), g.level()); i > 0; i--) 1244 1244 { … … 1271 1271 p = cf_getBigPrime( i ); 1272 1272 i--; 1273 while ( i >= 0 && mod( cl , p ) == 0 )1273 while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 ) 1274 1274 { 1275 1275 p = cf_getBigPrime( i ); … … 1278 1278 //printf("try p=%d\n",p); 1279 1279 setCharacteristic( p ); 1280 Dp = gcd_poly( mapinto( f ), mapinto( g ) ); 1281 //Dp = GCD_small_p (mapinto (f), mapinto (g), cofp, cogp); 1280 Dp = GCD_small_p (mapinto (f), mapinto (g), cofp, cogp); 1282 1281 Dp /=Dp.lc(); 1283 cofp= mapinto (f)/Dp;1284 cogp= mapinto (g)/Dp;1285 1282 cofp /= lc (cofp); 1286 1283 cogp /= lc (cogp); … … 1348 1345 equal= true; 1349 1346 //Dn /=vcontent(Dn,Variable(1)); 1350 if ((equal || q > b) || ((cofn*Dn==f) && (cogn*Dn==g))) 1347 if (((abs (cofn*Dn)==abs (f)) && (abs (cogn*Dn)==abs (g))) || 1348 ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g))) 1351 1349 { 1352 1350 //printf(" -> success\n"); -
factory/cf_gcd_smallp.cc
rcb7827 r597783 425 425 } 426 426 427 CanonicalForm 428 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 429 CanonicalForm& coF, CanonicalForm& coG, 430 Variable & alpha, CFList& l, bool& topLevel); 431 432 CanonicalForm 433 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 434 Variable & alpha, CFList& l, bool& topLevel) 435 { 436 CanonicalForm dummy1, dummy2; 437 CanonicalForm result= GCD_Fp_extension (F, G, dummy1, dummy2, alpha, l, 438 topLevel); 439 return result; 440 } 441 427 442 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ , 428 443 /// l and topLevel are only used internally, output is monic … … 431 446 CanonicalForm 432 447 GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 448 CanonicalForm& coF, CanonicalForm& coG, 433 449 Variable & alpha, CFList& l, bool& topLevel) 434 450 { 435 451 CanonicalForm A= F; 436 452 CanonicalForm B= G; 437 if (F.isZero() && degree(G) > 0) return G/Lc(G); 438 else if (G.isZero() && degree (F) > 0) return F/Lc(F); 439 else if (F.isZero() && G.isZero()) return F.genOne(); 440 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne(); 441 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F); 442 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); 443 if (F == G) return F/Lc(F); 453 if (F.isZero() && degree(G) > 0) 454 { 455 coF= 0; 456 coG= Lc (G); 457 return G/Lc(G); 458 } 459 else if (G.isZero() && degree (F) > 0) 460 { 461 coF= Lc (F); 462 coG= 0; 463 return F/Lc(F); 464 } 465 else if (F.isZero() && G.isZero()) 466 { 467 coF= coG= 0; 468 return F.genOne(); 469 } 470 if (F.inBaseDomain() || G.inBaseDomain()) 471 { 472 coF= F; 473 coG= G; 474 return F.genOne(); 475 } 476 if (F.isUnivariate() && fdivides(F, G, coG)) 477 { 478 coF= Lc (F); 479 return F/Lc(F); 480 } 481 if (G.isUnivariate() && fdivides(G, F, coF)) 482 { 483 coG= Lc (G); 484 return G/Lc(G); 485 } 486 if (F == G) 487 { 488 coF= coG= Lc (F); 489 return F/Lc(F); 490 } 444 491 445 492 CFMap M,N; 446 493 int best_level= myCompress (A, B, M, N, topLevel); 447 494 448 if (best_level == 0) return B.genOne(); 495 if (best_level == 0) 496 { 497 coF= F; 498 coG= G; 499 return B.genOne(); 500 } 449 501 450 502 A= M(A); … … 455 507 //univariate case 456 508 if (A.isUnivariate() && B.isUnivariate()) 457 return N (gcd(A,B)); 509 { 510 CanonicalForm result= gcd (A, B); 511 coF= N (A/result); 512 coG= N (B/result); 513 return N (result); 514 } 458 515 459 516 CanonicalForm cA, cB; // content of A and B … … 484 541 if (compressConvexDense) 485 542 { 543 CanonicalForm bufcA= cA; 544 CanonicalForm bufcB= cB; 486 545 cA= content (ppA, 1); 487 546 cB= content (ppB, 1); … … 489 548 ppB /= cB; 490 549 gcdcAcB *= gcd (cA, cB); 550 cA *= bufcA; 551 cB *= bufcB; 491 552 if (ppA.isUnivariate() || ppB.isUnivariate()) 492 553 { 493 554 if (ppA.level() == ppB.level()) 494 return N (gcd (ppA, ppB)*gcdcAcB); 555 { 556 CanonicalForm result= gcd (ppA, ppB); 557 coF= N ((ppA/result)*(cA/gcdcAcB)); 558 coG= N ((ppB/result)*(cB/gcdcAcB)); 559 return N (result*gcdcAcB); 560 } 495 561 else 562 { 563 coF= N (ppA*(cA/gcdcAcB)); 564 coG= N (ppB*(cB/gcdcAcB)); 496 565 return N (gcdcAcB); 566 } 497 567 } 498 568 … … 511 581 { 512 582 if (ppA.level() == ppB.level()) 513 return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB); 583 { 584 CanonicalForm result= gcd (ppA, ppB); 585 coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB)); 586 coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB)); 587 return N (decompress (result, MM, V)*gcdcAcB); 588 } 514 589 else 590 { 591 coF= N (decompress (ppA, MM, V)); 592 coG= N (decompress (ppB, MM, V)); 515 593 return N (gcdcAcB); 594 } 516 595 } 517 596 } … … 523 602 lcB= uni_lcoeff (ppB); 524 603 525 if (fdivides (lcA, lcB))604 /*if (fdivides (lcA, lcB)) 526 605 { 527 606 if (fdivides (A, B)) … … 532 611 if (fdivides (B, A)) 533 612 return G/Lc(G); 534 } 613 }*/ 535 614 536 615 gcdlcAlcB= gcd (lcA, lcB); … … 539 618 540 619 if (d == 0) 620 { 621 coF= N (ppA*(cA/gcdcAcB)); 622 coG= N (ppB*(cB/gcdcAcB)); 541 623 return N(gcdcAcB); 624 } 542 625 543 626 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); … … 545 628 d= d0; 546 629 if (d == 0) 630 { 631 coF= N (ppA*(cA/gcdcAcB)); 632 coG= N (ppB*(cB/gcdcAcB)); 547 633 return N(gcdcAcB); 634 } 548 635 549 636 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; 550 CanonicalForm newtonPoly; 637 CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m, 638 coG_m, ppCoF, ppCoG; 551 639 552 640 newtonPoly= 1; 553 641 m= gcdlcAlcB; 554 642 G_m= 0; 643 coF= 0; 644 coG= 0; 555 645 H= 0; 556 646 bool fail= false; … … 560 650 CanonicalForm prim_elem, im_prim_elem; 561 651 CFList source, dest; 562 int bound = tmin (degree (ppA, 1), degree (ppB, 1));563 int count= 0;652 int bound1= degree (ppA, 1); 653 int bound2= degree (ppB, 1); 564 654 do 565 655 { 566 random_element= randomElement (m , V_buf, l, fail);656 random_element= randomElement (m*lcA*lcB, V_buf, l, fail); 567 657 if (fail) 568 658 { … … 579 669 { 580 670 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 581 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 671 G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest); 672 coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest); 673 coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest); 582 674 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 583 675 source, dest); … … 586 678 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, 587 679 source, dest); 680 lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest); 681 lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest); 588 682 for (CFListIterator i= l; i.hasItem(); i++) 589 683 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, … … 605 699 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 606 700 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 701 coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 702 coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 607 703 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 608 704 source, dest); … … 611 707 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 612 708 source, dest); 709 lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 710 lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 613 711 614 712 fail= false; 615 random_element= randomElement (m , V_buf, l, fail );713 random_element= randomElement (m*lcA*lcB, V_buf, l, fail ); 616 714 DEBOUTLN (cerr, "fail= " << fail); 617 715 CFList list; 618 716 TIMING_START (gcd_recursion); 619 717 G_random_element= 620 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 718 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), 719 coF_random_element, coG_random_element, V_buf, 621 720 list, topLevel); 622 721 TIMING_END_AND_PRINT (gcd_recursion, … … 629 728 TIMING_START (gcd_recursion); 630 729 G_random_element= 631 GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf, 730 GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), 731 coF_random_element, coG_random_element, V_buf, 632 732 list, topLevel); 633 733 TIMING_END_AND_PRINT (gcd_recursion, … … 643 743 644 744 if (d0 == 0) 745 { 746 coF= N (ppA*(cA/gcdcAcB)); 747 coG= N (ppB*(cB/gcdcAcB)); 645 748 return N(gcdcAcB); 749 } 646 750 if (d0 > d) 647 751 { … … 654 758 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element)) 655 759 * G_random_element; 760 coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element)) 761 *coF_random_element; 762 coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element)) 763 *coG_random_element; 656 764 657 765 if (!G_random_element.inCoeffDomain()) … … 667 775 G_m= 0; 668 776 d= d0; 669 count= 0; 777 coF_m= 0; 778 coG_m= 0; 670 779 } 671 780 672 781 TIMING_START (newton_interpolation); 673 782 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); 783 coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x); 784 coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x); 674 785 TIMING_END_AND_PRINT (newton_interpolation, 675 786 "time for newton interpolation: "); 676 787 677 count++;678 788 //termination test 679 if (uni_lcoeff (H) == gcdlcAlcB) 680 { 681 cH= uni_content (H); 789 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H)) 790 { 791 if (gcdlcAlcB.isOne()) 792 cH= 1; 793 else 794 cH= uni_content (H); 682 795 ppH= H/cH; 796 CanonicalForm lcppH= gcdlcAlcB/cH; 797 CanonicalForm ccoF= lcA/lcppH; 798 ccoF /= Lc (ccoF); 799 CanonicalForm ccoG= lcB/lcppH; 800 ccoG /= Lc (ccoG); 801 ppCoF= coF/ccoF; 802 ppCoG= coG/ccoG; 683 803 if (inextension) 684 804 { 685 CFList u, v; 686 //maybe it's better to test if ppH is an element of F(\alpha) before 687 //mapping down 688 DEBOUTLN (cerr, "ppH before mapDown= " << ppH); 689 ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); 690 ppH /= Lc(ppH); 691 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 692 if ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB))) 693 { 805 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) && 806 (degree (ppCoG,1)+degree (ppH,1) == bound2) && 807 (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) || 808 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 809 { 810 CFList u, v; 694 811 DEBOUTLN (cerr, "ppH before mapDown= " << ppH); 695 812 ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); 813 ppCoF= mapDown (ppCoF, prim_elem, im_prim_elem, alpha, u, v); 814 ppCoF= mapDown (ppCoG, prim_elem, im_prim_elem, alpha, u, v); 696 815 ppH /= Lc(ppH); 697 816 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 698 817 if (compressConvexDense) 818 { 699 819 ppH= decompress (ppH, MM, V); 820 ppCoF= decompress (ppCoF, MM, V); 821 ppCoG= decompress (ppCoG, MM, V); 822 } 823 coF= N ((cA/gcdcAcB)*ppCoF); 824 coG= N ((cB/gcdcAcB)*ppCoG); 700 825 return N(gcdcAcB*ppH); 701 826 } 702 827 } 703 else if ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB))) 828 else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) && 829 (degree (ppCoG,1)+degree (ppH,1) == bound2) && 830 (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) || 831 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 704 832 { 705 833 if (compressConvexDense) 834 { 706 835 ppH= decompress (ppH, MM, V); 836 ppCoF= decompress (ppCoF, MM, V); 837 ppCoG= decompress (ppCoG, MM, V); 838 } 839 coF= N ((cA/gcdcAcB)*ppCoF); 840 coG= N ((cB/gcdcAcB)*ppCoG);; 707 841 return N(gcdcAcB*ppH); 708 842 } … … 710 844 711 845 G_m= H; 846 coF_m= coF; 847 coG_m= coG; 712 848 newtonPoly= newtonPoly*(x - random_element); 713 849 m= m*(x - random_element); … … 755 891 } 756 892 893 CanonicalForm 894 GCD_GF (const CanonicalForm& F, const CanonicalForm& G, 895 CanonicalForm& coF, CanonicalForm& coG, 896 CFList& l, bool& topLevel); 897 898 CanonicalForm 899 GCD_GF (const CanonicalForm& F, const CanonicalForm& G, CFList& l, 900 bool& topLevel) 901 { 902 CanonicalForm dummy1, dummy2; 903 CanonicalForm result= GCD_GF (F, G, dummy1, dummy2, l, topLevel); 904 return result; 905 } 906 757 907 /// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for 758 908 /// Computer Algebra" by Geddes, Czapor, Labahn … … 760 910 /// faster field arithmetics, however it might fail if the input is large since 761 911 /// the size of the base field is bounded by 2^16, output is monic 762 CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G, 912 CanonicalForm 913 GCD_GF (const CanonicalForm& F, const CanonicalForm& G, 914 CanonicalForm& coF, CanonicalForm& coG, 763 915 CFList& l, bool& topLevel) 764 916 { 765 917 CanonicalForm A= F; 766 918 CanonicalForm B= G; 767 if (F.isZero() && degree(G) > 0) return G/Lc(G); 768 else if (G.isZero() && degree (F) > 0) return F/Lc(F); 769 else if (F.isZero() && G.isZero()) return F.genOne(); 770 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne(); 771 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F); 772 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); 773 if (F == G) return F/Lc(F); 919 if (F.isZero() && degree(G) > 0) 920 { 921 coF= 0; 922 coG= Lc (G); 923 return G/Lc(G); 924 } 925 else if (G.isZero() && degree (F) > 0) 926 { 927 coF= Lc (F); 928 coG= 0; 929 return F/Lc(F); 930 } 931 else if (F.isZero() && G.isZero()) 932 { 933 coF= coG= 0; 934 return F.genOne(); 935 } 936 if (F.inBaseDomain() || G.inBaseDomain()) 937 { 938 coF= F; 939 coG= G; 940 return F.genOne(); 941 } 942 if (F.isUnivariate() && fdivides(F, G, coG)) 943 { 944 coF= Lc (F); 945 return F/Lc(F); 946 } 947 if (G.isUnivariate() && fdivides(G, F, coF)) 948 { 949 coG= Lc (G); 950 return G/Lc(G); 951 } 952 if (F == G) 953 { 954 coF= coG= Lc (F); 955 return F/Lc(F); 956 } 774 957 775 958 CFMap M,N; 776 959 int best_level= myCompress (A, B, M, N, topLevel); 777 960 778 if (best_level == 0) return B.genOne(); 961 if (best_level == 0) 962 { 963 coF= F; 964 coG= G; 965 return B.genOne(); 966 } 779 967 780 968 A= M(A); … … 785 973 //univariate case 786 974 if (A.isUnivariate() && B.isUnivariate()) 787 return N (gcd (A, B)); 975 { 976 CanonicalForm result= gcd (A, B); 977 coF= N (A/result); 978 coG= N (B/result); 979 return N (result); 980 } 788 981 789 982 CanonicalForm cA, cB; // content of A and B … … 814 1007 if (compressConvexDense) 815 1008 { 1009 CanonicalForm bufcA= cA; 1010 CanonicalForm bufcB= cB; 816 1011 cA= content (ppA, 1); 817 1012 cB= content (ppB, 1); … … 819 1014 ppB /= cB; 820 1015 gcdcAcB *= gcd (cA, cB); 1016 cA *= bufcA; 1017 cB *= bufcB; 821 1018 if (ppA.isUnivariate() || ppB.isUnivariate()) 822 1019 { 823 1020 if (ppA.level() == ppB.level()) 824 return N (gcd (ppA, ppB)*gcdcAcB); 1021 { 1022 CanonicalForm result= gcd (ppA, ppB); 1023 coF= N ((ppA/result)*(cA/gcdcAcB)); 1024 coG= N ((ppB/result)*(cB/gcdcAcB)); 1025 return N (result*gcdcAcB); 1026 } 825 1027 else 1028 { 1029 coF= N (ppA*(cA/gcdcAcB)); 1030 coG= N (ppB*(cB/gcdcAcB)); 826 1031 return N (gcdcAcB); 1032 } 827 1033 } 828 1034 … … 841 1047 { 842 1048 if (ppA.level() == ppB.level()) 843 return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB); 1049 { 1050 CanonicalForm result= gcd (ppA, ppB); 1051 coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB)); 1052 coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB)); 1053 return N (decompress (result, MM, V)*gcdcAcB); 1054 } 844 1055 else 1056 { 1057 coF= N (decompress (ppA, MM, V)); 1058 coG= N (decompress (ppB, MM, V)); 845 1059 return N (gcdcAcB); 1060 } 846 1061 } 847 1062 } … … 853 1068 lcB= uni_lcoeff (ppB); 854 1069 855 if (fdivides (lcA, lcB)) 856 { 857 if (fdivides (A, B)) 1070 /*if (fdivides (lcA, lcB)) 1071 { 1072 if (fdivides (ppA, ppB, coG)) 1073 { 1074 coF= 1; 1075 if (compressConvexDense) 1076 coG= decompress (coG, MM, V); 1077 coG= N (coG*(cB/gcdcAcB)); 858 1078 return F; 1079 } 859 1080 } 860 1081 if (fdivides (lcB, lcA)) 861 1082 { 862 if (fdivides (B, A)) 1083 if (fdivides (ppB, ppA, coF)) 1084 { 1085 coG= 1; 1086 if (compressConvexDense) 1087 coF= decompress (coF, MM, V); 1088 coF= N (coF*(cA/gcdcAcB)); 863 1089 return G; 864 } 1090 } 1091 }*/ 865 1092 866 1093 gcdlcAlcB= gcd (lcA, lcB); … … 868 1095 int d= totaldegree (ppA, Variable(2), Variable (ppA.level())); 869 1096 if (d == 0) 1097 { 1098 coF= N (ppA*(cA/gcdcAcB)); 1099 coG= N (ppB*(cB/gcdcAcB)); 870 1100 return N(gcdcAcB); 1101 } 871 1102 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level())); 872 1103 if (d0 < d) 873 1104 d= d0; 874 1105 if (d == 0) 1106 { 1107 coF= N (ppA*(cA/gcdcAcB)); 1108 coG= N (ppB*(cB/gcdcAcB)); 875 1109 return N(gcdcAcB); 1110 } 876 1111 877 1112 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; 878 CanonicalForm newtonPoly; 1113 CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m, 1114 coG_m, ppCoF, ppCoG; 879 1115 880 1116 newtonPoly= 1; 881 1117 m= gcdlcAlcB; 882 1118 G_m= 0; 1119 coF= 0; 1120 coG= 0; 883 1121 H= 0; 884 1122 bool fail= false; 885 topLevel= false;1123 //topLevel= false; 886 1124 bool inextension= false; 887 1125 int p=-1; … … 890 1128 int expon; 891 1129 char gf_name_buf= gf_name; 892 int bound = tmin (degree (ppA, 1), degree (ppB, 1));893 int count= 0;1130 int bound1= degree (ppA, 1); 1131 int bound2= degree (ppB, 1); 894 1132 do 895 1133 { 896 random_element= GFRandomElement (m , l, fail);1134 random_element= GFRandomElement (m*lcA*lcB, l, fail); 897 1135 if (fail) 898 1136 { … … 915 1153 G_m= GFMapUp (G_m, kk); 916 1154 newtonPoly= GFMapUp (newtonPoly, kk); 1155 coF_m= GFMapUp (coF_m, kk); 1156 coG_m= GFMapUp (coG_m, kk); 917 1157 ppA= GFMapUp (ppA, kk); 918 1158 ppB= GFMapUp (ppB, kk); 919 1159 gcdlcAlcB= GFMapUp (gcdlcAlcB, kk); 920 random_element= GFRandomElement (m, l, fail); 1160 lcA= GFMapUp (lcA, kk); 1161 lcB= GFMapUp (lcB, kk); 1162 random_element= GFRandomElement (m*lcA*lcB, l, fail); 921 1163 DEBOUTLN (cerr, "fail= " << fail); 922 1164 CFList list; 923 1165 TIMING_START (gcd_recursion); 924 1166 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 1167 coF_random_element, coG_random_element, 925 1168 list, topLevel); 926 1169 TIMING_END_AND_PRINT (gcd_recursion, … … 933 1176 TIMING_START (gcd_recursion); 934 1177 G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 1178 coF_random_element, coG_random_element, 935 1179 list, topLevel); 936 1180 TIMING_END_AND_PRINT (gcd_recursion, … … 949 1193 if (inextension) 950 1194 setCharacteristic (p, k, gf_name_buf); 1195 coF= N (ppA*(cA/gcdcAcB)); 1196 coG= N (ppB*(cB/gcdcAcB)); 951 1197 return N(gcdcAcB); 952 1198 } … … 961 1207 (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) * 962 1208 G_random_element; 1209 1210 coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element)) 1211 *coF_random_element; 1212 coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element)) 1213 *coG_random_element; 1214 963 1215 if (!G_random_element.inCoeffDomain()) 964 1216 d0= totaldegree (G_random_element, Variable(2), … … 973 1225 G_m= 0; 974 1226 d= d0; 975 count= 0; 1227 coF_m= 0; 1228 coG_m= 0; 976 1229 } 977 1230 978 1231 TIMING_START (newton_interpolation); 979 1232 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); 980 TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: "); 981 982 count++; 1233 coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x); 1234 coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x); 1235 TIMING_END_AND_PRINT (newton_interpolation, 1236 "time for newton interpolation: "); 1237 983 1238 //termination test 984 if (uni_lcoeff (H) == gcdlcAlcB) 985 { 986 cH= uni_content (H); 1239 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H)) 1240 { 1241 if (gcdlcAlcB.isOne()) 1242 cH= 1; 1243 else 1244 cH= uni_content (H); 987 1245 ppH= H/cH; 1246 CanonicalForm lcppH= gcdlcAlcB/cH; 1247 CanonicalForm ccoF= lcA/lcppH; 1248 ccoF /= Lc (ccoF); 1249 CanonicalForm ccoG= lcB/lcppH; 1250 ccoG /= Lc (ccoG); 1251 ppCoF= coF/ccoF; 1252 ppCoG= coG/ccoG; 988 1253 if (inextension) 989 1254 { 990 if ((count == bound) || (fdivides(ppH, ppA) && fdivides(ppH, ppB))) 1255 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) && 1256 (degree (ppCoG,1)+degree (ppH,1) == bound2) && 1257 (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) || 1258 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 991 1259 { 992 1260 DEBOUTLN (cerr, "ppH before mapDown= " << ppH); 993 1261 ppH= GFMapDown (ppH, k); 1262 ppCoF= GFMapDown (ppCoF, k); 1263 ppCoG= GFMapDown (ppCoG, k); 994 1264 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 995 1265 if (compressConvexDense) 1266 { 996 1267 ppH= decompress (ppH, MM, V); 1268 ppCoF= decompress (ppCoF, MM, V); 1269 ppCoG= decompress (ppCoG, MM, V); 1270 } 1271 coF= N ((cA/gcdcAcB)*ppCoF); 1272 coG= N ((cB/gcdcAcB)*ppCoG); 997 1273 setCharacteristic (p, k, gf_name_buf); 998 1274 return N(gcdcAcB*ppH); … … 1001 1277 else 1002 1278 { 1003 if (fdivides (ppH, ppA) && fdivides (ppH, ppB)) 1279 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) && 1280 (degree (ppCoG,1)+degree (ppH,1) == bound2) && 1281 (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) || 1282 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 1004 1283 { 1005 1284 if (compressConvexDense) 1285 { 1006 1286 ppH= decompress (ppH, MM, V); 1287 ppCoF= decompress (ppCoF, MM, V); 1288 ppCoG= decompress (ppCoG, MM, V); 1289 } 1290 coF= N ((cA/gcdcAcB)*ppCoF); 1291 coG= N ((cB/gcdcAcB)*ppCoG); 1007 1292 return N(gcdcAcB*ppH); 1008 1293 } … … 1011 1296 1012 1297 G_m= H; 1298 coF_m= coF; 1299 coG_m= coG; 1013 1300 newtonPoly= newtonPoly*(x - random_element); 1014 1301 m= m*(x - random_element); … … 1068 1355 } 1069 1356 1070 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1071 bool& topLevel, CFList& l) 1357 CanonicalForm 1358 GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1359 CanonicalForm& coF, CanonicalForm& coG, 1360 bool& topLevel, CFList& l); 1361 1362 CanonicalForm 1363 GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1364 bool& topLevel, CFList& l) 1365 { 1366 CanonicalForm dummy1, dummy2; 1367 CanonicalForm result= GCD_small_p (F, G, dummy1, dummy2, topLevel, l); 1368 return result; 1369 } 1370 1371 CanonicalForm 1372 GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, 1373 CanonicalForm& coF, CanonicalForm& coG, 1374 bool& topLevel, CFList& l) 1072 1375 { 1073 1376 CanonicalForm A= F; 1074 1377 CanonicalForm B= G; 1075 if (F.isZero() && degree(G) > 0) return G/Lc(G); 1076 else if (G.isZero() && degree (F) > 0) return F/Lc(F); 1077 else if (F.isZero() && G.isZero()) return F.genOne(); 1078 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne(); 1079 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F); 1080 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G); 1081 if (F == G) return F/Lc(F); 1082 1378 if (F.isZero() && degree(G) > 0) 1379 { 1380 coF= 0; 1381 coG= Lc (G); 1382 return G/Lc(G); 1383 } 1384 else if (G.isZero() && degree (F) > 0) 1385 { 1386 coF= Lc (F); 1387 coG= 0; 1388 return F/Lc(F); 1389 } 1390 else if (F.isZero() && G.isZero()) 1391 { 1392 coF= coG= 0; 1393 return F.genOne(); 1394 } 1395 if (F.inBaseDomain() || G.inBaseDomain()) 1396 { 1397 coF= F; 1398 coG= G; 1399 return F.genOne(); 1400 } 1401 if (F.isUnivariate() && fdivides(F, G, coG)) 1402 { 1403 coF= Lc (F); 1404 return F/Lc(F); 1405 } 1406 if (G.isUnivariate() && fdivides(G, F, coF)) 1407 { 1408 coG= Lc (G); 1409 return G/Lc(G); 1410 } 1411 if (F == G) 1412 { 1413 coF= coG= Lc (F); 1414 return F/Lc(F); 1415 } 1083 1416 CFMap M,N; 1084 1417 int best_level= myCompress (A, B, M, N, topLevel); 1085 1418 1086 if (best_level == 0) return B.genOne(); 1419 if (best_level == 0) 1420 { 1421 coF= F; 1422 coG= G; 1423 return B.genOne(); 1424 } 1087 1425 1088 1426 A= M(A); … … 1093 1431 //univariate case 1094 1432 if (A.isUnivariate() && B.isUnivariate()) 1095 return N (gcd (A, B)); 1433 { 1434 CanonicalForm result= gcd (A, B); 1435 coF= N (A/result); 1436 coG= N (B/result); 1437 return N (result); 1438 } 1096 1439 1097 1440 CanonicalForm cA, cB; // content of A and B … … 1122 1465 if (compressConvexDense) 1123 1466 { 1467 CanonicalForm bufcA= cA; 1468 CanonicalForm bufcB= cB; 1124 1469 cA= content (ppA, 1); 1125 1470 cB= content (ppB, 1); … … 1127 1472 ppB /= cB; 1128 1473 gcdcAcB *= gcd (cA, cB); 1474 cA *= bufcA; 1475 cB *= bufcB; 1129 1476 if (ppA.isUnivariate() || ppB.isUnivariate()) 1130 1477 { 1131 1478 if (ppA.level() == ppB.level()) 1132 return N (gcd (ppA, ppB)*gcdcAcB); 1479 { 1480 CanonicalForm result= gcd (ppA, ppB); 1481 coF= N ((ppA/result)*(cA/gcdcAcB)); 1482 coG= N ((ppB/result)*(cB/gcdcAcB)); 1483 return N (result*gcdcAcB); 1484 } 1133 1485 else 1486 { 1487 coF= N (ppA*(cA/gcdcAcB)); 1488 coG= N (ppB*(cB/gcdcAcB)); 1134 1489 return N (gcdcAcB); 1490 } 1135 1491 } 1136 1492 … … 1149 1505 { 1150 1506 if (ppA.level() == ppB.level()) 1151 return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB); 1507 { 1508 CanonicalForm result= gcd (ppA, ppB); 1509 coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB)); 1510 coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB)); 1511 return N (decompress (result, MM, V)*gcdcAcB); 1512 } 1152 1513 else 1514 { 1515 coF= N (decompress (ppA, MM, V)); 1516 coG= N (decompress (ppB, MM, V)); 1153 1517 return N (gcdcAcB); 1518 } 1154 1519 } 1155 1520 } … … 1160 1525 lcB= uni_lcoeff (ppB); 1161 1526 1162 if (fdivides (lcA, lcB))1527 /*if (fdivides (lcA, lcB)) 1163 1528 { 1164 1529 if (fdivides (A, B)) … … 1169 1534 if (fdivides (B, A)) 1170 1535 return G/Lc(G); 1171 } 1536 }*/ 1172 1537 1173 1538 gcdlcAlcB= gcd (lcA, lcB); … … 1177 1542 1178 1543 if (d == 0) 1544 { 1545 coF= N (ppA*(cA/gcdcAcB)); 1546 coG= N (ppB*(cB/gcdcAcB)); 1179 1547 return N(gcdcAcB); 1548 } 1180 1549 1181 1550 d0= totaldegree (ppB, Variable (2), Variable (ppB.level())); … … 1185 1554 1186 1555 if (d == 0) 1556 { 1557 coF= N (ppA*(cA/gcdcAcB)); 1558 coG= N (ppB*(cB/gcdcAcB)); 1187 1559 return N(gcdcAcB); 1560 } 1188 1561 1189 1562 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH; 1190 CanonicalForm newtonPoly= 1; 1563 CanonicalForm newtonPoly, coF_random_element, coG_random_element, 1564 coF_m, coG_m, ppCoF, ppCoG; 1565 1566 newtonPoly= 1; 1191 1567 m= gcdlcAlcB; 1192 1568 H= 0; 1569 coF= 0; 1570 coG= 0; 1193 1571 G_m= 0; 1194 1572 Variable alpha, V_buf; … … 1196 1574 bool inextension= false; 1197 1575 bool inextensionextension= false; 1198 topLevel= false;1576 bool topLevel2= false; 1199 1577 CFList source, dest; 1200 int bound = tmin (degree (ppA, 1), degree (ppB, 1));1201 int count= 0;1578 int bound1= degree (ppA, 1); 1579 int bound2= degree (ppB, 1); 1202 1580 do 1203 1581 { 1204 1582 if (inextension) 1205 random_element= randomElement (m , V_buf, l, fail);1583 random_element= randomElement (m*lcA*lcB, V_buf, l, fail); 1206 1584 else 1207 random_element= FpRandomElement (m , l, fail);1585 random_element= FpRandomElement (m*lcA*lcB, l, fail); 1208 1586 1209 1587 if (!fail && !inextension) … … 1212 1590 TIMING_START (gcd_recursion); 1213 1591 G_random_element= 1214 GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 1215 list); 1592 GCD_small_p (ppA (random_element,x), ppB (random_element,x), 1593 coF_random_element, coG_random_element, topLevel2, 1594 list); 1216 1595 TIMING_END_AND_PRINT (gcd_recursion, 1217 1596 "time for recursive call: "); … … 1223 1602 TIMING_START (gcd_recursion); 1224 1603 G_random_element= 1225 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1226 list, topLevel); 1604 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), 1605 coF_random_element, coG_random_element, alpha, 1606 list, topLevel2); 1227 1607 TIMING_END_AND_PRINT (gcd_recursion, 1228 1608 "time for recursive call: "); … … 1241 1621 inextension= true; 1242 1622 fail= false; 1243 random_element= randomElement (m , alpha, l, fail);1623 random_element= randomElement (m*lcA*lcB, alpha, l, fail); 1244 1624 deg++; 1245 1625 } while (fail); … … 1248 1628 TIMING_START (gcd_recursion); 1249 1629 G_random_element= 1250 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 1251 list, topLevel); 1630 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), 1631 coF_random_element, coG_random_element, alpha, 1632 list, topLevel2); 1252 1633 TIMING_END_AND_PRINT (gcd_recursion, 1253 1634 "time for recursive call: "); … … 1269 1650 { 1270 1651 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 1271 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest); 1652 G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest); 1653 coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest); 1654 coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest); 1272 1655 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, 1273 1656 source, dest); … … 1276 1659 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source, 1277 1660 dest); 1661 lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest); 1662 lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest); 1278 1663 for (CFListIterator i= l; i.hasItem(); i++) 1279 1664 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha, … … 1296 1681 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1297 1682 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1683 coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1684 coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1298 1685 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 1299 1686 source, dest); … … 1302 1689 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 1303 1690 source, dest); 1691 lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1692 lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest); 1304 1693 fail= false; 1305 random_element= randomElement (m , V_buf, l, fail );1694 random_element= randomElement (m*lcA*lcB, V_buf, l, fail ); 1306 1695 DEBOUTLN (cerr, "fail= " << fail); 1307 1696 CFList list; 1308 1697 TIMING_START (gcd_recursion); 1309 1698 G_random_element= 1310 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf, 1311 list, topLevel); 1699 GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), 1700 coF_random_element, coG_random_element, V_buf, 1701 list, topLevel2); 1312 1702 TIMING_END_AND_PRINT (gcd_recursion, 1313 1703 "time for recursive call: "); … … 1322 1712 1323 1713 if (d0 == 0) 1714 { 1715 coF= N (ppA*(cA/gcdcAcB)); 1716 coG= N (ppB*(cB/gcdcAcB)); 1324 1717 return N(gcdcAcB); 1718 } 1719 1325 1720 if (d0 > d) 1326 1721 { … … 1333 1728 *G_random_element; 1334 1729 1730 coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element)) 1731 *coF_random_element; 1732 coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element)) 1733 *coG_random_element; 1335 1734 1336 1735 if (!G_random_element.inCoeffDomain()) … … 1346 1745 G_m= 0; 1347 1746 d= d0; 1348 count= 0; 1747 coF_m= 0; 1748 coG_m= 0; 1349 1749 } 1350 1750 1351 1751 TIMING_START (newton_interpolation); 1352 1752 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x); 1753 coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x); 1754 coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x); 1353 1755 TIMING_END_AND_PRINT (newton_interpolation, 1354 1756 "time for newton_interpolation: "); 1355 1757 1356 count++;1357 1758 //termination test 1358 if (uni_lcoeff (H) == gcdlcAlcB) 1359 { 1360 cH= uni_content (H); 1759 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H)) 1760 { 1761 if (gcdlcAlcB.isOne()) 1762 cH= 1; 1763 else 1764 cH= uni_content (H); 1361 1765 ppH= H/cH; 1362 1766 ppH /= Lc (ppH); 1767 CanonicalForm lcppH= gcdlcAlcB/cH; 1768 CanonicalForm ccoF= lcppH/Lc (lcppH); 1769 CanonicalForm ccoG= lcppH/Lc (lcppH); 1770 ppCoF= coF/ccoF; 1771 ppCoG= coG/ccoG; 1363 1772 DEBOUTLN (cerr, "ppH= " << ppH); 1364 if ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB))) 1773 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) && 1774 (degree (ppCoG,1)+degree (ppH,1) == bound2) && 1775 (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) || 1776 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 1365 1777 { 1366 1778 if (compressConvexDense) 1779 { 1367 1780 ppH= decompress (ppH, MM, V); 1781 ppCoF= decompress (ppCoF, MM, V); 1782 ppCoG= decompress (ppCoG, MM, V); 1783 } 1784 coF= N ((cA/gcdcAcB)*ppCoF); 1785 coG= N ((cB/gcdcAcB)*ppCoG); 1368 1786 return N(gcdcAcB*ppH); 1369 1787 } … … 1371 1789 1372 1790 G_m= H; 1791 coF_m= coF; 1792 coG_m= coG; 1373 1793 newtonPoly= newtonPoly*(x - random_element); 1374 1794 m= m*(x - random_element); … … 4452 4872 { 4453 4873 cand = DD[2] / content( DD[2], Variable(1) ); 4454 gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); 4874 gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); //TODO use multiplication instead of fdivides 4455 4875 4456 4876 if (passToGF && gcdfound) -
factory/cf_gcd_smallp.h
rcb7827 r597783 44 44 bool& top_level, CFList& l); 45 45 46 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& coF, CanonicalForm& coG, 47 bool& topLevel, CFList& l); 48 46 49 ///GCD of A and B over \f$ F_{p} \f$ 47 50 static inline CanonicalForm GCD_small_p (const CanonicalForm& A, const CanonicalForm& B) … … 50 53 bool top_level= true; 51 54 return GCD_small_p (A, B, top_level, list); 55 } 56 57 static inline CanonicalForm GCD_small_p (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& coA, CanonicalForm& coB) 58 { 59 CFList list; 60 bool top_level= true; 61 return GCD_small_p (A, B, coA, coB, top_level, list); 52 62 } 53 63
Note: See TracChangeset
for help on using the changeset viewer.