songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Moved code to compute subresultants and sturm-habicht-sequences into Polynomial-package

This commit is contained in:
Michael Kerber 2008-06-27 10:03:35 +00:00
parent 6794adacb8
commit 5af98f25dc
7 changed files with 1807 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitattributes vendored
View File

@ -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

3
Polynomial/include/CGAL/Polynomial/bezout_matrix.h
View File

@ -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

27
Polynomial/include/CGAL/Polynomial/fwd.h
View File

@ -73,8 +73,35 @@ inline int square_free_factorization_for_regular_polynomial(const Polynomial<Coe
template <class NT> inline bool may_have_multiple_factor( const Polynomial<NT>&);
template <class NT> inline bool may_have_common_factor(const Polynomial<NT>&,const Polynomial<NT>&);
template< class Coeff >
struct Simple_matrix;
template<class NT>
CGALi::Simple_matrix<NT> polynomial_subresultant_matrix(
CGAL::Polynomial<NT> f,
CGAL::Polynomial<NT> g,
int d=0);
template <typename OutputIterator, typename NT> inline
OutputIterator polynomial_subresultants(CGAL::Polynomial<NT> A, CGAL::Polynomial<NT> B,
OutputIterator out);
template <typename OutputIterator, typename NT> inline
OutputIterator principal_subresultants(CGAL::Polynomial<NT> A, CGAL::Polynomial<NT> B,
OutputIterator out);
template <typename OutputIterator, typename NT > inline
OutputIterator principal_sturm_habicht_sequence(CGAL::Polynomial<NT> A,
OutputIterator out);
template<typename OutputIterator, typename NT> OutputIterator
sturm_habicht_sequence(CGAL::Polynomial<NT> P, OutputIterator out);
} // namespace CGALi
} // namespace CGAL

396
Polynomial/include/CGAL/Polynomial/sturm_habicht_sequence.h Normal file
View File

