From 5af98f25dcee361fe8a431660eaaccc93b035a57 Mon Sep 17 00:00:00 2001 From: Michael Kerber Date: Fri, 27 Jun 2008 10:03:35 +0000 Subject: [PATCH] Moved code to compute subresultants and sturm-habicht-sequences into Polynomial-package --- .gitattributes | 4 + .../include/CGAL/Polynomial/bezout_matrix.h | 3 +- Polynomial/include/CGAL/Polynomial/fwd.h | 27 + .../CGAL/Polynomial/sturm_habicht_sequence.h | 396 +++++++++++++ .../include/CGAL/Polynomial/subresultants.h | 560 ++++++++++++++++++ .../Polynomial/sturm_habicht_sequence.cpp | 347 +++++++++++ Polynomial/test/Polynomial/subresultants.cpp | 472 +++++++++++++++ 7 files changed, 1807 insertions(+), 2 deletions(-) create mode 100644 Polynomial/include/CGAL/Polynomial/sturm_habicht_sequence.h create mode 100644 Polynomial/include/CGAL/Polynomial/subresultants.h create mode 100644 Polynomial/test/Polynomial/sturm_habicht_sequence.cpp create mode 100644 Polynomial/test/Polynomial/subresultants.cpp diff --git a/.gitattributes b/.gitattributes index fd6df64d893..38e14cd64d0 100644 --- a/.gitattributes +++ b/.gitattributes @@ -3145,7 +3145,11 @@ Polynomial/include/CGAL/Polynomial/may_have_common_factor.h -text Polynomial/include/CGAL/Polynomial/modular_filter.h -text Polynomial/include/CGAL/Polynomial/resultant.h -text Polynomial/include/CGAL/Polynomial/square_free_factorization.h -text +Polynomial/include/CGAL/Polynomial/sturm_habicht_sequence.h -text +Polynomial/include/CGAL/Polynomial/subresultants.h -text Polynomial/include/CGAL/Polynomial/univariate_polynomial_utils.h -text +Polynomial/test/Polynomial/sturm_habicht_sequence.cpp -text +Polynomial/test/Polynomial/subresultants.cpp -text Polytope_distance_d/doc_tex/Polytope_distance_d/dist.png -text Polytope_distance_d/doc_tex/Polytope_distance_d/polydist.eps -text svneol=unset#application/postscript Polytope_distance_d/doc_tex/Polytope_distance_d/polydist.gif -text svneol=unset#image/gif diff --git a/Polynomial/include/CGAL/Polynomial/bezout_matrix.h b/Polynomial/include/CGAL/Polynomial/bezout_matrix.h index 463e4e69fc0..9c1bac420c6 100644 --- a/Polynomial/include/CGAL/Polynomial/bezout_matrix.h +++ b/Polynomial/include/CGAL/Polynomial/bezout_matrix.h @@ -452,8 +452,6 @@ namespace CGALi { } return it; } -}// namespace CGALi - /*! \ingroup CGAL_resultant_matrix * \brief computes the coefficients of the polynomial subresultant sequence @@ -554,6 +552,7 @@ typename CGALi::Simple_matrix< NT> polynomial_subresultant_matrix( return Ret; } +} CGAL_END_NAMESPACE diff --git a/Polynomial/include/CGAL/Polynomial/fwd.h b/Polynomial/include/CGAL/Polynomial/fwd.h index 1cfaf0c7380..31ea2545dc2 100644 --- a/Polynomial/include/CGAL/Polynomial/fwd.h +++ b/Polynomial/include/CGAL/Polynomial/fwd.h @@ -73,8 +73,35 @@ inline int square_free_factorization_for_regular_polynomial(const Polynomial inline bool may_have_multiple_factor( const Polynomial&); template inline bool may_have_common_factor(const Polynomial&,const Polynomial&); +template< class Coeff > +struct Simple_matrix; + +template +CGALi::Simple_matrix polynomial_subresultant_matrix( + CGAL::Polynomial f, + CGAL::Polynomial g, + int d=0); + + +template inline + OutputIterator polynomial_subresultants(CGAL::Polynomial A, CGAL::Polynomial B, + OutputIterator out); +template inline + OutputIterator principal_subresultants(CGAL::Polynomial A, CGAL::Polynomial B, + OutputIterator out); + +template inline + OutputIterator principal_sturm_habicht_sequence(CGAL::Polynomial A, + OutputIterator out); + +template OutputIterator + sturm_habicht_sequence(CGAL::Polynomial P, OutputIterator out); + } // namespace CGALi + + + } // namespace CGAL diff --git a/Polynomial/include/CGAL/Polynomial/sturm_habicht_sequence.h b/Polynomial/include/CGAL/Polynomial/sturm_habicht_sequence.h new file mode 100644 index 00000000000..fb7f116772c --- /dev/null +++ b/Polynomial/include/CGAL/Polynomial/sturm_habicht_sequence.h @@ -0,0 +1,396 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Michael Kerber +// +// ============================================================================ + +#ifndef CGAL_ACK_STURM_HABICHT +#define CGAL_ACK_STURM_HABICHT 1 + +#include +#include +#include +#include + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +/*! + * \brief compute the leading coefficients of the Sturm-Habicht sequence of + * the polynomial P + * + * The principal Sturm-Habicht sequence is obtained by computing the scalar + * subresultant sequence of P and its derivative, extended + * by P and P' and some sign changes. + * + * For details, see: Gonzalez-Vega,Lombardi,Recio,Roy: Determinants and Real + * Roots of Univariate Polynomials. Texts and Monographs in Symbolic + * Computation. Springer (1999) 300-316. + * Only the special case Q=1 is implemented + */ + template + OutputIterator prs_principal_sturm_habicht_sequence(CGAL::Polynomial P, + OutputIterator out) { + std::vector > stha; + CGAL::CGALi::sturm_habicht_sequence(P,std::back_inserter(stha)); + for(int i=0; i(stha.size()); i++) { + int d = stha[i].degree(); + CGAL_assertion(d<=i); + if(dP + * + * The principal Sturm-Habicht sequence is obtained by computing the scalar + * subresultant sequence of P and its derivative, extended + * by P and P' and some sign changes. + * + * For details, see: Gonzalez-Vega,Lombardi,Recio,Roy: Determinants and Real + * Roots of Univariate Polynomials. Texts and Monographs in Symbolic + * Computation. Springer (1999) 300-316. + * Only the special case Q=1 is implemented + */ + template + OutputIterator bezout_principal_sturm_habicht_sequence + (CGAL::Polynomial P, + OutputIterator out) { + + CGAL::Polynomial Px(P); + Px.diff(); + CGAL::CGALi::Simple_matrix M + = CGAL::CGALi::polynomial_subresultant_matrix(P,Px,1); + int n = static_cast(M.row_dimension()); + for(int i=0; i + void prs_first_two_sturm_habicht_coefficients(CGAL::Polynomial P, + OutputIterator1 pstha, + OutputIterator2 copstha) { + + std::vector > stha; + int n = P.degree(); + + sturm_habicht_sequence(P,std::back_inserter(stha)); + CGAL_assertion(static_cast(stha.size())==n+1); + for(int i=0;i<=n;i++) { + int d = stha[i].degree(); + CGAL_assertion(d<=i); + if(d + void bezout_first_two_sturm_habicht_coefficients(CGAL::Polynomial P, + OutputIterator1 pstha, + OutputIterator2 copstha) { + CGAL::Polynomial Px(P); + Px.diff(); + CGAL::CGALi::Simple_matrix M + = CGAL::CGALi::polynomial_subresultant_matrix(P,Px,2); + int n = static_cast(M.row_dimension()); + for(int i=0; i inline + void + first_two_sturm_habicht_coefficients_ + (CGAL::Polynomial A, + OutputIterator1 pstha, + OutputIterator2 copstha, + CGAL::Integral_domain_without_division_tag){ + + bezout_first_two_sturm_habicht_coefficients(A,pstha,copstha); + + } + + // the general function for CGAL::Unique_factorization_domain_tag + template inline + void + first_two_sturm_habicht_coefficients_ + (CGAL::Polynomial A, + OutputIterator1 pstha, + OutputIterator2 copstha, + CGAL::Unique_factorization_domain_tag) { + + return prs_first_two_sturm_habicht_coefficients(A,pstha,copstha); + + } + + template inline + void + first_two_sturm_habicht_coefficients_(CGAL::Polynomial A, + OutputIterator1 pstha, + OutputIterator2 copstha) { + typedef typename + CGAL::Algebraic_structure_traits::Algebraic_category + Algebraic_category; + first_two_sturm_habicht_coefficients_(A,pstha,copstha, + Algebraic_category()); + } + + // the general function for CGAL::Integral_domain_without_division_tag + template inline + OutputIterator + principal_sturm_habicht_sequence_ + (CGAL::Polynomial A, + OutputIterator out, + CGAL::Integral_domain_without_division_tag) { + + return bezout_principal_sturm_babicht_sequence(A,out); + } + + // the specialization for CGAL::Unique_factorization_domain_tag + template inline + OutputIterator + principal_sturm_habicht_sequence_ + (CGAL::Polynomial A, + OutputIterator out, + CGAL::Unique_factorization_domain_tag) { + + return prs_principal_sturm_habicht_sequence(A,out); + } + + template inline + OutputIterator principal_sturm_habicht_sequence_(CGAL::Polynomial A, + OutputIterator out) { + typedef typename + CGAL::Algebraic_structure_traits::Algebraic_category + Algebraic_category; + return principal_sturm_habicht_sequence_(A,out,Algebraic_category()); + } + + + + /*! \ingroup NiX_resultant_matrix + * \brief compute the sequence of + * principal Sturm-Habicht coefficients + */ + template inline + OutputIterator + principal_sturm_habicht_sequence(CGAL::Polynomial A, + OutputIterator out){ + + return CGAL::CGALi::principal_sturm_habicht_sequence_(A,out); + } + + /*! \ingroup NiX_resultant_matrix + * \brief computes the first two coefficients of each polynomial of + * the Sturm-Habicht sequence. + * + * This function is needed in Curve_analysis_2 for certain genericity checks + */ + template inline + void first_two_sturm_habicht_coefficients(CGAL::Polynomial A, + OutputIterator1 pstha, + OutputIterator2 copstha){ + + return CGAL::CGALi::first_two_sturm_habicht_coefficients_ + (A,pstha,copstha); + } + + /*! \ingroup NiX_resultant_matrix + * \brief compute the Sturm-Habicht sequence + */ + template OutputIterator + sturm_habicht_sequence(CGAL::Polynomial P, OutputIterator out) { + int p = P.degree(); + + CGAL::Polynomial P_x(P); + P_x.diff(); + + std::vector > stha; + + CGAL::CGALi::polynomial_subresultants(P,P_x,std::back_inserter(stha)); + stha.push_back(P); + + CGAL_assertion(static_cast(stha.size())==p+1); + + for(int i=0;i<=p; i++) { + if((p-i)%4==0 || (p-i)%4==1) { + *out++ = stha[i]; + } else { + *out++ = -stha[i]; + } + } + + return out; + } + + /*! \ingroup NiX_resultant_matrix + * \brief compute the Sturm-Habicht sequence with cofactors + */ + template OutputIterator1 + sturm_habicht_sequence_with_cofactors(CGAL::Polynomial P, + OutputIterator1 out_stha, + OutputIterator2 out_f, + OutputIterator3 out_fx) { + int p = P.degree(); + + CGAL::Polynomial P_x(P); + P_x.diff(); + + std::vector > stha,co_f,co_fx; + + CGAL::CGALi::prs_subresultants_with_cofactors(P,P_x, + std::back_inserter(stha), + std::back_inserter(co_f), + std::back_inserter(co_fx)); + stha.push_back(P); + co_f.push_back(CGAL::Polynomial(1)); + co_fx.push_back(CGAL::Polynomial(0)); + + CGAL_assertion(static_cast(stha.size())==p+1); + + for(int i=0;i<=p; i++) { + if((p-i)%4==0 || (p-i)%4==1) { + *out_stha++ = stha[i]; + *out_f++ = co_f[i]; + *out_fx++ = co_fx[i]; + } else { + *out_stha++ = -stha[i]; + *out_f++ = -co_f[i]; + *out_fx++ = -co_fx[i]; + } + } + + return out_stha; + } + + + /*! + * \brief returns the number of roots of a polynomial with given + * principal Sturm-Habicht sequence (counted without multiplicity) + */ + template + int stha_count_number_of_real_roots(InputIterator start,InputIterator end) { + if(start==end) { + return 0; + } + int m = 0; + + CGAL::Sign last_non_zero=CGAL::ZERO; //marks the starting point + CGAL::Sign curr_sign; + int k; + InputIterator el=start; + + //std::cout << "Sign of." << (*el) << std::endl; + + curr_sign=CGAL::sign(*el); + + while(curr_sign==CGAL::ZERO && el!=end) { + el++; + curr_sign=CGAL::sign(*el); + } + if(el==end) return 0; + + last_non_zero=curr_sign; + k=0; + el++; + + while(el!=end) { + curr_sign=CGAL::sign(*el); + + el++; + + if(curr_sign==CGAL::ZERO) { + k++; + } + else { + if(k%2==0) { // k is even + k=k/2; + int pm_one = (curr_sign==last_non_zero ? 1 : -1); + pm_one = (k%2==1) ? -pm_one : pm_one; + m+=pm_one; + } + k=0; + last_non_zero=curr_sign; + } + + } + return m; + } + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif // CGAL_ACK_STURM_HABICHT diff --git a/Polynomial/include/CGAL/Polynomial/subresultants.h b/Polynomial/include/CGAL/Polynomial/subresultants.h new file mode 100644 index 00000000000..308dfa2da23 --- /dev/null +++ b/Polynomial/include/CGAL/Polynomial/subresultants.h @@ -0,0 +1,560 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Michael Kerber +// +// ============================================================================ +#ifndef CGAL_ACK_SUBRESULTANTS_H +#define CGAL_ACK_SUBRESULTANTS_H + +#include + +#include +#include + +CGAL_BEGIN_NAMESPACE + + + namespace CGALi { + + // Intern function needed for Ducos algorithm + + template void lazard_optimization(NT y, + double n, + CGAL::Polynomial B, + CGAL::Polynomial& C) { + CGAL_precondition(n>0); + NT x = B.lcoeff(); + double a = pow(2.,std::floor(log(n)/log(2.))); + NT c = x; + n -= a; + while(a!=1) { + a/=2; + c=CGAL::integral_division(c*c,y); + if(n>=a) { + c=CGAL::integral_division(c*x,y); + n-=a; + } + } + C=c*B/y; + } + + template + void lickteig_roy_optimization(CGAL::Polynomial A, + CGAL::Polynomial B, + CGAL::Polynomial C, + NT s, + CGAL::Polynomial& D) { + typedef CGAL::Polynomial Poly; // for convenience + int d = A.degree(), e = B.degree(); + CGAL_precondition(d>=e); + std::vector H(d+1); + std::list initial; + initial.push_front(C.lcoeff()); + for(int i=0;i=e ? H[i][e] : NT(0); + H[i]-=(h_i_e*B)/B.lcoeff(); + initial.clear(); + std::copy(H[i].begin(),H[i].end(),std::back_inserter(initial)); + initial.push_front(NT(0)); + } + H[d]=Poly(initial.begin(),initial.end()); + D=Poly(0); + for(int i=0;i=e ? H[d][e] : NT(0); + D=(B.lcoeff()*(H[d]+D)-Hde*B)/s; + if((d-e)%2==0) { + D=-D; + } + return; + } + + template NT resultant_for_constant_polynomial + (CGAL::Polynomial P, CGAL::Polynomial Q) { + CGAL_assertion(P.degree() < 1 || Q.degree() < 1); + if(P.is_zero() || Q.is_zero() ) { + return NT(0); + } + if(P.degree()==0) { + return CGAL::ipower(P.lcoeff(),Q.degree()); + } else { + return CGAL::ipower(Q.lcoeff(),P.degree()); + } + } + + + /*! + * \brief Compute the sequence of subresultants with pseudo-division + * + * This is an implementation of Ducos' algorithm. It improves on the + * classical methods for subresultant computation by reducing the + * swell-up of intermediate results. For all details, see + * L.Ducos: Optimazations of the Subresultant algorithm. Journal of Pure + * and Applied Algebra 145 (2000) 149--163 + */ + template inline + OutputIterator prs_polynomial_subresultants(CGAL::Polynomial P, + CGAL::Polynomial Q, + OutputIterator out) { + + if(P.degree() < 1 || Q.degree() < 1) { + *out++ = CGAL::CGALi::resultant_for_constant_polynomial(P,Q); + return out; + } + + bool poly_swapped = (P.degree() < Q.degree()); + + if(poly_swapped) { + std::swap(P,Q); + } + + CGAL::Polynomial zero_pol(NT(0)); + std::vector > sres; + + int deg_diff=P.degree()-Q.degree(); + + if(deg_diff==0) { + sres.push_back(Q); + } else { + sres.push_back(CGAL::ipower(Q.lcoeff(),deg_diff-1)*Q); + } + + CGAL::Polynomial A,B,C,D,dummy_pol; + NT s,dummy_nt; + int delta,d,e; + + A=Q; + + s=CGAL::ipower(Q.lcoeff(),deg_diff); + + CGAL::Polynomial::pseudo_division(P, -Q, dummy_pol, B, dummy_nt); + + while(true) { + d=A.degree(); + e=B.degree(); + if(B.is_zero()) { + for(int i=0;i1) { + CGAL::CGALi::lazard_optimization(s,double(delta-1),B,C); + //C=CGAL::ipower(CGAL::integral_division(B.lcoeff(),s),delta-1)*B; + for(int i=0;i::pseudo_division(A, -B, dummy_pol, D, dummy_nt); + //B= D / (CGAL::ipower(s,delta)*A.lcoeff()); + A=C; + s=A.lcoeff(); + } + + CGAL_assertion(static_cast(sres.size()) + == Q.degree()+1); + + // If P and Q were swapped, correct the signs + if(poly_swapped) { + int p = P.degree(); + int q = Q.degree(); + for(int i=0;i<=q;i++) { + if((p-i)*(q-i) % 2 == 1) { + sres[q-i]=-sres[q-i]; + } + } + } + + // Now, reverse the entries + return std::copy(sres.rbegin(),sres.rend(),out); + } + + + /*! + * \brief Computes the polynomial subresultants + * as minors of the Bezout matrix + */ + template inline + OutputIterator bezout_polynomial_subresultants(CGAL::Polynomial P, + CGAL::Polynomial Q, + OutputIterator out) { + if(P.degree() < 1 || Q.degree() < 1) { + *out++ = CGAL::CGALi::resultant_for_constant_polynomial(P,Q); + return out; + } + + typedef CGAL::CGALi::Simple_matrix Matrix; + Matrix M = CGAL::CGALi::polynomial_subresultant_matrix(P,Q); + + int r = static_cast(M.row_dimension()); + + for(int i = 0;i < r; i++) { + std::vector curr_row; + std::copy(M[r-1-i].begin(), + M[r-1-i].end(), + std::back_inserter(curr_row)); + //std::reverse(curr_row.begin(),curr_row.end()); + *out++=CGAL::Polynomial(curr_row.rbegin(),curr_row.rend()); + } + int deg_diff=P.degree()-Q.degree(); + if(deg_diff==0) { + *out++=Q; + } else if(deg_diff>0) { + *out++=CGAL::ipower(Q.lcoeff(),deg_diff-1)*Q; + } else { + *out++=CGAL::ipower(P.lcoeff(),-deg_diff-1)*P; + } + + return out; + + } + + /*! + * \brief Compute the sequence of principal subresultants + * with pseudo-division + * + * Uses Ducos algorithm for the polynomial subresultant, and + * returns the formal leading coefficients. + */ + template inline + OutputIterator prs_principal_subresultants(CGAL::Polynomial P, + CGAL::Polynomial Q, + OutputIterator out) { + + std::vector > sres; + int q = std::min(Q.degree(),P.degree()); + + CGAL::CGALi::prs_polynomial_subresultants(P,Q,std::back_inserter(sres)); + CGAL_assertion(static_cast(sres.size()) == q+1); + for(int i=0; i <= q; i++) { + int d = sres[i].degree(); + CGAL_assertion(d<=i); + if(d inline + OutputIterator bezout_principal_subresultants(CGAL::Polynomial P, + CGAL::Polynomial Q, + OutputIterator out) { + if(P.degree() < 1 || Q.degree() < 1) { + *out++ = CGAL::CGALi::resultant_for_constant_polynomial(P,Q); + return out; + } + + typedef CGAL::CGALi::Simple_matrix Matrix; + Matrix M = CGAL::CGALi::polynomial_subresultant_matrix(P,Q,1); + + int r = static_cast(M.row_dimension()); + + for(int i = r - 1;i >=0; i--) { + *out++=M[i][i]; + } + int deg_diff=P.degree()-Q.degree(); + if(deg_diff==0) { + *out++=NT(1); + } else if(deg_diff>0) { + *out++=CGAL::ipower(Q.lcoeff(),deg_diff); + } else { + *out++=CGAL::ipower(P.lcoeff(),-deg_diff); + } + return out; + + } + + /*! + * \brief Computes the subresultants together with the according cofactors + */ + template + OutputIterator1 prs_subresultants_with_cofactors(CGAL::Polynomial P, + CGAL::Polynomial Q, + OutputIterator1 sres_out, + OutputIterator2 coP_out, + OutputIterator3 coQ_out) { + + + if(P.degree() < 1 || Q.degree() < 1) { + *sres_out++ = CGAL::CGALi::resultant_for_constant_polynomial(P,Q); + *coP_out++ = Q.lcoeff(); + *coQ_out++ = P.lcoeff(); + return sres_out; + } + + bool poly_swapped = (P.degree() < Q.degree()); + + if(poly_swapped) { + std::swap(P,Q); + } + + CGAL::Polynomial zero_pol(NT(0)); + std::vector > sres, coP, coQ; + + int deg_diff=P.degree()-Q.degree(); + + if(deg_diff==0) { + sres.push_back(Q); + } else { + sres.push_back(CGAL::ipower(Q.lcoeff(),deg_diff-1)*Q); + } + + + CGAL::Polynomial A,B,C,D,Quo, coPA, coPB, coQA, coQB, coPC, coQC; + NT s,m; + int delta,d,e; + + coPA = CGAL::Polynomial(NT(0)); + coQA = CGAL::Polynomial(CGAL::ipower(Q.lcoeff(),deg_diff-1)); + + coP.push_back(coPA); + coQ.push_back(coQA); + + A=Q; + + s=CGAL::ipower(Q.lcoeff(),deg_diff); + + CGAL::Polynomial::pseudo_division(P, -Q, Quo, B, m); + + coPB = CGAL::Polynomial(m); + coQB = Quo; + + + while(true) { + d=A.degree(); + e=B.degree(); + if(B.is_zero()) { + for(int i=0;i1) { + C=CGAL::ipower(B.lcoeff(),delta-1)*B / CGAL::ipower(s,delta-1); + + coPC = CGAL::ipower(B.lcoeff(),delta-1)*coPB / + CGAL::ipower(s,delta-1); + coQC = CGAL::ipower(B.lcoeff(),delta-1)*coQB / + CGAL::ipower(s,delta-1); + for(int i=0;i::pseudo_division(A, -B, Quo, D, m); + coPB = (m*coPA + Quo*coPB) / denominator; + coQB = (m*coQA + Quo*coQB) / denominator; + B = D / denominator; + A = C; + coPA = coPC; + coQA = coQC; + s = A.lcoeff(); + } + + CGAL_assertion(static_cast(sres.size()) + == Q.degree()+1); + + // If P and Q were swapped, correct the signs + if(poly_swapped) { + int p = P.degree(); + int q = Q.degree(); + for(int i=0;i<=q;i++) { + if((p-i)*(q-i) % 2 == 1) { + sres[q-i] = -sres[q-i]; + coP[q-i] = -coP[q-i]; + coQ[q-i] = -coQ[q-i]; + } + } + for(int i=0;i<=q;i++) { + // Swap coP and coQ: + CGAL::Polynomial help = coP[i]; + coP[i] = coQ[i]; + coQ[i] = help; + } + } + + // Now, reverse the entries + std::copy(coP.rbegin(),coP.rend(),coP_out); + std::copy(coQ.rbegin(),coQ.rend(),coQ_out); + return std::copy(sres.rbegin(),sres.rend(),sres_out); + + } + + // the general function for CGAL::Integral_domain_without_division_tag + template inline + OutputIterator + polynomial_subresultants_(CGAL::Polynomial A, + CGAL::Polynomial B, + OutputIterator out, + CGAL::Integral_domain_without_division_tag){ + + return bezout_polynomial_subresultants(A,B,out); + + } + + + // the specialization for CGAL::Unique_factorization_domain_tag + template inline + OutputIterator + polynomial_subresultants_(CGAL::Polynomial A, + CGAL::Polynomial B, + OutputIterator out, + CGAL::Unique_factorization_domain_tag){ + + return prs_polynomial_subresultants(A,B,out); + + } + + template inline + OutputIterator polynomial_subresultants_(CGAL::Polynomial A, + CGAL::Polynomial B, + OutputIterator out) { + typedef typename + CGAL::Algebraic_structure_traits::Algebraic_category + Algebraic_category; + return polynomial_subresultants_(A,B,out,Algebraic_category()); + } + + // the general function for CGAL::Integral_domain_without_division_tag + template inline + OutputIterator + principal_subresultants_(CGAL::Polynomial A, + CGAL::Polynomial B, + OutputIterator out, + CGAL::Integral_domain_without_division_tag){ + + return bezout_principal_subresultants(A,B,out); + + } + + // the specialization for CGAL::Unique_factorization_domain_tag + template inline + OutputIterator + principal_subresultants_(CGAL::Polynomial A, + CGAL::Polynomial B, + OutputIterator out, + CGAL::Unique_factorization_domain_tag){ + + return prs_principal_subresultants(A,B,out); + + } + + template inline + OutputIterator principal_subresultants_(CGAL::Polynomial A, + CGAL::Polynomial B, + OutputIterator out) { + typedef typename + CGAL::Algebraic_structure_traits::Algebraic_category + Algebraic_category; + return principal_subresultants_(A,B,out,Algebraic_category()); + } + + + /*! \relates CGAL::Polynomial + * \brief compute the polynomial subresultants of the polynomials + * \c A and \c B + * + * If \c n is the degree of A, the routine returns a sequence + * of length n+1, the (polynomial) subresultants of \c A and \c B. + * It starts with the resultant of \c A and \c B. + * The ith polynomial has degree + * at most i, and the last polynomial is \c A itself. + * + * The way the subresultants are computed depends on the Algebra_type. + * In general the subresultant will be computed by the function + * CGAL::bezout_polynomial_subresultants, but if possible the function + * CGAL::prs_polynomial_subresultants is used. + */ + template inline + OutputIterator polynomial_subresultants(CGAL::Polynomial A, CGAL::Polynomial B, + OutputIterator out) { + return CGAL::CGALi::polynomial_subresultants_(A, B, out); + } + + /*! \relates CGAL::Polynomial + * \brief compute the principal subresultants of the polynomials + * \c A and \c B + * + * If \c q is the degree of B, the routine returns a sequence + * of length q+1, the (principal) subresultants of \c A and \c B. + * It starts with the resultant of \c A and \c B, + * and ends with the leading coefficient of \c B. + * + * The way the subresultants are computed depends on the Algebra_type. + * In general the subresultant will be computed by the function + * CGAL::bezout_principal_subresultants, but if possible the function + * CGAL::prs_principal_subresultants is used. + */ + template inline + OutputIterator principal_subresultants(CGAL::Polynomial A, CGAL::Polynomial B, + OutputIterator out) { + return CGAL::CGALi::principal_subresultants_(A, B, out); + } + + } // namespace CGALi + +CGAL_END_NAMESPACE +#endif// CGAL_ACK_SUBRESULTANTS_H diff --git a/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp b/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp new file mode 100644 index 00000000000..7f05e5cff9c --- /dev/null +++ b/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp @@ -0,0 +1,347 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Michael Kerber +// +// ============================================================================ + +#include + +#include +#include +#include + +#include +#include + + +template void test_routine() { + + // All results are computed with the sres_stha.mpl worksheet + typedef typename ArithmeticKernel::Integer NT; + typedef CGAL::Polynomial PNT_1; + typedef CGAL::Polynomial PNT_2; + + { // test for the principal Sturm-Habicht sequence + // f=x^3*(x+1)*(x-1)*(x-2)^2 + PNT_1 f(NT(0),NT(0),NT(0),NT(-4),NT(4),NT(3),NT(-4),NT(1)); + std::vector stha; + CGAL::CGALi::principal_sturm_habicht_sequence(f,std::back_inserter(stha)); + + assert(stha.size()==8); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(0)); + assert(stha[2] == NT(0)); + assert(stha[3] == NT(864)); + assert(stha[4] == NT(324)); + assert(stha[5] == NT(54)); + assert(stha[6] == NT(7)); + assert(stha[7] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots + (stha.begin(),stha.end()) == 4); + + stha.clear(); + + CGAL::CGALi::prs_principal_sturm_habicht_sequence(f,std::back_inserter(stha)); + + assert(stha.size()==8); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(0)); + assert(stha[2] == NT(0)); + assert(stha[3] == NT(864)); + assert(stha[4] == NT(324)); + assert(stha[5] == NT(54)); + assert(stha[6] == NT(7)); + assert(stha[7] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots + (stha.begin(),stha.end()) == 4); + + stha.clear(); + + CGAL::CGALi::bezout_principal_sturm_habicht_sequence(f,std::back_inserter(stha)); + + assert(stha.size()==8); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(0)); + assert(stha[2] == NT(0)); + assert(stha[3] == NT(864)); + assert(stha[4] == NT(324)); + assert(stha[5] == NT(54)); + assert(stha[6] == NT(7)); + assert(stha[7] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots + (stha.begin(),stha.end()) == 4); + + stha.clear(); + + typename + CGAL::Polynomial_traits_d::Principal_sturm_habicht_sequence() + (f,std::back_inserter(stha)); + + assert(stha.size()==8); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(0)); + assert(stha[2] == NT(0)); + assert(stha[3] == NT(864)); + assert(stha[4] == NT(324)); + assert(stha[5] == NT(54)); + assert(stha[6] == NT(7)); + assert(stha[7] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots + (stha.begin(),stha.end()) == 4); + } + + { // Example for the defective case + // f:=x^8+x^2 + PNT_1 f(NT(0),NT(0),NT(1),NT(0),NT(0),NT(0),NT(0),NT(0),NT(1)); + + std::vector stha; + CGAL::CGALi::principal_sturm_habicht_sequence (f,std::back_inserter(stha)); + + assert(stha.size()==9); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(-93312)); + assert(stha[2] == NT(-62208)); + assert(stha[3] == NT(0)); + assert(stha[4] == NT(0)); + assert(stha[5] == NT(0)); + assert(stha[6] == NT(0)); + assert(stha[7] == NT(8)); + assert(stha[8] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots(stha.begin(),stha.end()) == 1); + + stha.clear(); + + CGAL::CGALi::prs_principal_sturm_habicht_sequence (f,std::back_inserter(stha)); + + assert(stha.size()==9); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(-93312)); + assert(stha[2] == NT(-62208)); + assert(stha[3] == NT(0)); + assert(stha[4] == NT(0)); + assert(stha[5] == NT(0)); + assert(stha[6] == NT(0)); + assert(stha[7] == NT(8)); + assert(stha[8] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots(stha.begin(),stha.end()) == 1); + + stha.clear(); + + CGAL::CGALi::bezout_principal_sturm_habicht_sequence (f,std::back_inserter(stha)); + + assert(stha.size()==9); + assert(stha[0] == NT(0)); + assert(stha[1] == NT(-93312)); + assert(stha[2] == NT(-62208)); + assert(stha[3] == NT(0)); + assert(stha[4] == NT(0)); + assert(stha[5] == NT(0)); + assert(stha[6] == NT(0)); + assert(stha[7] == NT(8)); + assert(stha[8] == NT(1)); + assert(CGAL::CGALi::stha_count_number_of_real_roots(stha.begin(),stha.end()) == 1); + + } + + + { // test for the complete Sturm-Habicht sequence + // f:=(2*x+1)^3*(x+1)*(x-1) + PNT_1 f(NT(-1),NT(-6),NT(-11),NT(-2),NT(12),NT(8)); + PNT_1 fx(f); + fx.diff(); + + std::vector stha; + + CGAL::CGALi::sturm_habicht_sequence(f,std::back_inserter(stha)); + + assert(stha.size()==6); + + assert(stha[5]==f); + assert(stha[4]==fx); + assert(stha[3]==PNT_1(NT(1024),NT(5568),NT(9984),NT(5888))); + assert(stha[2]==PNT_1(NT(55296),NT(221184),NT(221184))); + assert(stha[1].is_zero()); + assert(stha[0].is_zero()); + + stha.clear(); + + typename CGAL::Polynomial_traits_d::Sturm_habicht_sequence() + (f,std::back_inserter(stha)); + + assert(stha.size()==6); + + assert(stha[5]==f); + assert(stha[4]==fx); + assert(stha[3]==PNT_1(NT(1024),NT(5568),NT(9984),NT(5888))); + assert(stha[2]==PNT_1(NT(55296),NT(221184),NT(221184))); + assert(stha[1].is_zero()); + assert(stha[0].is_zero()); + + } + + { // test for cofactors + // f:=(2*x+1)^3*(x+1)*(x-1) + PNT_1 f(NT(-1),NT(-6),NT(-11),NT(-2),NT(12),NT(8)); + PNT_1 fx = CGAL::diff(f); + + std::vector stha_check,stha,co_f,co_fx; + + CGAL::CGALi::sturm_habicht_sequence(f,std::back_inserter(stha_check)); + CGAL::CGALi::sturm_habicht_sequence_with_cofactors(f, + std::back_inserter(stha), + std::back_inserter(co_f), + std::back_inserter(co_fx) + ); + + int n = static_cast(stha.size()); + assert(n == static_cast(stha_check.size())); + for(int i = 0; i < n; i++) { + assert(stha[i]==stha_check[i]); + assert(stha[i]==co_f[i]*f + co_fx[i]*fx); + } + + } + + { // Cofactors, defective case + // f:=x^8+x^2 + PNT_1 f(NT(0),NT(0),NT(1),NT(0),NT(0),NT(0),NT(0),NT(0),NT(1)); + PNT_1 fx = CGAL::diff(f); + + std::vector stha_check,stha,co_f,co_fx; + + + CGAL::CGALi::sturm_habicht_sequence(f,std::back_inserter(stha_check)); + CGAL::CGALi::sturm_habicht_sequence_with_cofactors(f, + std::back_inserter(stha), + std::back_inserter(co_f), + std::back_inserter(co_fx) + ); + + int n = static_cast(stha.size()); + assert(n == static_cast(stha_check.size())); + for(int i = 0; i < n; i++) { + assert(stha[i]==stha_check[i]); + assert(stha[i]==co_f[i]*f + co_fx[i]*fx); + } + + + } + + { // Test for bivariate polynomials + // f:=y^4 + 6*y^3 + 9*y^2 + (9*x^3 + 3*x^2 + 2*x + 5)*y + (2*x^2 + 5*x + 7) + std::stringstream ss("P[4(0,P[2(0,7)(1,5)(2,2)])(1,P[3(0,5)(1,2)(2,3)(3,9)])(2,P[0(0,9)])(3,P[0(0,6)])(4,P[0(0,1)])]"); + PNT_2 f; + ss >> f; + + { + std::vector stha,costha; + CGAL::CGALi::first_two_sturm_habicht_coefficients(f,std::back_inserter(stha), + std::back_inserter(costha)); + + assert(stha.size()==5); + assert(costha.size()==4); + + NT c1[13] = { NT(-81383),NT(-113448),NT(-123840),NT(-223936),NT(-272142), + NT(-14700),NT(419954),NT(156816),NT(-247131),NT(-498636), + NT(-275562),NT(-236196),NT(-177147) }; + PNT_1 p1(c1,c1+13); + assert(stha[0] == p1); + assert(stha[1] == PNT_1(NT(-828),NT(-1008),NT(-864),NT(-1728),NT(-1620), + NT(-1944),NT(-2916)) ); + assert(stha[2] == PNT_1(NT(36)) ); + assert(stha[3] == PNT_1(NT(4)) ); + assert(stha[4] == PNT_1(NT(1)) ); + + + + assert(costha[0] == PNT_1(NT(-2742),NT(-2592),NT(-1788),NT(-2100), + NT(-1638),NT(108),NT(1458)) ); + assert(costha[1] == PNT_1(NT(48),NT(-24),NT(-36),NT(-108)) ); + assert(costha[2] == PNT_1(NT(18)) ); + assert(costha[3] == PNT_1(NT(6)) ); + } + { + std::vector stha,costha; + CGAL::CGALi::prs_first_two_sturm_habicht_coefficients(f,std::back_inserter(stha), + std::back_inserter(costha)); + + assert(stha.size()==5); + assert(costha.size()==4); + + NT c1[13] = { NT(-81383),NT(-113448),NT(-123840),NT(-223936),NT(-272142), + NT(-14700),NT(419954),NT(156816),NT(-247131),NT(-498636), + NT(-275562),NT(-236196),NT(-177147) }; + PNT_1 p1(c1,c1+13); + assert(stha[0] == p1); + assert(stha[1] == PNT_1(NT(-828),NT(-1008),NT(-864),NT(-1728),NT(-1620), + NT(-1944),NT(-2916)) ); + assert(stha[2] == PNT_1(NT(36)) ); + assert(stha[3] == PNT_1(NT(4)) ); + assert(stha[4] == PNT_1(NT(1)) ); + + + + assert(costha[0] == PNT_1(NT(-2742),NT(-2592),NT(-1788),NT(-2100), + NT(-1638),NT(108),NT(1458)) ); + assert(costha[1] == PNT_1(NT(48),NT(-24),NT(-36),NT(-108)) ); + assert(costha[2] == PNT_1(NT(18)) ); + assert(costha[3] == PNT_1(NT(6)) ); + } + { + std::vector stha,costha; + CGAL::CGALi::bezout_first_two_sturm_habicht_coefficients + (f,std::back_inserter(stha), + std::back_inserter(costha)); + + assert(stha.size()==5); + assert(costha.size()==4); + + NT c1[13] = { NT(-81383),NT(-113448),NT(-123840),NT(-223936),NT(-272142), + NT(-14700),NT(419954),NT(156816),NT(-247131),NT(-498636), + NT(-275562),NT(-236196),NT(-177147) }; + PNT_1 p1(c1,c1+13); + assert(stha[0] == p1); + assert(stha[1] == PNT_1(NT(-828),NT(-1008),NT(-864),NT(-1728),NT(-1620), + NT(-1944),NT(-2916)) ); + assert(stha[2] == PNT_1(NT(36)) ); + assert(stha[3] == PNT_1(NT(4)) ); + assert(stha[4] == PNT_1(NT(1)) ); + + + + assert(costha[0] == PNT_1(NT(-2742),NT(-2592),NT(-1788),NT(-2100), + NT(-1638),NT(108),NT(1458)) ); + assert(costha[1] == PNT_1(NT(48),NT(-24),NT(-36),NT(-108)) ); + assert(costha[2] == PNT_1(NT(18)) ); + assert(costha[3] == PNT_1(NT(6)) ); + } + + + } + + + +} + +int main() { + +#ifdef CGAL_USE_LEDA + test_routine(); +#else + std::cout << "LEDA tests skipped" << std::endl; +#endif + +#ifdef CGAL_USE_CORE + test_routine(); +#else + std::cout << "CORE tests skipped" << std::endl; +#endif + return 0; +} diff --git a/Polynomial/test/Polynomial/subresultants.cpp b/Polynomial/test/Polynomial/subresultants.cpp new file mode 100644 index 00000000000..4b48e699a05 --- /dev/null +++ b/Polynomial/test/Polynomial/subresultants.cpp @@ -0,0 +1,472 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Michael Kerber +// +// ============================================================================ +#include +#include +#include + +template Poly_ from_string(const char* s) { + std::stringstream ss(s); + Poly_ f; + ss >> f; + return f; +} + +template +void test_routine() { + typedef ArithmeticKernel Arithmetic_kernel; + + typedef typename Arithmetic_kernel::Rational Rational; + typedef typename Arithmetic_kernel::Integer Integer; + + typedef CGAL::Polynomial Poly_int1; + typedef CGAL::Polynomial Poly_int2; + + typedef CGAL::Polynomial Poly_rat1; + { + //Example for the regular case: + Poly_int1 f(-5,-2,3,-6,-7,3,-2,4); + Poly_int1 g(0,-1,5,-7,5); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Integer(25)*g); + CGAL_assertion(sres[3]==Poly_int1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_int1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_int1(Integer(602925),Integer(657683))); + CGAL_assertion(sres[0]==Poly_int1(Integer(4474810))); + sres.clear(); + CGAL::CGALi::bezout_polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Integer(25)*g); + CGAL_assertion(sres[3]==Poly_int1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_int1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_int1(Integer(602925),Integer(657683))); + CGAL_assertion(sres[0]==Poly_int1(Integer(4474810))); + sres.clear(); + CGAL::CGALi::polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Integer(25)*g); + CGAL_assertion(sres[3]==Poly_int1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_int1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_int1(Integer(602925),Integer(657683))); + CGAL_assertion(sres[0]==Poly_int1(Integer(4474810))); + sres.clear(); + typename CGAL::Polynomial_traits_d::Polynomial_subresultants() + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Integer(25)*g); + CGAL_assertion(sres[3]==Poly_int1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_int1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_int1(Integer(602925),Integer(657683))); + CGAL_assertion(sres[0]==Poly_int1(Integer(4474810))); + std::vector psres; + CGAL::CGALi::prs_principal_subresultants(f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + psres.clear(); + CGAL::CGALi::bezout_principal_subresultants(f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + psres.clear(); + CGAL::CGALi::principal_subresultants(f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + psres.clear(); + typename CGAL::Polynomial_traits_d::Principal_subresultants() + (f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + } + + { + // Defective exampke + Poly_int1 f(1,-2,0,0,0,0,1); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants + (f,CGAL::diff(f),std::back_inserter(sres)); + CGAL_assertion(sres.size()==6); + CGAL_assertion(sres[5]==CGAL::diff(f)); + CGAL_assertion(sres[4]==Poly_int1(Integer(36),Integer(-60))); + CGAL_assertion(sres[3]==Poly_int1(0)); + CGAL_assertion(sres[2]==Poly_int1(0)); + CGAL_assertion(sres[1]==Poly_int1(Integer(-36000),Integer(60000))); + CGAL_assertion(sres[0]==Poly_int1(Integer(-153344))); + sres.clear(); + CGAL::CGALi::bezout_polynomial_subresultants + (f,CGAL::diff(f),std::back_inserter(sres)); + CGAL_assertion(sres.size()==6); + CGAL_assertion(sres[5]==CGAL::diff(f)); + CGAL_assertion(sres[4]==Poly_int1(Integer(36),Integer(-60))); + CGAL_assertion(sres[3]==Poly_int1(0)); + CGAL_assertion(sres[2]==Poly_int1(0)); + CGAL_assertion(sres[1]==Poly_int1(Integer(-36000),Integer(60000))); + CGAL_assertion(sres[0]==Poly_int1(Integer(-153344))); + std::vector psres; + CGAL::CGALi::prs_principal_subresultants(f,CGAL::diff(f), + std::back_inserter(psres)); + CGAL_assertion(psres.size()==6); + CGAL_assertion(psres[5]==6); + CGAL_assertion(psres[4]==0); + CGAL_assertion(psres[3]==0); + CGAL_assertion(psres[2]==0); + CGAL_assertion(psres[1]==60000); + CGAL_assertion(psres[0]==-153344); + psres.clear(); + CGAL::CGALi::bezout_principal_subresultants(f,CGAL::diff(f), + std::back_inserter(psres)); + CGAL_assertion(psres.size()==6); + CGAL_assertion(psres[5]==6); + CGAL_assertion(psres[4]==0); + CGAL_assertion(psres[3]==0); + CGAL_assertion(psres[2]==0); + CGAL_assertion(psres[1]==60000); + CGAL_assertion(psres[0]==-153344); + + + } + { + //Example for rational values: + Poly_rat1 f(-5,-2,3,-6,-7,3,-2,4); + Poly_rat1 g(0,-1,5,-7,5); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Rational(25)*g); + CGAL_assertion(sres[3]==Poly_rat1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_rat1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_rat1(Rational(602925),Rational(657683))); + CGAL_assertion(sres[0]==Poly_rat1(Rational(4474810))); + sres.clear(); + CGAL::CGALi::bezout_polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Rational(25)*g); + CGAL_assertion(sres[3]==Poly_rat1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_rat1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_rat1(Rational(602925),Rational(657683))); + CGAL_assertion(sres[0]==Poly_rat1(Rational(4474810))); + sres.clear(); + CGAL::CGALi::polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + CGAL_assertion(sres[4]==Rational(25)*g); + CGAL_assertion(sres[3]==Poly_rat1(-3125,-1768,4970,-9451)); + CGAL_assertion(sres[2]==Poly_rat1(206535,-262340,252423)); + CGAL_assertion(sres[1]==Poly_rat1(Rational(602925),Rational(657683))); + CGAL_assertion(sres[0]==Poly_rat1(Rational(4474810))); + std::vector psres; + CGAL::CGALi::prs_principal_subresultants(f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + psres.clear(); + CGAL::CGALi::bezout_principal_subresultants(f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + psres.clear(); + CGAL::CGALi::principal_subresultants(f,g,std::back_inserter(psres)); + CGAL_assertion(psres.size()==5); + CGAL_assertion(psres[4]==125); + CGAL_assertion(psres[3]==-9451); + CGAL_assertion(psres[2]==252423); + CGAL_assertion(psres[1]==657683); + CGAL_assertion(psres[0]==4474810); + psres.clear(); + } + + { + // Another test included on 10 Jan 2007 - situation comes from arrangement + // computation of curves: + + // f:=y^4 + (x)*y^2 + (2*x^4 + (-1)*x^3); + Poly_int2 f=from_string("P[4(0,P[4(3,-1)(4,2)])(2,P[1(1,1)])(4,P[0(0,1)])]"); + // g:=g:=eval(10000*f,[y=2*y,x=x-1/10]); + Poly_int2 g=from_string("P[4(0,P[4(0,12)(1,-380)(2,4200)(3,-18000)(4,20000)])(2,P[1(0,-4000)(1,40000)])(4,P[0(0,160000)])]"); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + // Computed with MAPLE + CGAL_assertion(sres[4]==g); + CGAL_assertion(sres[3]==from_string("P[2(0,P[4(0,12)(1,-380)(2,4200)(3,142000)(4,-300000)])(2,P[1(0,-4000)(1,-120000)])]")); + CGAL_assertion(sres[2]==from_string("P[2(0,P[5(0,-48000)(1,80000)(2,28800000)(3,-1072000000)(4,-15840000000)(5,36000000000)])(2,P[2(0,16000000)(1,960000000)(2,14400000000)])]")); + CGAL_assertion(sres[1]==from_string("P[0(0,P[9(0,576000)(1,172800000)(2,5326400000)(3,-158512000000)(4,-5164000000000)(5,24705600000000)(6,615472000000000)(7,912480000000000)(8,-9864000000000000)(9,10800000000000000)])]")); + CGAL_assertion(sres[0]==from_string("P[0(0,P[16(0,20736)(1,11197440)(2,1559232000)(3,5760000)(4,-3426163040000)(5,-9736288000000)(6,2377866144000000)(7,-1780931200000000)(8,-427278798400000000)(9,-507616640000000000)(10,31454608000000000000)(11,83909222400000000000)(12,-697197584000000000000)(13,-919113600000000000000)(14,9138960000000000000000)(15,-15336000000000000000000)(16,8100000000000000000000)])]")); + sres.clear(); + CGAL::CGALi::bezout_polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + // Computed with MAPLE + CGAL_assertion(sres[4]==g); + CGAL_assertion(sres[3]==from_string("P[2(0,P[4(0,12)(1,-380)(2,4200)(3,142000)(4,-300000)])(2,P[1(0,-4000)(1,-120000)])]")); + CGAL_assertion(sres[2]==from_string("P[2(0,P[5(0,-48000)(1,80000)(2,28800000)(3,-1072000000)(4,-15840000000)(5,36000000000)])(2,P[2(0,16000000)(1,960000000)(2,14400000000)])]")); + CGAL_assertion(sres[1]==from_string("P[0(0,P[9(0,576000)(1,172800000)(2,5326400000)(3,-158512000000)(4,-5164000000000)(5,24705600000000)(6,615472000000000)(7,912480000000000)(8,-9864000000000000)(9,10800000000000000)])]")); + CGAL_assertion(sres[0]==from_string("P[0(0,P[16(0,20736)(1,11197440)(2,1559232000)(3,5760000)(4,-3426163040000)(5,-9736288000000)(6,2377866144000000)(7,-1780931200000000)(8,-427278798400000000)(9,-507616640000000000)(10,31454608000000000000)(11,83909222400000000000)(12,-697197584000000000000)(13,-919113600000000000000)(14,9138960000000000000000)(15,-15336000000000000000000)(16,8100000000000000000000)])]")); + sres.clear(); + CGAL::CGALi::polynomial_subresultants(f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==5); + // Computed with MAPLE + CGAL_assertion(sres[4]==g); + CGAL_assertion(sres[3]==from_string("P[2(0,P[4(0,12)(1,-380)(2,4200)(3,142000)(4,-300000)])(2,P[1(0,-4000)(1,-120000)])]")); + CGAL_assertion(sres[2]==from_string("P[2(0,P[5(0,-48000)(1,80000)(2,28800000)(3,-1072000000)(4,-15840000000)(5,36000000000)])(2,P[2(0,16000000)(1,960000000)(2,14400000000)])]")); + CGAL_assertion(sres[1]==from_string("P[0(0,P[9(0,576000)(1,172800000)(2,5326400000)(3,-158512000000)(4,-5164000000000)(5,24705600000000)(6,615472000000000)(7,912480000000000)(8,-9864000000000000)(9,10800000000000000)])]")); + CGAL_assertion(sres[0]==from_string("P[0(0,P[16(0,20736)(1,11197440)(2,1559232000)(3,5760000)(4,-3426163040000)(5,-9736288000000)(6,2377866144000000)(7,-1780931200000000)(8,-427278798400000000)(9,-507616640000000000)(10,31454608000000000000)(11,83909222400000000000)(12,-697197584000000000000)(13,-919113600000000000000)(14,9138960000000000000000)(15,-15336000000000000000000)(16,8100000000000000000000)])]")); + } + + { // Test for f.degree() < g.degree() + Poly_int1 f = from_string("P[10(2,-4)(3,-73)(6,97)(7,-62)(10,-56)]"); + Poly_int1 g = from_string("P[15(5,-75)(6,-17)(8,71)(11,-44)(14,80)(15,-82)]"); + + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==11); + CGAL_assertion(sres[10]==from_string("P[10(2,-39337984)(3,-717918208)(6,953946112)(7,-609738752)(10,-550731776)]")); + CGAL_assertion(sres[9]==from_string("P[9(2,-305262755840)(3,-4966105776128)(4,10840151891968)(5,-5962969628672)(6,6702091010048)(7,-22436989575168)(8,19712814514176)(9,-3099911815168)]")); + CGAL_assertion(sres[8]==from_string("P[8(2,9680208765108224)(3,156728928112275456)(4,-360058057701515264)(5,152421565222248448)(6,-176106510157656064)(7,746064242175737856)(8,-579306262853885952)]")); + CGAL_assertion(sres[7]==from_string("P[7(2,-1486764415342623952896)(3,-26708453688140199103488)(4,8054371431842323176960)(5,3491478546579489027072)(6,38677404054075120663552)(7,-43611415672712033527296)]")); + CGAL_assertion(sres[6]==from_string("P[6(2,9979660953345863413381568)(3,193906324691706110661441648)(4,213219021894852075965695648)(5,362883385829743140658019008)(6,-93331232763027671680558576)]")); + CGAL_assertion(sres[5]==from_string("P[5(2,-70907443741718133356344565824)(3,-1389117178389007568663227054224)(4,-1747565962047585993090805107776)(5,-2771071359404167222865326731904)]")); + CGAL_assertion(sres[4]==from_string("P[4(2,-715876861001937085090535662991104)(3,-13330409138559036348802401660667328)(4,-5258395198219675575175450753910144)]")); + CGAL_assertion(sres[3]==from_string("P[3(2,-2842398393513781389004679512270277696)(3,-51817649494648561412362002061296745104)]")); + CGAL_assertion(sres[2]==from_string("P[2(2,-46088453959546967380998470482217286656)]")); + CGAL_assertion(sres[1]==Poly_int1(Integer(0))); + CGAL_assertion(sres[0]==Poly_int1(Integer(0))); + sres.clear(); + + CGAL::CGALi::bezout_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==11); + CGAL_assertion(sres[10]==from_string("P[10(2,-39337984)(3,-717918208)(6,953946112)(7,-609738752)(10,-550731776)]")); + CGAL_assertion(sres[9]==from_string("P[9(2,-305262755840)(3,-4966105776128)(4,10840151891968)(5,-5962969628672)(6,6702091010048)(7,-22436989575168)(8,19712814514176)(9,-3099911815168)]")); + CGAL_assertion(sres[8]==from_string("P[8(2,9680208765108224)(3,156728928112275456)(4,-360058057701515264)(5,152421565222248448)(6,-176106510157656064)(7,746064242175737856)(8,-579306262853885952)]")); + CGAL_assertion(sres[7]==from_string("P[7(2,-1486764415342623952896)(3,-26708453688140199103488)(4,8054371431842323176960)(5,3491478546579489027072)(6,38677404054075120663552)(7,-43611415672712033527296)]")); + CGAL_assertion(sres[6]==from_string("P[6(2,9979660953345863413381568)(3,193906324691706110661441648)(4,213219021894852075965695648)(5,362883385829743140658019008)(6,-93331232763027671680558576)]")); + CGAL_assertion(sres[5]==from_string("P[5(2,-70907443741718133356344565824)(3,-1389117178389007568663227054224)(4,-1747565962047585993090805107776)(5,-2771071359404167222865326731904)]")); + CGAL_assertion(sres[4]==from_string("P[4(2,-715876861001937085090535662991104)(3,-13330409138559036348802401660667328)(4,-5258395198219675575175450753910144)]")); + CGAL_assertion(sres[3]==from_string("P[3(2,-2842398393513781389004679512270277696)(3,-51817649494648561412362002061296745104)]")); + CGAL_assertion(sres[2]==from_string("P[2(2,-46088453959546967380998470482217286656)]")); + CGAL_assertion(sres[1]==Poly_int1(Integer(0))); + CGAL_assertion(sres[0]==Poly_int1(Integer(0))); + sres.clear(); + } + { // Another example where sign switches are necessary + Poly_int1 f = from_string("P[3(0,42)(1,-40)(2,-7)(3,-10)]"); + Poly_int1 g = from_string("P[5(0,74)(1,6)(2,-92)(3,75)(4,23)(5,-50)]"); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==4); + CGAL_assertion(sres[3]==from_string("P[3(0,-420)(1,400)(2,70)(3,100)]")); + CGAL_assertion(sres[2]==from_string("P[2(0,-1058480)(1,688000)(2,698080)]")); + CGAL_assertion(sres[1]==from_string("P[1(0,-22577275200)(1,28253151360)]")); + CGAL_assertion(sres[0]==from_string("P[0(0,-103066942158720)]")); + sres.clear(); + CGAL::CGALi::bezout_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==4); + CGAL_assertion(sres[3]==from_string("P[3(0,-420)(1,400)(2,70)(3,100)]")); + CGAL_assertion(sres[2]==from_string("P[2(0,-1058480)(1,688000)(2,698080)]")); + CGAL_assertion(sres[1]==from_string("P[1(0,-22577275200)(1,28253151360)]")); + CGAL_assertion(sres[0]==from_string("P[0(0,-103066942158720)]")); + } + { // Test cases with constant polynomials + { + Poly_int1 f(Integer(8)); + Poly_int1 g=from_string("P[3(0,43)(1,40)(2,-19)(3,2)]"); + std::vector sres; + + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0].degree()==0); + CGAL_assertion(sres[0][0]==Integer(8*8*8)); + sres.clear(); + + CGAL::CGALi::bezout_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0].degree()==0); + CGAL_assertion(sres[0][0]==Integer(8*8*8)); + std::vector psres; + + CGAL::CGALi::prs_principal_subresultants + (g,f,std::back_inserter(psres)); + CGAL_assertion(psres.size()==1); + CGAL_assertion(psres[0]==Integer(8*8*8)); + psres.clear(); + + CGAL::CGALi::bezout_principal_subresultants + (g,f,std::back_inserter(psres)); + CGAL_assertion(psres.size()==1); + CGAL_assertion(psres[0]==Integer(8*8*8)); + } + { + Poly_int1 f(Integer(0)); + Poly_int1 g=from_string("P[1(0,1)(1,1)]"); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0].is_zero()); + sres.clear(); + + CGAL::CGALi::bezout_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0].is_zero()); + } + { + Poly_int1 f(Integer(7)); + Poly_int1 g(Integer(-12)); + + std::vector sres; + CGAL::CGALi::prs_principal_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0]==Integer(1)); + sres.clear(); + + CGAL::CGALi::bezout_principal_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0]==Integer(1)); + } + { + Poly_int2 f=from_string("P[0(0,P[0(0,122)])]"); + Poly_int2 g=from_string("P[0(0,P[0(0,0)])]"); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0].is_zero()); + sres.clear(); + + CGAL::CGALi::bezout_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==1); + CGAL_assertion(sres[0].is_zero()); + } + + } + { // Test for two linear polynomials: + Poly_int1 f = from_string("P[1(0,-5)(1,6)]"); + Poly_int1 g = from_string("P[1(0,-3)(1,4)]"); + std::vector sres; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==2); + CGAL_assertion(sres[1]==g); + CGAL_assertion(sres[0]==from_string("P[0(0,2)]")); + sres.clear(); + CGAL::CGALi::bezout_polynomial_subresultants + (f,g,std::back_inserter(sres)); + CGAL_assertion(sres.size()==2); + CGAL_assertion(sres[1]==g); + CGAL_assertion(sres[0]==from_string("P[0(0,2)]")); + } + { // Test for cofactors, univariate + Poly_int1 f = from_string("P[6(0,246)(1,100)(2,197)(3,14)(4,136)(5,191)(6,207)]"); + Poly_int1 g = from_string("P[5(0,100)(1,394)(2,42)(3,544)(4,955)(5,1242)]"); + std::vector sres_check,sres,coP,coQ; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres_check)); + CGAL::CGALi::prs_subresultants_with_cofactors + (f,g, + std::back_inserter(sres), + std::back_inserter(coP), + std::back_inserter(coQ)); + CGAL_assertion(sres.size()==sres_check.size()); + CGAL_assertion(sres.size()==coP.size()); + CGAL_assertion(sres.size()==coQ.size()); + for(int i=0;i < static_cast(sres.size()); i++) { + CGAL_assertion(sres[i]==sres_check[i]); + CGAL_assertion(sres[i]==coP[i]*f + coQ[i]*g); + } + } + { // Test for trivariate + typedef CGAL::Polynomial Poly_int3; + Poly_int3 f = from_string("P[6(0,P[6(0,P[6(2,3)(4,-3)(6,1)])(1,P[1(1,-2)])(2,P[4(0,3)(2,-5)(4,3)])(4,P[2(0,-3)(2,3)])(6,P[0(0,1)])])(3,P[1(0,P[0(0,2)])(1,P[1(1,-2)])])(6,P[0(0,P[0(0,1)])])]"); + Poly_int3 g = CGAL::diff(f); + std::vector sres_check,sres,coP,coQ; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres_check)); + CGAL::CGALi::prs_subresultants_with_cofactors + (f,g, + std::back_inserter(sres), + std::back_inserter(coP), + std::back_inserter(coQ)); + CGAL_assertion(sres.size()==sres_check.size()); + CGAL_assertion(sres.size()==coP.size()); + CGAL_assertion(sres.size()==coQ.size()); + for(int i=0;i < static_cast(sres.size()); i++) { + CGAL_assertion(sres[i]==sres_check[i]); + CGAL_assertion(sres[i]==coP[i]*f + coQ[i]*g); + } + } + { // Test for cofactors, bivariate + Poly_int2 f = from_string("P[7(0,P[7(0,72)(1,113)(2,238)(3,75)(4,77)(5,149)(6,34)(7,75)])(1,P[6(0,113)(1,69)(2,94)(3,171)(4,148)(5,103)(6,42)])(2,P[5(0,103)(1,233)(2,16)(3,131)(4,2)(5,156)])(3,P[4(0,249)(1,194)(2,156)(3,81)(4,176)])(4,P[3(0,39)(1,140)(2,134)(3,190)])(5,P[2(0,117)(1,249)(2,158)])(6,P[1(0,109)(1,46)])(7,P[0(0,236)])]"); + Poly_int2 g = from_string("P[6(0,P[6(0,113)(1,69)(2,94)(3,171)(4,148)(5,103)(6,42)])(1,P[5(0,206)(1,466)(2,32)(3,262)(4,4)(5,312)])(2,P[4(0,747)(1,582)(2,468)(3,243)(4,528)])(3,P[3(0,156)(1,560)(2,536)(3,760)])(4,P[2(0,585)(1,1245)(2,790)])(5,P[1(0,654)(1,276)])(6,P[0(0,1652)])]"); + std::vector sres_check,sres,coP,coQ; + CGAL::CGALi::prs_polynomial_subresultants + (f,g,std::back_inserter(sres_check)); + CGAL::CGALi::prs_subresultants_with_cofactors + (f,g, + std::back_inserter(sres), + std::back_inserter(coP), + std::back_inserter(coQ)); + CGAL_assertion(sres.size()==sres_check.size()); + CGAL_assertion(sres.size()==coP.size()); + CGAL_assertion(sres.size()==coQ.size()); + for(int i=0;i < static_cast(sres.size()); i++) { + CGAL_assertion(sres[i]==sres_check[i]); + CGAL_assertion(sres[i]==coP[i]*f + coQ[i]*g); + } + } + return; +} + + +int main(){ + +#ifdef CGAL_USE_LEDA + test_routine(); +#else + std::cerr << "LEDA tests skipped!" << std::endl; +#endif +#ifdef CGAL_USE_CORE + test_routine(); +#else + std::cerr << "CORE tests skipped!" << std::endl; +#endif + return 0; +}