Changeset eb3abbb in git
 Timestamp:
 May 23, 2013, 7:20:14 PM (10 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a657104b677b4c461d018cbf3204d72d34ad66a9')
 Children:
 02d3d641de579bc515b51a2769fe96f9149cd06755150e5fbf40e9e0fa5d2c2fe5efd309f4ab2e12
 Parents:
 e2c9b2fd1bc99d892a2a72d1c9e5046a794ac330aba90e858d998d23e9f343db0210fff2a167df2f
 Location:
 factory
 Files:

 5 edited
Legend:
 Unmodified
 Added
 Removed

factory/cfNewtonPolygon.cc
re2c9b2 reb3abbb 499 499 } 500 500 501 #ifdef HAVE_NTL 502 void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A) 501 void mpz_mat_mul (const mpz_t* N, mpz_t*& M) 502 { 503 mpz_t * tmp= new mpz_t[4]; 504 505 mpz_init_set (tmp[0], N[0]); 506 mpz_mul (tmp[0], tmp[0], M[0]); 507 mpz_addmul (tmp[0], N[1], M[2]); 508 509 mpz_init_set (tmp[1], N[0]); 510 mpz_mul (tmp[1], tmp[1], M[1]); 511 mpz_addmul (tmp[1], N[1], M[3]); 512 513 mpz_init_set (tmp[2], N[2]); 514 mpz_mul (tmp[2], tmp[2], M[0]); 515 mpz_addmul (tmp[2], N[3], M[2]); 516 517 mpz_init_set (tmp[3], N[2]); 518 mpz_mul (tmp[3], tmp[3], M[1]); 519 mpz_addmul (tmp[3], N[3], M[3]); 520 521 mpz_set (M[0], tmp[0]); 522 mpz_set (M[1], tmp[1]); 523 mpz_set (M[2], tmp[2]); 524 mpz_set (M[3], tmp[3]); 525 526 mpz_clear (tmp[0]); 527 mpz_clear (tmp[1]); 528 mpz_clear (tmp[2]); 529 mpz_clear (tmp[3]); 530 531 delete [] tmp; 532 } 533 534 void mpz_mat_inv (mpz_t*& M) 535 { 536 mpz_t det; 537 mpz_init_set (det, M[0]); 538 mpz_mul (det, det, M[3]); 539 mpz_submul (det, M[1], M[2]); 540 541 mpz_t tmp; 542 mpz_init_set (tmp, M[0]); 543 mpz_divexact (tmp, tmp, det); 544 mpz_set (M[0], M[3]); 545 mpz_divexact (M[0], M[0], det); 546 mpz_set (M[3], tmp); 547 548 mpz_neg (M[1], M[1]); 549 mpz_divexact (M[1], M[1], det); 550 mpz_neg (M[2], M[2]); 551 mpz_divexact (M[2], M[2], det); 552 553 mpz_clear (det); 554 mpz_clear (tmp); 555 } 556 557 void convexDense(int** points, int sizePoints, mpz_t*& M, mpz_t*& A) 503 558 { 504 559 if (sizePoints < 3) … … 506 561 if (sizePoints == 2) 507 562 { 508 int maxX= (points [1] [1] < points [0] [1])?points [0] [1]:points [1] [1]; 509 int maxY= (points [1] [0] < points [0] [0])?points [0] [0]:points [1] [0]; 510 long u,v,g; 511 XGCD (g, u, v, maxX, maxY); 512 M.SetDims (2,2); 513 A.SetLength (2); 563 mpz_t u,v,g,maxX,maxY; 564 mpz_init (u); 565 mpz_init (v); 566 mpz_init (g); 567 mpz_init_set_si (maxX, 568 (points[1][1] < points[0][1])?points[0][1]:points[1][1]); 569 mpz_init_set_si (maxY, 570 (points[1][0] < points[0][0])?points[0][0]:points[1][0]); 571 mpz_gcdext (g, u, v, maxX, maxY); 514 572 if (points [0] [1] != points [0] [0] && points [1] [0] != points [1] [1]) 515 573 { 516 M (1,1)= u; 517 M (1,2)= v; 518 M (2,1)= maxY/g; 519 M (2,2)= maxX/g; 520 A (1)= u*maxX; 521 A (2)= (maxY/g)*maxX; 574 mpz_set (A[0], u); 575 mpz_mul (A[0], A[0], maxX); 576 mpz_set (M[2], maxY); 577 mpz_divexact (M[2], M[2], g); 578 mpz_set (A[1], M[2]); 579 mpz_neg (A[1], A[1]); 580 mpz_mul (A[1], A[1], maxX); 581 mpz_neg (u, u); 582 mpz_set (M[0], u); 583 mpz_set (M[1], v); 584 mpz_set (M[3], maxX); 585 mpz_divexact (M[3], M[3], g); 522 586 } 523 587 else 524 588 { 525 M (1,1)= u; 526 M (1,2)= v; 527 M (2,1)= maxY/g; 528 M (2,2)= maxX/g; 529 A (1)= to_ZZ (0); 530 A (2)= to_ZZ (0); 589 mpz_set (M[0], u); 590 mpz_set (M[1], v); 591 mpz_set (M[2], maxY); 592 mpz_divexact (M[2], M[2], g); 593 mpz_neg (M[2], M[2]); 594 mpz_set (M[3], maxX); 595 mpz_divexact (M[3], M[3], g); 531 596 } 597 mpz_clear (u); 598 mpz_clear (v); 599 mpz_clear (g); 600 mpz_clear (maxX); 601 mpz_clear (maxY); 532 602 } 533 603 else if (sizePoints == 1) 534 604 { 535 ident (M, 2); 536 A.SetLength (2); 537 A (1)= to_ZZ (0); 538 A (2)= to_ZZ (0); 605 mpz_set_si (M[0], 1); 606 mpz_set_si (M[3], 1); 539 607 } 540 608 return; 541 609 } 542 A.SetLength (2);543 A (1)= to_ZZ (0);544 A (2)= to_ZZ (0); 545 ident (M, 2);546 m at_ZZ Mu;547 Mu.SetDims (2, 2);548 Mu (2,1)= to_ZZ (1);549 Mu (1,2)= to_ZZ (1);550 Mu (1,1)= to_ZZ (0); 551 Mu (2,2)= to_ZZ (0);552 m at_ZZ Lambda;553 Lambda.SetDims (2, 2);554 Lambda (1,1)= to_ZZ (1);555 Lambda (1,2)= to_ZZ (1);556 Lambda (2,2)= to_ZZ (1); 557 Lambda (2,1)= to_ZZ (0);558 m at_ZZ InverseLambda;559 InverseLambda.SetDims (2,2);560 InverseLambda (1,1)= to_ZZ (1);561 InverseLambda (1,2)= to_ZZ (1);562 InverseLambda (2,2)= to_ZZ (1); 563 InverseLambda (2,1)= to_ZZ (0);564 ZZ tmp;610 mpz_set_si (M[0], 1); 611 mpz_set_si (M[3], 1); 612 613 mpz_t * Mu= new mpz_t[4]; 614 mpz_init_set_si (Mu[1], 1); 615 mpz_init_set_si (Mu[2], 1); 616 mpz_init (Mu[0]); 617 mpz_init (Mu[3]); 618 619 mpz_t * Lambda= new mpz_t[4]; 620 mpz_init_set_si (Lambda[0], 1); 621 mpz_init_set_si (Lambda[1], 1); 622 mpz_init_set_si (Lambda[3], 1); 623 mpz_init (Lambda[2]); 624 625 mpz_t * InverseLambda= new mpz_t[4]; 626 mpz_init_set_si (InverseLambda[0], 1); 627 mpz_init_set_si (InverseLambda[1], 1); 628 mpz_init_set_si (InverseLambda[3], 1); 629 mpz_init (InverseLambda[2]); 630 631 mpz_t tmp; 632 mpz_init (tmp); 565 633 int minDiff, minSum, maxDiff, maxSum, maxX, maxY, b, d, f, h; 566 getMaxMin 634 getMaxMin(points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY); 567 635 do 568 636 { … … 570 638 { 571 639 mu (points, sizePoints); 572 M= Mu*M; 573 tmp= A (1); 574 A (1)= A (2); 575 A (2)= tmp; 576 } 577 getMaxMin (points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY); 640 641 mpz_mat_mul (Mu, M); 642 643 mpz_set (tmp, A[0]); 644 mpz_set (A[0], A[1]); 645 mpz_set (A[1], tmp); 646 } 647 getMaxMin(points, sizePoints, minDiff, minSum, maxDiff, maxSum, maxX, maxY); 578 648 b= maxX  maxDiff; 579 649 d= maxX + maxY  maxSum; … … 584 654 lambda (points, sizePoints); 585 655 tau (points, sizePoints, maxY  f); 586 M= Lambda*M; 587 A [0] += (maxYf); 656 657 mpz_mat_mul (Lambda, M); 658 659 if (maxYf > 0) 660 mpz_add_ui (A[0], A[0], maxYf); 661 else 662 mpz_add_ui (A[0], A[0], fmaxY); 588 663 maxX= maxX + maxY  b  f; 589 664 } … … 592 667 lambdaInverse (points, sizePoints); 593 668 tau (points, sizePoints, h); 594 M= InverseLambda*M; 595 A [0] += (h); 669 670 mpz_mat_mul (InverseLambda, M); 671 672 if (h < 0) 673 mpz_add_ui (A[0], A[0], h); 674 else 675 mpz_sub_ui (A[0], A[0], h); 596 676 maxX= maxX + maxY  d  h; 597 677 } 598 678 else 679 { 680 mpz_clear (tmp); 681 mpz_clear (Mu[0]); 682 mpz_clear (Mu[1]); 683 mpz_clear (Mu[2]); 684 mpz_clear (Mu[3]); 685 delete [] Mu; 686 687 mpz_clear (Lambda[0]); 688 mpz_clear (Lambda[1]); 689 mpz_clear (Lambda[2]); 690 mpz_clear (Lambda[3]); 691 delete [] Lambda; 692 693 mpz_clear (InverseLambda[0]); 694 mpz_clear (InverseLambda[1]); 695 mpz_clear (InverseLambda[2]); 696 mpz_clear (InverseLambda[3]); 697 delete [] InverseLambda; 698 599 699 return; 700 } 600 701 } while (1); 601 702 } 602 703 603 704 CanonicalForm 604 compress (const CanonicalForm& F, m at_ZZ& M, vec_ZZ& A, bool computeMA)705 compress (const CanonicalForm& F, mpz_t*& M, mpz_t*& A, bool computeMA) 605 706 { 606 707 int n; … … 611 712 convexDense (newtonPolyg, n, M, A); 612 713 } 714 613 715 CanonicalForm result= 0; 614 ZZ expX, expY;615 716 Variable x= Variable (1); 616 717 Variable y= Variable (2); 617 718 618 ZZ minExpX, minExpY; 719 mpz_t expX, expY, minExpX, minExpY; 720 mpz_init (expX); 721 mpz_init (expY); 722 mpz_init (minExpX); 723 mpz_init (minExpY); 619 724 620 725 int k= 0; … … 624 729 if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha)) 625 730 { 626 expX= i.exp()*M (1,2) + A (1); 627 expY= i.exp()*M (2,2) + A (2); 731 mpz_set (expX, A[0]); 732 mpz_set (expY, A[1]); 733 mpz_addmul_ui (expX, M[1], i.exp()); 734 mpz_addmul_ui (expY, M[3], i.exp()); 735 628 736 if (k == 0) 629 737 { 630 m inExpY= expY;631 m inExpX= expX;738 mpz_set (minExpX, expX); 739 mpz_set (minExpY, expY); 632 740 k= 1; 633 741 } 634 742 else 635 743 { 636 minExpY= (minExpY > expY) ? expY : minExpY; 637 minExpX= (minExpX > expX) ? expX : minExpX; 744 if (mpz_cmp (minExpY, expY) > 0) 745 mpz_set (minExpY, expY); 746 if (mpz_cmp (minExpX, expX) > 0) 747 mpz_set (minExpX, expX); 638 748 } 639 result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 749 result += i.coeff()*power (x, mpz_get_si (expX))* 750 power (y, mpz_get_si (expY)); 640 751 continue; 641 752 } … … 643 754 if (k == 0) 644 755 { 645 expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1); 646 expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2); 647 minExpX= expX; 648 minExpY= expY; 649 result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 756 mpz_set (expX, A[0]); 757 mpz_addmul_ui (expX, M[1], i.exp()); 758 mpz_addmul_ui (expX, M[0], j.exp()); 759 760 mpz_set (expY, A[1]); 761 mpz_addmul_ui (expY, M[3], i.exp()); 762 mpz_addmul_ui (expY, M[2], j.exp()); 763 764 mpz_set (minExpX, expX); 765 mpz_set (minExpY, expY); 766 result += j.coeff()*power (x, mpz_get_si (expX))* 767 power (y, mpz_get_si (expY)); 650 768 j++; 651 769 k= 1; … … 654 772 for (; j.hasTerms(); j++) 655 773 { 656 expX= j.exp()*M (1,1) + i.exp()*M (1,2) + A (1); 657 expY= j.exp()*M (2,1) + i.exp()*M (2,2) + A (2); 658 result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 659 minExpY= (minExpY > expY) ? expY : minExpY; 660 minExpX= (minExpX > expX) ? expX : minExpX; 661 } 662 } 663 664 if (to_long (minExpX) < 0) 665 { 666 result *= power (x,to_long(minExpX)); 774 mpz_set (expX, A[0]); 775 mpz_addmul_ui (expX, M[1], i.exp()); 776 mpz_addmul_ui (expX, M[0], j.exp()); 777 778 mpz_set (expY, A[1]); 779 mpz_addmul_ui (expY, M[3], i.exp()); 780 mpz_addmul_ui (expY, M[2], j.exp()); 781 782 result += j.coeff()*power (x, mpz_get_si (expX))* 783 power (y, mpz_get_si (expY)); 784 if (mpz_cmp (minExpY, expY) > 0) 785 mpz_set (minExpY, expY); 786 if (mpz_cmp (minExpX, expX) > 0) 787 mpz_set (minExpX, expX); 788 } 789 } 790 791 if (mpz_sgn (minExpX) < 0) 792 { 793 result *= power (x,mpz_get_si(minExpX)); 667 794 result /= CanonicalForm (x, 0); 668 795 } 669 796 else 670 result /= power (x, to_long(minExpX));671 672 if ( to_long(minExpY) < 0)673 { 674 result *= power (y, to_long(minExpY));797 result /= power (x,mpz_get_si(minExpX)); 798 799 if (mpz_sgn (minExpY) < 0) 800 { 801 result *= power (y,mpz_get_si(minExpY)); 675 802 result /= CanonicalForm (y, 0); 676 803 } 677 804 else 678 result /= power (y, to_long(minExpY));805 result /= power (y,mpz_get_si(minExpY)); 679 806 680 807 CanonicalForm tmp= LC (result); … … 692 819 delete [] newtonPolyg [i]; 693 820 delete [] newtonPolyg; 694 M= inv (M); 695 } 821 mpz_mat_inv (M); 822 } 823 824 mpz_clear (expX); 825 mpz_clear (expY); 826 mpz_clear (minExpX); 827 mpz_clear (minExpY); 696 828 697 829 return result; … … 699 831 700 832 CanonicalForm 701 decompress (const CanonicalForm& F, const m at_ZZ& inverseM, const vec_ZZ&A)833 decompress (const CanonicalForm& F, const mpz_t* inverseM, const mpz_t * A) 702 834 { 703 835 CanonicalForm result= 0; 704 ZZ expX, expY;705 836 Variable x= Variable (1); 706 837 Variable y= Variable (2); 707 ZZ minExpX, minExpY; 838 839 mpz_t expX, expY, minExpX, minExpY; 840 mpz_init (expX); 841 mpz_init (expY); 842 mpz_init (minExpX); 843 mpz_init (minExpY); 844 708 845 if (F.isUnivariate() && F.level() == 1) 709 846 { 710 847 CFIterator i= F; 711 expX= (i.exp()  A (1))*inverseM (1,1) + (A (2))*inverseM (1,2); 712 expY= (i.exp()  A (1))*inverseM (2,1) + (A (2))*inverseM (2,2); 713 minExpX= expX; 714 minExpY= expY; 715 result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 848 849 mpz_set_si (expX, i.exp()); 850 mpz_sub (expX, expX, A[0]); 851 mpz_mul (expX, expX, inverseM[0]); 852 mpz_submul (expX, inverseM[1], A[1]); 853 854 mpz_set_si (expY, i.exp()); 855 mpz_sub (expY, expY, A[0]); 856 mpz_mul (expY, expY, inverseM[2]); 857 mpz_submul (expY, inverseM[3], A[1]); 858 859 mpz_set (minExpX, expX); 860 mpz_set (minExpY, expY); 861 result += i.coeff()*power (x, mpz_get_si (expX))* 862 power (y, mpz_get_si (expY)); 716 863 i++; 717 864 for (; i.hasTerms(); i++) 718 865 { 719 expX= (i.exp()  A (1))*inverseM (1,1) + (A (2))*inverseM (1,2); 720 expY= (i.exp()  A (1))*inverseM (2,1) + (A (2))*inverseM (2,2); 721 result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 722 minExpY= (minExpY > expY) ? expY : minExpY; 723 minExpX= (minExpX > expX) ? expX : minExpX; 724 } 725 726 if (to_long (minExpX) < 0) 727 { 728 result *= power (x,to_long(minExpX)); 866 mpz_set_si (expX, i.exp()); 867 mpz_sub (expX, expX, A[0]); 868 mpz_mul (expX, expX, inverseM[0]); 869 mpz_submul (expX, inverseM[1], A[1]); 870 871 mpz_set_si (expY, i.exp()); 872 mpz_sub (expY, expY, A[0]); 873 mpz_mul (expY, expY, inverseM[2]); 874 mpz_submul (expY, inverseM[3], A[1]); 875 876 result += i.coeff()*power (x, mpz_get_si (expX))* 877 power (y, mpz_get_si (expY)); 878 if (mpz_cmp (minExpY, expY) > 0) 879 mpz_set (minExpY, expY); 880 if (mpz_cmp (minExpX, expX) > 0) 881 mpz_set (minExpX, expX); 882 } 883 884 if (mpz_sgn (minExpX) < 0) 885 { 886 result *= power (x,mpz_get_si(minExpX)); 729 887 result /= CanonicalForm (x, 0); 730 888 } 731 889 else 732 result /= power (x, to_long(minExpX));733 734 if ( to_long(minExpY) < 0)735 { 736 result *= power (y, to_long(minExpY));890 result /= power (x,mpz_get_si(minExpX)); 891 892 if (mpz_sgn (minExpY) < 0) 893 { 894 result *= power (y,mpz_get_si(minExpY)); 737 895 result /= CanonicalForm (y, 0); 738 896 } 739 897 else 740 result /= power (y,to_long(minExpY)); 898 result /= power (y,mpz_get_si(minExpY)); 899 900 mpz_clear (expX); 901 mpz_clear (expY); 902 mpz_clear (minExpX); 903 mpz_clear (minExpY); 741 904 742 905 return result/ Lc (result); //normalize 743 906 } 744 907 908 mpz_t tmp; 909 mpz_init (tmp); 745 910 int k= 0; 746 911 Variable alpha; … … 749 914 if (i.coeff().inCoeffDomain() && hasFirstAlgVar (i.coeff(), alpha)) 750 915 { 751 expX= A(1)*inverseM (1,1) + (i.exp()  A (2))*inverseM (1,2); 752 expY= A(1)*inverseM (2,1) + (i.exp()  A (2))*inverseM (2,2); 916 mpz_set_si (expX, i.exp()); 917 mpz_sub (expX, expX, A[1]); 918 mpz_mul (expX, expX, inverseM[1]); 919 mpz_submul (expX, A[0], inverseM[0]); 920 921 mpz_set_si (expY, i.exp()); 922 mpz_sub (expY, expY, A[1]); 923 mpz_mul (expY, expY, inverseM[3]); 924 mpz_submul (expY, A[0], inverseM[2]); 925 753 926 if (k == 0) 754 927 { 755 m inExpY= expY;756 m inExpX= expX;928 mpz_set (minExpX, expX); 929 mpz_set (minExpY, expY); 757 930 k= 1; 758 931 } 759 932 else 760 933 { 761 minExpY= (minExpY > expY) ? expY : minExpY; 762 minExpX= (minExpX > expX) ? expX : minExpX; 934 if (mpz_cmp (minExpY, expY) > 0) 935 mpz_set (minExpY, expY); 936 if (mpz_cmp (minExpX, expX) > 0) 937 mpz_set (minExpX, expX); 763 938 } 764 result += i.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 939 result += i.coeff()*power (x, mpz_get_si (expX))* 940 power (y, mpz_get_si (expY)); 765 941 continue; 766 942 } … … 768 944 if (k == 0) 769 945 { 770 expX= (j.exp()  A (1))*inverseM (1,1) + (i.exp()  A (2))*inverseM (1,2); 771 expY= (j.exp()  A (1))*inverseM (2,1) + (i.exp()  A (2))*inverseM (2,2); 772 minExpX= expX; 773 minExpY= expY; 774 result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 946 mpz_set_si (expX, j.exp()); 947 mpz_sub (expX, expX, A[0]); 948 mpz_mul (expX, expX, inverseM[0]); 949 mpz_set_si (tmp, i.exp()); 950 mpz_sub (tmp, tmp, A[1]); 951 mpz_addmul (expX, tmp, inverseM[1]); 952 953 mpz_set_si (expY, j.exp()); 954 mpz_sub (expY, expY, A[0]); 955 mpz_mul (expY, expY, inverseM[2]); 956 mpz_set_si (tmp, i.exp()); 957 mpz_sub (tmp, tmp, A[1]); 958 mpz_addmul (expY, tmp, inverseM[3]); 959 960 mpz_set (minExpX, expX); 961 mpz_set (minExpY, expY); 962 result += j.coeff()*power (x, mpz_get_si (expX))* 963 power (y, mpz_get_si (expY)); 775 964 j++; 776 965 k= 1; … … 779 968 for (; j.hasTerms(); j++) 780 969 { 781 expX= (j.exp()  A (1))*inverseM (1,1) + (i.exp()  A (2))*inverseM (1,2); 782 expY= (j.exp()  A (1))*inverseM (2,1) + (i.exp()  A (2))*inverseM (2,2); 783 result += j.coeff()*power (x, to_long (expX))*power (y, to_long (expY)); 784 minExpY= (minExpY > expY) ? expY : minExpY; 785 minExpX= (minExpX > expX) ? expX : minExpX; 786 } 787 } 788 789 if (to_long (minExpX) < 0) 790 { 791 result *= power (x,to_long(minExpX)); 970 mpz_set_si (expX, j.exp()); 971 mpz_sub (expX, expX, A[0]); 972 mpz_mul (expX, expX, inverseM[0]); 973 mpz_set_si (tmp, i.exp()); 974 mpz_sub (tmp, tmp, A[1]); 975 mpz_addmul (expX, tmp, inverseM[1]); 976 977 mpz_set_si (expY, j.exp()); 978 mpz_sub (expY, expY, A[0]); 979 mpz_mul (expY, expY, inverseM[2]); 980 mpz_set_si (tmp, i.exp()); 981 mpz_sub (tmp, tmp, A[1]); 982 mpz_addmul (expY, tmp, inverseM[3]); 983 984 result += j.coeff()*power (x, mpz_get_si (expX))* 985 power (y, mpz_get_si (expY)); 986 if (mpz_cmp (minExpY, expY) > 0) 987 mpz_set (minExpY, expY); 988 if (mpz_cmp (minExpX, expX) > 0) 989 mpz_set (minExpX, expX); 990 } 991 } 992 993 if (mpz_sgn (minExpX) < 0) 994 { 995 result *= power (x,mpz_get_si(minExpX)); 792 996 result /= CanonicalForm (x, 0); 793 997 } 794 998 else 795 result /= power (x, to_long(minExpX));796 797 if ( to_long(minExpY) < 0)798 { 799 result *= power (y, to_long(minExpY));999 result /= power (x,mpz_get_si(minExpX)); 1000 1001 if (mpz_sgn (minExpY) < 0) 1002 { 1003 result *= power (y,mpz_get_si(minExpY)); 800 1004 result /= CanonicalForm (y, 0); 801 1005 } 802 1006 else 803 result /= power (y,to_long(minExpY)); 1007 result /= power (y,mpz_get_si(minExpY)); 1008 1009 mpz_clear (expX); 1010 mpz_clear (expY); 1011 mpz_clear (minExpX); 1012 mpz_clear (minExpY); 1013 mpz_clear (tmp); 804 1014 805 1015 return result/Lc (result); //normalize 806 1016 } 807 #endif808 1017 809 1018 //assumes the input is a Newton polygon of a bivariate polynomial which is 
factory/cfNewtonPolygon.h
re2c9b2 reb3abbb 16 16 17 17 // #include "config.h" 18 19 #ifdef HAVE_NTL20 #include "NTLconvert.h"21 #endif22 18 23 19 /// compute a polygon … … 104 100 ); 105 101 106 107 #ifdef HAVE_NTL108 102 /// Algorithm 5 as described in ConvexDense Bivariate Polynomial Factorization 109 103 /// by Berthomieu, Lecerf … … 111 105 ///< M (points)+A 112 106 int sizePoints,///< [in] size of points 113 m at_ZZ& M, ///< [in,out] returns an invertiblematrix114 vec_ZZ& A ///< [in,out] returns translation107 mpz_t*& M, ///< [in,out] returns an invertible 2x2 matrix 108 mpz_t*& A ///< [in,out] returns translation 115 109 ); 116 110 … … 122 116 compress (const CanonicalForm& F, ///< [in] compressed, i.e. F.level()==2, 123 117 ///< bivariate poly 124 m at_ZZ& inverseM, ///< [in,out] returns the inverse of M,118 mpz_t*& inverseM, ///< [in,out] returns the inverse of M, 125 119 ///< if computeMA==true, M otherwise 126 vec_ZZ& A, ///< [in,out] returns translation120 mpz_t*& A, ///< [in,out] returns translation 127 121 bool computeMA= true ///< [in] whether to compute M and A 128 122 ); … … 135 129 decompress (const CanonicalForm& F,///< [in] compressed, i.e. F.level()<= 2, 136 130 ///< uni or bivariate poly 137 const m at_ZZ&M, ///< [in] matrix M obtained from compress138 const vec_ZZ&A ///< [in] vector A obtained from compress131 const mpz_t* M, ///< [in] matrix M obtained from compress 132 const mpz_t* A ///< [in] vector A obtained from compress 139 133 ); 140 #endif141 134 142 135 #endif 
factory/cf_gcd_smallp.cc
re2c9b2 reb3abbb 537 537 ppB= B/cB; 538 538 539 int sizeNewtonPolyg;540 int ** newtonPolyg= NULL;541 mat_ZZ MM;542 vec_ZZ V;543 bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);544 if (compressConvexDense)545 {546 CanonicalForm bufcA= cA;547 CanonicalForm bufcB= cB;548 cA= content (ppA, 1);549 cB= content (ppB, 1);550 ppA /= cA;551 ppB /= cB;552 gcdcAcB *= gcd (cA, cB);553 cA *= bufcA;554 cB *= bufcB;555 if (ppA.isUnivariate()  ppB.isUnivariate())556 {557 if (ppA.level() == ppB.level())558 {559 CanonicalForm result= gcd (ppA, ppB);560 coF= N ((ppA/result)*(cA/gcdcAcB));561 coG= N ((ppB/result)*(cB/gcdcAcB));562 return N (result*gcdcAcB);563 }564 else565 {566 coF= N (ppA*(cA/gcdcAcB));567 coG= N (ppB*(cB/gcdcAcB));568 return N (gcdcAcB);569 }570 }571 572 newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);573 convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);574 575 for (int i= 0; i < sizeNewtonPolyg; i++)576 delete [] newtonPolyg[i];577 delete [] newtonPolyg;578 579 ppA= compress (ppA, MM, V, false);580 ppB= compress (ppB, MM, V, false);581 MM= inv (MM);582 583 if (ppA.isUnivariate() && ppB.isUnivariate())584 {585 if (ppA.level() == ppB.level())586 {587 CanonicalForm result= gcd (ppA, ppB);588 coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));589 coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));590 return N (decompress (result, MM, V)*gcdcAcB);591 }592 else593 {594 coF= N (decompress (ppA, MM, V));595 coG= N (decompress (ppB, MM, V));596 return N (gcdcAcB);597 }598 }599 }600 601 539 CanonicalForm lcA, lcB; // leading coefficients of A and B 602 540 CanonicalForm gcdlcAlcB; … … 604 542 lcA= uni_lcoeff (ppA); 605 543 lcB= uni_lcoeff (ppB); 606 607 /*if (fdivides (lcA, lcB))608 {609 if (fdivides (A, B))610 return F/Lc(F);611 }612 if (fdivides (lcB, lcA))613 {614 if (fdivides (B, A))615 return G/Lc(G);616 }*/617 544 618 545 gcdlcAlcB= gcd (lcA, lcB); … … 817 744 ppCoG= mapDown (ppCoG, prim_elem, im_prim_elem, alpha, u, v); 818 745 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 819 if (compressConvexDense)820 {821 ppH= decompress (ppH, MM, V);822 ppCoF= decompress (ppCoF, MM, V);823 ppCoG= decompress (ppCoG, MM, V);824 }825 746 coF= N ((cA/gcdcAcB)*ppCoF); 826 747 coG= N ((cB/gcdcAcB)*ppCoG); … … 835 756 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 836 757 { 837 if (compressConvexDense)838 {839 ppH= decompress (ppH, MM, V);840 ppCoF= decompress (ppCoF, MM, V);841 ppCoG= decompress (ppCoG, MM, V);842 }843 758 coF= N ((cA/gcdcAcB)*ppCoF); 844 759 coG= N ((cB/gcdcAcB)*ppCoG); … … 998 913 ppB= B/cB; 999 914 1000 int sizeNewtonPolyg;1001 int ** newtonPolyg= NULL;1002 mat_ZZ MM;1003 vec_ZZ V;1004 bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);1005 if (compressConvexDense)1006 {1007 CanonicalForm bufcA= cA;1008 CanonicalForm bufcB= cB;1009 cA= content (ppA, 1);1010 cB= content (ppB, 1);1011 ppA /= cA;1012 ppB /= cB;1013 gcdcAcB *= gcd (cA, cB);1014 cA *= bufcA;1015 cB *= bufcB;1016 if (ppA.isUnivariate()  ppB.isUnivariate())1017 {1018 if (ppA.level() == ppB.level())1019 {1020 CanonicalForm result= gcd (ppA, ppB);1021 coF= N ((ppA/result)*(cA/gcdcAcB));1022 coG= N ((ppB/result)*(cB/gcdcAcB));1023 return N (result*gcdcAcB);1024 }1025 else1026 {1027 coF= N (ppA*(cA/gcdcAcB));1028 coG= N (ppB*(cB/gcdcAcB));1029 return N (gcdcAcB);1030 }1031 }1032 1033 newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);1034 convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);1035 1036 for (int i= 0; i < sizeNewtonPolyg; i++)1037 delete [] newtonPolyg[i];1038 delete [] newtonPolyg;1039 1040 ppA= compress (ppA, MM, V, false);1041 ppB= compress (ppB, MM, V, false);1042 MM= inv (MM);1043 1044 if (ppA.isUnivariate() && ppB.isUnivariate())1045 {1046 if (ppA.level() == ppB.level())1047 {1048 CanonicalForm result= gcd (ppA, ppB);1049 coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));1050 coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));1051 return N (decompress (result, MM, V)*gcdcAcB);1052 }1053 else1054 {1055 coF= N (decompress (ppA, MM, V));1056 coG= N (decompress (ppB, MM, V));1057 return N (gcdcAcB);1058 }1059 }1060 }1061 1062 915 CanonicalForm lcA, lcB; // leading coefficients of A and B 1063 916 CanonicalForm gcdlcAlcB; … … 1065 918 lcA= uni_lcoeff (ppA); 1066 919 lcB= uni_lcoeff (ppB); 1067 1068 /*if (fdivides (lcA, lcB))1069 {1070 if (fdivides (ppA, ppB, coG))1071 {1072 coF= 1;1073 if (compressConvexDense)1074 coG= decompress (coG, MM, V);1075 coG= N (coG*(cB/gcdcAcB));1076 return F;1077 }1078 }1079 if (fdivides (lcB, lcA))1080 {1081 if (fdivides (ppB, ppA, coF))1082 {1083 coG= 1;1084 if (compressConvexDense)1085 coF= decompress (coF, MM, V);1086 coF= N (coF*(cA/gcdcAcB));1087 return G;1088 }1089 }*/1090 920 1091 921 gcdlcAlcB= gcd (lcA, lcB); … … 1119 949 H= 0; 1120 950 bool fail= false; 1121 //topLevel= false;951 topLevel= false; 1122 952 bool inextension= false; 1123 953 int p=1; … … 1265 1095 ppCoG= GFMapDown (ppCoG, k); 1266 1096 DEBOUTLN (cerr, "ppH after mapDown= " << ppH); 1267 if (compressConvexDense)1268 {1269 ppH= decompress (ppH, MM, V);1270 ppCoF= decompress (ppCoF, MM, V);1271 ppCoG= decompress (ppCoG, MM, V);1272 }1273 1097 coF= N ((cA/gcdcAcB)*ppCoF); 1274 1098 coG= N ((cB/gcdcAcB)*ppCoG); … … 1286 1110 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 1287 1111 { 1288 if (compressConvexDense)1289 {1290 ppH= decompress (ppH, MM, V);1291 ppCoF= decompress (ppCoF, MM, V);1292 ppCoG= decompress (ppCoG, MM, V);1293 }1294 1112 coF= N ((cA/gcdcAcB)*ppCoF); 1295 1113 coG= N ((cB/gcdcAcB)*ppCoG); … … 1459 1277 ppB= B/cB; 1460 1278 1461 int sizeNewtonPolyg;1462 int ** newtonPolyg= NULL;1463 mat_ZZ MM;1464 vec_ZZ V;1465 bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);1466 if (compressConvexDense)1467 {1468 CanonicalForm bufcA= cA;1469 CanonicalForm bufcB= cB;1470 cA= content (ppA, 1);1471 cB= content (ppB, 1);1472 ppA /= cA;1473 ppB /= cB;1474 gcdcAcB *= gcd (cA, cB);1475 cA *= bufcA;1476 cB *= bufcB;1477 if (ppA.isUnivariate()  ppB.isUnivariate())1478 {1479 if (ppA.level() == ppB.level())1480 {1481 CanonicalForm result= gcd (ppA, ppB);1482 coF= N ((ppA/result)*(cA/gcdcAcB));1483 coG= N ((ppB/result)*(cB/gcdcAcB));1484 return N (result*gcdcAcB);1485 }1486 else1487 {1488 coF= N (ppA*(cA/gcdcAcB));1489 coG= N (ppB*(cB/gcdcAcB));1490 return N (gcdcAcB);1491 }1492 }1493 1494 newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);1495 convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);1496 1497 for (int i= 0; i < sizeNewtonPolyg; i++)1498 delete [] newtonPolyg[i];1499 delete [] newtonPolyg;1500 1501 ppA= compress (ppA, MM, V, false);1502 ppB= compress (ppB, MM, V, false);1503 MM= inv (MM);1504 1505 if (ppA.isUnivariate() && ppB.isUnivariate())1506 {1507 if (ppA.level() == ppB.level())1508 {1509 CanonicalForm result= gcd (ppA, ppB);1510 coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));1511 coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));1512 return N (decompress (result, MM, V)*gcdcAcB);1513 }1514 else1515 {1516 coF= N (decompress (ppA, MM, V));1517 coG= N (decompress (ppB, MM, V));1518 return N (gcdcAcB);1519 }1520 }1521 }1522 1523 1279 CanonicalForm lcA, lcB; // leading coefficients of A and B 1524 1280 CanonicalForm gcdlcAlcB; 1525 1281 lcA= uni_lcoeff (ppA); 1526 1282 lcB= uni_lcoeff (ppB); 1527 1528 /*if (fdivides (lcA, lcB))1529 {1530 if (fdivides (A, B))1531 return F/Lc(F);1532 }1533 if (fdivides (lcB, lcA))1534 {1535 if (fdivides (B, A))1536 return G/Lc(G);1537 }*/1538 1283 1539 1284 gcdlcAlcB= gcd (lcA, lcB); … … 1776 1521 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG))) 1777 1522 { 1778 if (compressConvexDense)1779 {1780 ppH= decompress (ppH, MM, V);1781 ppCoF= decompress (ppCoF, MM, V);1782 ppCoG= decompress (ppCoG, MM, V);1783 }1784 1523 coF= N ((cA/gcdcAcB)*ppCoF); 1785 1524 coG= N ((cB/gcdcAcB)*ppCoG); 
factory/facBivar.h
re2c9b2 reb3abbb 83 83 return result; 84 84 } 85 mat_ZZ M; 86 vec_ZZ S; 85 86 mpz_t * M=new mpz_t [4]; 87 mpz_init (M[0]); 88 mpz_init (M[1]); 89 mpz_init (M[2]); 90 mpz_init (M[3]); 91 92 mpz_t * S=new mpz_t [2]; 93 mpz_init (S[0]); 94 mpz_init (S[1]); 95 87 96 F= compress (F, M, S); 88 97 CFList result= biFactorize (F, v); … … 98 107 result.insert (Lc (G)); 99 108 } 109 110 mpz_clear (M[0]); 111 mpz_clear (M[1]); 112 mpz_clear (M[2]); 113 mpz_clear (M[3]); 114 delete [] M; 115 116 mpz_clear (S[0]); 117 mpz_clear (S[1]); 118 delete [] S; 119 100 120 return result; 101 121 } … … 197 217 return result; 198 218 } 199 mat_ZZ M; 200 vec_ZZ S; 219 220 mpz_t * M=new mpz_t [4]; 221 mpz_init (M[0]); 222 mpz_init (M[1]); 223 mpz_init (M[2]); 224 mpz_init (M[3]); 225 226 mpz_t * S=new mpz_t [2]; 227 mpz_init (S[0]); 228 mpz_init (S[1]); 229 201 230 F= compress (F, M, S); 202 231 TIMING_START (fac_bi_sqrf); … … 233 262 result.insert (CFFactor (LcF, 1)); 234 263 } 264 265 mpz_clear (M[0]); 266 mpz_clear (M[1]); 267 mpz_clear (M[2]); 268 mpz_clear (M[3]); 269 delete [] M; 270 271 mpz_clear (S[0]); 272 mpz_clear (S[1]); 273 delete [] S; 274 235 275 return result; 236 276 } 
factory/facFqBivar.h
re2c9b2 reb3abbb 94 94 return result; 95 95 } 96 mat_ZZ M; 97 vec_ZZ S; 96 mpz_t * M=new mpz_t [4]; 97 mpz_init (M[0]); 98 mpz_init (M[1]); 99 mpz_init (M[2]); 100 mpz_init (M[3]); 101 102 mpz_t * S=new mpz_t [2]; 103 mpz_init (S[0]); 104 mpz_init (S[1]); 105 98 106 F= compress (F, M, S); 99 107 CFList result= biFactorize (F, info); … … 106 114 normalize (result); 107 115 result.insert (Lc(G)); 116 117 mpz_clear (M[0]); 118 mpz_clear (M[1]); 119 mpz_clear (M[2]); 120 mpz_clear (M[3]); 121 delete [] M; 122 123 mpz_clear (S[0]); 124 mpz_clear (S[1]); 125 delete [] S; 126 108 127 return result; 109 128 } … … 228 247 return result; 229 248 } 230 mat_ZZ M; 231 vec_ZZ S; 249 mpz_t * M=new mpz_t [4]; 250 mpz_init (M[0]); 251 mpz_init (M[1]); 252 mpz_init (M[2]); 253 mpz_init (M[3]); 254 255 mpz_t * S=new mpz_t [2]; 256 mpz_init (S[0]); 257 mpz_init (S[1]); 258 232 259 F= compress (F, M, S); 233 260 … … 254 281 normalize (result); 255 282 result.insert (CFFactor (LcF, 1)); 283 284 mpz_clear (M[0]); 285 mpz_clear (M[1]); 286 mpz_clear (M[2]); 287 mpz_clear (M[3]); 288 delete [] M; 289 290 mpz_clear (S[0]); 291 mpz_clear (S[1]); 292 delete [] S; 293 256 294 return result; 257 295 } … … 335 373 return result; 336 374 } 337 mat_ZZ M; 338 vec_ZZ S; 375 376 mpz_t * M=new mpz_t [4]; 377 mpz_init (M[0]); 378 mpz_init (M[1]); 379 mpz_init (M[2]); 380 mpz_init (M[3]); 381 382 mpz_t * S=new mpz_t [2]; 383 mpz_init (S[0]); 384 mpz_init (S[1]); 385 339 386 F= compress (F, M, S); 340 387 … … 361 408 normalize (result); 362 409 result.insert (CFFactor (LcF, 1)); 410 411 mpz_clear (M[0]); 412 mpz_clear (M[1]); 413 mpz_clear (M[2]); 414 mpz_clear (M[3]); 415 delete [] M; 416 417 mpz_clear (S[0]); 418 mpz_clear (S[1]); 419 delete [] S; 420 363 421 return result; 364 422 } … … 443 501 return result; 444 502 } 445 mat_ZZ M; 446 vec_ZZ S; 503 504 mpz_t * M=new mpz_t [4]; 505 mpz_init (M[0]); 506 mpz_init (M[1]); 507 mpz_init (M[2]); 508 mpz_init (M[3]); 509 510 mpz_t * S=new mpz_t [2]; 511 mpz_init (S[0]); 512 mpz_init (S[1]); 513 447 514 F= compress (F, M, S); 448 515 … … 469 536 normalize (result); 470 537 result.insert (CFFactor (LcF, 1)); 538 539 mpz_clear (M[0]); 540 mpz_clear (M[1]); 541 mpz_clear (M[2]); 542 mpz_clear (M[3]); 543 delete [] M; 544 545 mpz_clear (S[0]); 546 mpz_clear (S[1]); 547 delete [] S; 548 471 549 return result; 472 550 }
Note: See TracChangeset
for help on using the changeset viewer.