@ -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 <mkerber@mpi-inf.mpg.de>
//
// ============================================================================
#ifndef CGAL_ACK_STURM_HABICHT
#define CGAL_ACK_STURM_HABICHT 1
#include <vector>
#include <algorithm>
#include <CGAL/Polynomial/bezout_matrix.h>
#include <CGAL/Polynomial/subresultants.h>
CGAL_BEGIN_NAMESPACE
namespace CGALi {
/*!
* \brief compute the leading coefficients of the Sturm-Habicht sequence of
* the polynomial <I>P</I>
*
* The principal Sturm-Habicht sequence is obtained by computing the scalar
* subresultant sequence of <I>P</I> and its derivative, extended
* by <I>P</I> and <I>P'</I> 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<typename NT,typename OutputIterator>
OutputIterator prs_principal_sturm_habicht_sequence(CGAL::Polynomial<NT> P,
OutputIterator out) {
std::vector<CGAL::Polynomial<NT> > stha;
CGAL::CGALi::sturm_habicht_sequence(P,std::back_inserter(stha));
for(int i=0; i<static_cast<int>(stha.size()); i++) {
int d = stha[i].degree();
CGAL_assertion(d<=i);
if(d<i) {
*out++ = NT(0);
} else {
*out++ = stha[i][i];
}
}
return out;
}
/*!
* \brief compute the leading coefficients of the Sturm-Habicht sequence of
* the polynomial <I>P</I>
*
* The principal Sturm-Habicht sequence is obtained by computing the scalar
* subresultant sequence of <I>P</I> and its derivative, extended
* by <I>P</I> and <I>P'</I> 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<typename NT,typename OutputIterator>
OutputIterator bezout_principal_sturm_habicht_sequence
(CGAL::Polynomial<NT> P,
OutputIterator out) {
CGAL::Polynomial<NT> Px(P);
Px.diff();
CGAL::CGALi::Simple_matrix<NT> M
= CGAL::CGALi::polynomial_subresultant_matrix(P,Px,1);
int n = static_cast<int>(M.row_dimension());
for(int i=0; i<n; i++) {
if((n-1-i)%4==0 || (n-1-i)%4==1) {
*out++ = -M[n-1-i][n-1-i];
} else {
*out++ = M[n-1-i][n-1-i];
}
}
*out++=Px.lcoeff();
*out++=P.lcoeff();
return out;
}
/*! \ingroup NiX_resultant_matrix
* \brief compute the principal and coprincipal Sturm-Habicht sequence
*/
template<typename NT,typename OutputIterator1,typename OutputIterator2>
void prs_first_two_sturm_habicht_coefficients(CGAL::Polynomial<NT> P,
OutputIterator1 pstha,
OutputIterator2 copstha) {
std::vector<CGAL::Polynomial<NT> > stha;
int n = P.degree();
sturm_habicht_sequence(P,std::back_inserter(stha));
CGAL_assertion(static_cast<int>(stha.size())==n+1);
for(int i=0;i<=n;i++) {
int d = stha[i].degree();
CGAL_assertion(d<=i);
if(d<i) {
*pstha++ = NT(0);
} else {
*pstha++ = stha[i][i];
}
}
for(int i=1;i<=n;i++) {
int d = stha[i].degree();
CGAL_assertion(d<=i);
if(d<i-1) {
*copstha++ = NT(0);
} else {
*copstha++ = stha[i][i-1];
}
}
return;
}
/*! \brief compute the principal and coprincipal Sturm-Habicht sequence
*/
template<typename NT,typename OutputIterator1,typename OutputIterator2>
void bezout_first_two_sturm_habicht_coefficients(CGAL::Polynomial<NT> P,
OutputIterator1 pstha,
OutputIterator2 copstha) {
CGAL::Polynomial<NT> Px(P);
Px.diff();
CGAL::CGALi::Simple_matrix<NT> M
= CGAL::CGALi::polynomial_subresultant_matrix(P,Px,2);
int n = static_cast<int>(M.row_dimension());
for(int i=0; i<n; i++) {
if((n-1-i)%4==0 || (n-1-i)%4==1) {
*pstha++ = -M[n-1-i][n-1-i];
} else {
*pstha++ = M[n-1-i][n-1-i];
}
}
*pstha++ = Px.lcoeff();
*pstha++ = P.lcoeff();
for(int i=1; i<n; i++) {
if(n-i-1%4==0 || n-i-1%4==1) {
*copstha++ = -M[n-i-1][n-i];
} else {
*copstha++ = M[n-1-i][n-i];
}
}
*copstha++ = Px[Px.degree()-1];
*copstha++ = P[P.degree()-1];
}
// the general function for CGAL::Integral_domain_without_division_tag
template <typename OutputIterator1,
typename OutputIterator2,
typename NT> inline
void
first_two_sturm_habicht_coefficients_
(CGAL::Polynomial<NT> 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 <typename OutputIterator1,
typename OutputIterator2,
typename NT> inline
void
first_two_sturm_habicht_coefficients_
(CGAL::Polynomial<NT> A,
OutputIterator1 pstha,
OutputIterator2 copstha,
CGAL::Unique_factorization_domain_tag) {
return prs_first_two_sturm_habicht_coefficients(A,pstha,copstha);
}
template <typename OutputIterator1,
typename OutputIterator2,
typename NT > inline
void
first_two_sturm_habicht_coefficients_(CGAL::Polynomial<NT> A,
OutputIterator1 pstha,
OutputIterator2 copstha) {
typedef typename
CGAL::Algebraic_structure_traits<NT>::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 <typename OutputIterator, typename NT> inline
OutputIterator
principal_sturm_habicht_sequence_
(CGAL::Polynomial<NT> 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 <typename OutputIterator, typename NT> inline
OutputIterator
principal_sturm_habicht_sequence_
(CGAL::Polynomial<NT> A,
OutputIterator out,
CGAL::Unique_factorization_domain_tag) {
return prs_principal_sturm_habicht_sequence(A,out);
}
template <typename OutputIterator, typename NT > inline
OutputIterator principal_sturm_habicht_sequence_(CGAL::Polynomial<NT> A,
OutputIterator out) {
typedef typename
CGAL::Algebraic_structure_traits<NT>::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 <typename OutputIterator, typename NT> inline
OutputIterator
principal_sturm_habicht_sequence(CGAL::Polynomial<NT> 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 <typename OutputIterator1,
typename OutputIterator2,
typename NT> inline
void first_two_sturm_habicht_coefficients(CGAL::Polynomial<NT> 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<typename OutputIterator, typename NT> OutputIterator
sturm_habicht_sequence(CGAL::Polynomial<NT> P, OutputIterator out) {
int p = P.degree();
CGAL::Polynomial<NT> P_x(P);
P_x.diff();
std::vector<CGAL::Polynomial<NT> > stha;
CGAL::CGALi::polynomial_subresultants(P,P_x,std::back_inserter(stha));
stha.push_back(P);
CGAL_assertion(static_cast<int>(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<typename OutputIterator1,
typename OutputIterator2,
typename OutputIterator3,
typename NT> OutputIterator1
sturm_habicht_sequence_with_cofactors(CGAL::Polynomial<NT> P,
OutputIterator1 out_stha,
OutputIterator2 out_f,
OutputIterator3 out_fx) {
int p = P.degree();
CGAL::Polynomial<NT> P_x(P);
P_x.diff();
std::vector<CGAL::Polynomial<NT> > 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<NT>(1));
co_fx.push_back(CGAL::Polynomial<NT>(0));
CGAL_assertion(static_cast<int>(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<typename NT,typename InputIterator>
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

560
Polynomial/include/CGAL/Polynomial/subresultants.h Normal file
View File

@ -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 <mkerber@mpi-inf.mpg.de>
//
// ============================================================================
#ifndef CGAL_ACK_SUBRESULTANTS_H
#define CGAL_ACK_SUBRESULTANTS_H
#include <list>
#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial/bezout_matrix.h>
CGAL_BEGIN_NAMESPACE
namespace CGALi {
// Intern function needed for Ducos algorithm
template<typename NT> void lazard_optimization(NT y,
double n,
CGAL::Polynomial<NT> B,
CGAL::Polynomial<NT>& 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<typename NT>
void lickteig_roy_optimization(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> B,
CGAL::Polynomial<NT> C,
NT s,
CGAL::Polynomial<NT>& D) {
typedef CGAL::Polynomial<NT> Poly; // for convenience
int d = A.degree(), e = B.degree();
CGAL_precondition(d>=e);
std::vector<Poly> H(d+1);
std::list<NT> initial;
initial.push_front(C.lcoeff());
for(int i=0;i<e;i++) {
H[i] = Poly(initial.begin(),initial.end());
initial.push_front(NT(0));
}
H[e]=Poly(initial.begin(),initial.end())-C;
CGAL_assertion(H[e].degree()<e);
initial.clear();
std::copy(H[e].begin(),H[e].end(),std::back_inserter(initial));
initial.push_front(NT(0));
for(int i=e+1;i<d;i++) {
H[i]=Poly(initial.begin(),initial.end());
NT h_i_e=H[i].degree()>=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<d;i++) {
D+=A[i]*H[i];
}
D/=A.lcoeff();
NT Hde = H[d].degree()>=e ? H[d][e] : NT(0);
D=(B.lcoeff()*(H[d]+D)-Hde*B)/s;
if((d-e)%2==0) {
D=-D;
}
return;
}
template<typename NT> NT resultant_for_constant_polynomial
(CGAL::Polynomial<NT> P, CGAL::Polynomial<NT> 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. <i>Journal of Pure
* and Applied Algebra</i> <b>145</b> (2000) 149--163
*/
template <typename NT,typename OutputIterator> inline
OutputIterator prs_polynomial_subresultants(CGAL::Polynomial<NT> P,
CGAL::Polynomial<NT> 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<NT> zero_pol(NT(0));
std::vector<CGAL::Polynomial<NT> > 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<NT> 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<NT>::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;i<d;i++) {
sres.push_back(zero_pol);
}
break;
}
sres.push_back(B);
delta=d-e;
if(delta>1) {
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<delta-2;i++) {
sres.push_back(zero_pol);
}
sres.push_back(C);
}
else {
C=B;
}
if(e==0) {
break;
}
CGAL::CGALi::lickteig_roy_optimization(A,B,C,s,D);
B=D;
//CGAL::Polynomial<NT>::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<int>(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 <typename NT,typename OutputIterator> inline
OutputIterator bezout_polynomial_subresultants(CGAL::Polynomial<NT> P,
CGAL::Polynomial<NT> 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<NT> Matrix;
Matrix M = CGAL::CGALi::polynomial_subresultant_matrix(P,Q);
int r = static_cast<int>(M.row_dimension());
for(int i = 0;i < r; i++) {
std::vector<NT> 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<NT>(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 <typename NT,typename OutputIterator> inline
OutputIterator prs_principal_subresultants(CGAL::Polynomial<NT> P,
CGAL::Polynomial<NT> Q,
OutputIterator out) {
std::vector<CGAL::Polynomial<NT> > sres;
int q = std::min(Q.degree(),P.degree());
CGAL::CGALi::prs_polynomial_subresultants(P,Q,std::back_inserter(sres));
CGAL_assertion(static_cast<int>(sres.size()) == q+1);
for(int i=0; i <= q; i++) {
int d = sres[i].degree();
CGAL_assertion(d<=i);
if(d<i) {
*out++ = NT(0);
} else {
*out++ = sres[i][i];
}
}
return out;
}
/*!
* \brief Compute the sequence of principal subresultants
* with minors of the Bezout matrix
*
*/
template <typename NT,typename OutputIterator> inline
OutputIterator bezout_principal_subresultants(CGAL::Polynomial<NT> P,
CGAL::Polynomial<NT> 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<NT> Matrix;
Matrix M = CGAL::CGALi::polynomial_subresultant_matrix(P,Q,1);
int r = static_cast<int>(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<typename NT,
typename OutputIterator1,
typename OutputIterator2,
typename OutputIterator3>
OutputIterator1 prs_subresultants_with_cofactors(CGAL::Polynomial<NT> P,
CGAL::Polynomial<NT> 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<NT> zero_pol(NT(0));
std::vector<CGAL::Polynomial<NT> > 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<NT> A,B,C,D,Quo, coPA, coPB, coQA, coQB, coPC, coQC;
NT s,m;
int delta,d,e;
coPA = CGAL::Polynomial<NT>(NT(0));
coQA = CGAL::Polynomial<NT>(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<NT>::pseudo_division(P, -Q, Quo, B, m);
coPB = CGAL::Polynomial<NT>(m);
coQB = Quo;
while(true) {
d=A.degree();
e=B.degree();
if(B.is_zero()) {
for(int i=0;i<d;i++) {
sres.push_back(zero_pol);
coP.push_back(zero_pol);
coQ.push_back(zero_pol);
}
break;
}
sres.push_back(B);
coP.push_back(coPB);
coQ.push_back(coQB);
delta=d-e;
if(delta>1) {
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<delta-2;i++) {
sres.push_back(zero_pol);
coP.push_back(zero_pol);
coQ.push_back(zero_pol);
}
sres.push_back(C);
coP.push_back(coPC);
coQ.push_back(coQC);
}
else {
C=B;
coPC = coPB;
coQC = coQB;
}
if(e==0) {
break;
}
NT denominator = CGAL::ipower(s,delta)*A.lcoeff();
CGAL::Polynomial<NT>::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<int>(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<NT> 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 <typename OutputIterator, typename NT> inline
OutputIterator
polynomial_subresultants_(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> 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 <typename OutputIterator, typename NT> inline
OutputIterator
polynomial_subresultants_(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> B,
OutputIterator out,
CGAL::Unique_factorization_domain_tag){
return prs_polynomial_subresultants(A,B,out);
}
template <typename OutputIterator, typename NT > inline
OutputIterator polynomial_subresultants_(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> B,
OutputIterator out) {
typedef typename
CGAL::Algebraic_structure_traits<NT>::Algebraic_category
Algebraic_category;
return polynomial_subresultants_(A,B,out,Algebraic_category());
}
// the general function for CGAL::Integral_domain_without_division_tag
template <typename OutputIterator, typename NT> inline
OutputIterator
principal_subresultants_(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> 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 <typename OutputIterator, typename NT> inline
OutputIterator
principal_subresultants_(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> B,
OutputIterator out,
CGAL::Unique_factorization_domain_tag){
return prs_principal_subresultants(A,B,out);
}
template <typename OutputIterator, typename NT > inline
OutputIterator principal_subresultants_(CGAL::Polynomial<NT> A,
CGAL::Polynomial<NT> B,
OutputIterator out) {
typedef typename
CGAL::Algebraic_structure_traits<NT>::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 <tt>i</tt>th 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 <typename OutputIterator, typename NT> inline
OutputIterator polynomial_subresultants(CGAL::Polynomial<NT> A, CGAL::Polynomial<NT> 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 <typename OutputIterator, typename NT> inline
OutputIterator principal_subresultants(CGAL::Polynomial<NT> A, CGAL::Polynomial<NT> B,
OutputIterator out) {
return CGAL::CGALi::principal_subresultants_(A, B, out);
}
} // namespace CGALi
CGAL_END_NAMESPACE
#endif// CGAL_ACK_SUBRESULTANTS_H

347
Polynomial/test/Polynomial/sturm_habicht_sequence.cpp Normal file
View File

@ -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 <mkerber@mpi-inf.mpg.de>
//
// ============================================================================
#include <CGAL/basic.h>
#include <vector>
#include <sstream>
#include <cstddef>
#include <CGAL/Polynomial/sturm_habicht_sequence.h>
#include <CGAL/Arithmetic_kernel.h>
template<typename ArithmeticKernel> void test_routine() {
// All results are computed with the sres_stha.mpl worksheet
typedef typename ArithmeticKernel::Integer NT;
typedef CGAL::Polynomial<NT> PNT_1;
typedef CGAL::Polynomial<PNT_1> 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<NT> 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<NT>
(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<NT>
(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<NT>
(stha.begin(),stha.end()) == 4);
stha.clear();
typename
CGAL::Polynomial_traits_d<PNT_1>::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<NT>
(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<NT> 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<NT>(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<NT>(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<NT>(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<PNT_1> 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<PNT_1>::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<PNT_1> 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<int>(stha.size());
assert(n == static_cast<int>(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<PNT_1> 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<int>(stha.size());
assert(n == static_cast<int>(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<PNT_1> 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<PNT_1> 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<PNT_1> 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<CGAL::LEDA_arithmetic_kernel>();
#else
std::cout << "LEDA tests skipped" << std::endl;
#endif
#ifdef CGAL_USE_CORE
test_routine<CGAL::CORE_arithmetic_kernel>();
#else
std::cout << "CORE tests skipped" << std::endl;
#endif
return 0;
}

472
Polynomial/test/Polynomial/subresultants.cpp Normal file
View File

@ -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 <mkerber@mpi-inf.mpg.de>
//
// ============================================================================
#include <CGAL/basic.h>
#include <CGAL/Polynomial/subresultants.h>
#include <CGAL/Arithmetic_kernel.h>
template<typename Poly_> Poly_ from_string(const char* s) {
std::stringstream ss(s);
Poly_ f;
ss >> f;
return f;
}
template<typename ArithmeticKernel>
void test_routine() {
typedef ArithmeticKernel Arithmetic_kernel;
typedef typename Arithmetic_kernel::Rational Rational;
typedef typename Arithmetic_kernel::Integer Integer;
typedef CGAL::Polynomial<Integer> Poly_int1;
typedef CGAL::Polynomial<Poly_int1> Poly_int2;
typedef CGAL::Polynomial<Rational> 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<Poly_int1> 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<Poly_int1>::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<Integer> 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<Poly_int1>::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<Poly_int1> 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<Integer> 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<Poly_rat1> 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<Rational> 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<Poly_int2>("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<Poly_int2>("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<Poly_int2> 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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int2>("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<Poly_int1>("P[10(2,-4)(3,-73)(6,97)(7,-62)(10,-56)]");
Poly_int1 g = from_string<Poly_int1>("P[15(5,-75)(6,-17)(8,71)(11,-44)(14,80)(15,-82)]");
std::vector<Poly_int1> sres;
CGAL::CGALi::prs_polynomial_subresultants
(f,g,std::back_inserter(sres));
CGAL_assertion(sres.size()==11);
CGAL_assertion(sres[10]==from_string<Poly_int1>("P[10(2,-39337984)(3,-717918208)(6,953946112)(7,-609738752)(10,-550731776)]"));
CGAL_assertion(sres[9]==from_string<Poly_int1>("P[9(2,-305262755840)(3,-4966105776128)(4,10840151891968)(5,-5962969628672)(6,6702091010048)(7,-22436989575168)(8,19712814514176)(9,-3099911815168)]"));
CGAL_assertion(sres[8]==from_string<Poly_int1>("P[8(2,9680208765108224)(3,156728928112275456)(4,-360058057701515264)(5,152421565222248448)(6,-176106510157656064)(7,746064242175737856)(8,-579306262853885952)]"));
CGAL_assertion(sres[7]==from_string<Poly_int1>("P[7(2,-1486764415342623952896)(3,-26708453688140199103488)(4,8054371431842323176960)(5,3491478546579489027072)(6,38677404054075120663552)(7,-43611415672712033527296)]"));
CGAL_assertion(sres[6]==from_string<Poly_int1>("P[6(2,9979660953345863413381568)(3,193906324691706110661441648)(4,213219021894852075965695648)(5,362883385829743140658019008)(6,-93331232763027671680558576)]"));
CGAL_assertion(sres[5]==from_string<Poly_int1>("P[5(2,-70907443741718133356344565824)(3,-1389117178389007568663227054224)(4,-1747565962047585993090805107776)(5,-2771071359404167222865326731904)]"));
CGAL_assertion(sres[4]==from_string<Poly_int1>("P[4(2,-715876861001937085090535662991104)(3,-13330409138559036348802401660667328)(4,-5258395198219675575175450753910144)]"));
CGAL_assertion(sres[3]==from_string<Poly_int1>("P[3(2,-2842398393513781389004679512270277696)(3,-51817649494648561412362002061296745104)]"));
CGAL_assertion(sres[2]==from_string<Poly_int1>("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<Poly_int1>("P[10(2,-39337984)(3,-717918208)(6,953946112)(7,-609738752)(10,-550731776)]"));
CGAL_assertion(sres[9]==from_string<Poly_int1>("P[9(2,-305262755840)(3,-4966105776128)(4,10840151891968)(5,-5962969628672)(6,6702091010048)(7,-22436989575168)(8,19712814514176)(9,-3099911815168)]"));
CGAL_assertion(sres[8]==from_string<Poly_int1>("P[8(2,9680208765108224)(3,156728928112275456)(4,-360058057701515264)(5,152421565222248448)(6,-176106510157656064)(7,746064242175737856)(8,-579306262853885952)]"));
CGAL_assertion(sres[7]==from_string<Poly_int1>("P[7(2,-1486764415342623952896)(3,-26708453688140199103488)(4,8054371431842323176960)(5,3491478546579489027072)(6,38677404054075120663552)(7,-43611415672712033527296)]"));
CGAL_assertion(sres[6]==from_string<Poly_int1>("P[6(2,9979660953345863413381568)(3,193906324691706110661441648)(4,213219021894852075965695648)(5,362883385829743140658019008)(6,-93331232763027671680558576)]"));
CGAL_assertion(sres[5]==from_string<Poly_int1>("P[5(2,-70907443741718133356344565824)(3,-1389117178389007568663227054224)(4,-1747565962047585993090805107776)(5,-2771071359404167222865326731904)]"));
CGAL_assertion(sres[4]==from_string<Poly_int1>("P[4(2,-715876861001937085090535662991104)(3,-13330409138559036348802401660667328)(4,-5258395198219675575175450753910144)]"));
CGAL_assertion(sres[3]==from_string<Poly_int1>("P[3(2,-2842398393513781389004679512270277696)(3,-51817649494648561412362002061296745104)]"));
CGAL_assertion(sres[2]==from_string<Poly_int1>("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<Poly_int1>("P[3(0,42)(1,-40)(2,-7)(3,-10)]");
Poly_int1 g = from_string<Poly_int1>("P[5(0,74)(1,6)(2,-92)(3,75)(4,23)(5,-50)]");
std::vector<Poly_int1> sres;
CGAL::CGALi::prs_polynomial_subresultants
(f,g,std::back_inserter(sres));
CGAL_assertion(sres.size()==4);
CGAL_assertion(sres[3]==from_string<Poly_int1>("P[3(0,-420)(1,400)(2,70)(3,100)]"));
CGAL_assertion(sres[2]==from_string<Poly_int1>("P[2(0,-1058480)(1,688000)(2,698080)]"));
CGAL_assertion(sres[1]==from_string<Poly_int1>("P[1(0,-22577275200)(1,28253151360)]"));
CGAL_assertion(sres[0]==from_string<Poly_int1>("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<Poly_int1>("P[3(0,-420)(1,400)(2,70)(3,100)]"));
CGAL_assertion(sres[2]==from_string<Poly_int1>("P[2(0,-1058480)(1,688000)(2,698080)]"));
CGAL_assertion(sres[1]==from_string<Poly_int1>("P[1(0,-22577275200)(1,28253151360)]"));
CGAL_assertion(sres[0]==from_string<Poly_int1>("P[0(0,-103066942158720)]"));
}
{ // Test cases with constant polynomials
{
Poly_int1 f(Integer(8));
Poly_int1 g=from_string<Poly_int1>("P[3(0,43)(1,40)(2,-19)(3,2)]");
std::vector<Poly_int1> 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<Integer> 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<Poly_int1>("P[1(0,1)(1,1)]");
std::vector<Poly_int1> 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<Poly_int1> 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<Poly_int2>("P[0(0,P[0(0,122)])]");
Poly_int2 g=from_string<Poly_int2>("P[0(0,P[0(0,0)])]");
std::vector<Poly_int2> 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<Poly_int1>("P[1(0,-5)(1,6)]");
Poly_int1 g = from_string<Poly_int1>("P[1(0,-3)(1,4)]");
std::vector<Poly_int1> 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<Poly_int1>("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<Poly_int1>("P[0(0,2)]"));
}
{ // Test for cofactors, univariate
Poly_int1 f = from_string<Poly_int1>("P[6(0,246)(1,100)(2,197)(3,14)(4,136)(5,191)(6,207)]");
Poly_int1 g = from_string<Poly_int1>("P[5(0,100)(1,394)(2,42)(3,544)(4,955)(5,1242)]");
std::vector<Poly_int1> 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<int>(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_int2> Poly_int3;
Poly_int3 f = from_string<Poly_int3>("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<Poly_int3> 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<int>(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<Poly_int2>("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<Poly_int2>("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<Poly_int2> 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<int>(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<CGAL::LEDA_arithmetic_kernel>();
#else
std::cerr << "LEDA tests skipped!" << std::endl;
#endif
#ifdef CGAL_USE_CORE
test_routine<CGAL::CORE_arithmetic_kernel>();
#else
std::cerr << "CORE tests skipped!" << std::endl;
#endif
return 0;
}