User Manual

Author

Michael Hemmer

Fundamentals

Note that this is just a very brief introduction to polynomials. For a quick reference we refer to the Wikipedia or for a more elaborate introduction to any class book on elementary algebra.

A polynomial is either zero, or can be written as the sum of one or more non-zero terms. The number of terms is finite. A term consist of a constant coefficient and a monomial, that is, the product of zero or more variables. Each variable may have an exponent that is a non-negative integer. The exponent on a variable in a term is equal to the degree of that variable in that term. A term with no variables is called a constant term. The degree of a constant term is 0.

For example, \( -7x^3y\) is a term. The coefficient is \( -7\), the monomial is \( x^3y\), comprised of the variables \( x\) and \( y\), the degree of \( x\) is three, and the degree of \( y\) is one. The total degree of the entire term is the sum of the degrees in each variable. In the example above, the degree is \( 3 + 1 = 4\).

A one-variable (univariate) polynomial \( f\) of degree \( n\) has the following form:

\[ f = a_nx^n + a_{n-1}x^{n-1} + ... + a_2x^2 + a_1x + a_0 \]

The coefficient \( a_0\) is called the constant coefficient, \( a_n\) is called the leading coefficient. If \( f\) is not the zero polynomial the leading coefficient is not zero. The polynomial is called monic if \( a_n = 1\). In case the coefficient domain of \( f\) possess a greatest common divisor (gcd) the content of \( f\) is the gcd of all coefficients of \( f\). For instance, the content of \( 12 x^3 + 6\) is \( 6\).

A multivariate polynomial is a polynomial in more than one variable. According to the number of variables it is possible to further classify multivariate polynomials as bivariate, trivariate etc. In contrast to univariate polynomials the terms of a multivariate polynomial are not completely ordered by their total degree. However, given a certain order on the variables it is possible to define a lexicographic order on the terms. Given this order the leading coefficient of a multivariate polynomial is defined as the coefficient of the highest term. For instance the leading coefficient of the multivariate polynomial \( p = 5 x^3y + 7xy^2\) is \( 7\), given that \( y\) has an higher order than \( x\).

However, it is also possible to interpret a multivariate polynomial as a univariate polynomial in that variable. For instance the trivariate polynomial

\[ q = x^5 + 7x^2y^1z^2 + 13x^1y^2z^2 \in \Z[x,y,z] \]

may be interpreted as a univariate polynomial in \( z\), that is, \( q\) is interpreted as an element of \( R[z]\), with \( R=\Z[x,y]\).

\[ q = (13x^1y^2 + 7x^2y^1)z^2 + x^5z^0 \in R[z] \]

In this case the leading coefficient of \( q\) with respect to \( z\) is \( 13x^1y^2 + 7x^2y^1\) and \( x^5\) becomes the 'constant' term.

A homogeneous polynomial is a polynomial whose terms do all have the same total degree. For example, \( h = x^5 + 7x^2y^1z^2 + 13x^1y^2z^2\) is a homogeneous polynomial of degree \( 5\), in three variables.

General Design

The package introduces a concept Polynomial_d, a concept for multivariate polynomials in \( d\) variables. Though the concept is written for an arbitrary number of variables, the number of variables is considered as fixed for a particular model of Polynomial_d. The concept also allows univariate polynomials.

First of all a model of Polynomial_d is considered as an algebraic structure, that is, the ring operations \( \{+, -, \cdot\}\) are provided due to the fact that Polynomial_d refines at least the concept IntegralDomainWithoutDivision. However, a model of Polynomial_d has to be accompanied by a traits class Polynomial_traits_d<Polynomial_d> being a model of PolynomialTraits_d. This traits class provides all further functionalities on polynomials.

Given a \( d\)-variate polynomial over some base ring \( R\) there are at least two different possible views on such a polynomial.

• The recursive or univariate view: In this view, a polynomial is considered as an element of \( R[x_0,\dots,x_{d-2}][x_{d-1}]\). That is, the polynomial is treated as a univariate polynomial over the ring \( R[x_0,\dots,x_{d-2}]\).
• The symmetric or multivariate view: This view is almost symmetric with respect to all variables. It considers the polynomial as an element of \( R [x_0,\dots,x_{d-1}]\).

According to these two different views the traits class is required to provide two different coefficient types:

• Polynomial_traits_d::Coefficient_type representing \( R[x_0,\dots,x_{d-2}]\).
• Polynomial_traits_d::Innermost_coefficient_type representing the base ring \( R\).

Another important type which is introduced by this package is Exponent_vector. It is derived from std::vector<int> and used to identify multivariate monomials. For instance the exponent vector containing the sequence \( [3,2,4]\) corresponds to the trivariate monomial \( x_0^3x_1^2x_2^4\). Note that a vector with negative exponents is considered as invalid. However, we allow negative exponents as they may appear as intermediate results, in particular we did not derive from std::vector<unsigned int>.

Constructing a Multivariate Polynomial

First of all the concept Polynomial_d requires that the model is constructible from int. This is due to the fact that Polynomial_d refines IntegralDomainWithoutDivision which in turn refines FromIntConstructible. Of course this allows only the construction of constant polynomials.

In general a polynomial is constructed using the functor Polynomial_traits_d::Construct_polynomial a model of PolynomialTraits_d::ConstructPolynomial. Basically there are two options:

• The polynomial is constructed from an iterator range with value type Polynomial_traits_d::Coefficient_type, where the begin iterator refers to the constant term (constant with respect to the outermost variable).
• The polynomial is constructed from an iterator range with value type std::pair<Exponent_vector, Polynomial_traits_d::Innermost_coefficient_type>, where each pair defines the coefficient for the monomial defined by the exponent vector.

However, in some cases it might be more convenient to just construct the polynomials representing the different variables and to obtain the final polynomial using algebraic expressions. The most elegant way to construct a certain variable is Polynomial_traits_d::Shift being a model of PolynomialTraits_d::Shift.

Example

The following example illustrates different ways to construct a bivariate polynomial:
File Polynomial/construction.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
 
  typedef CGAL::Polynomial_type_generator<int,2>::Type Poly_2;
  typedef CGAL::Polynomial_traits_d<Poly_2>            PT_2;
  typedef PT_2::Coefficient_type                       Poly_1;
  typedef PT_2::Innermost_coefficient_type             Integer;
 
  PT_2::Construct_polynomial construct_polynomial;
  Poly_2 dc;
 
 
  // constructing a constant polynomial from int
  Poly_2 two(2); // = 2
  std::cout << "A constant polynomial: " << two << std::endl;
 
 
  // construction from an iterator range of univariate polynomials
 
  std::list<Poly_1> univariate_coeffs;
  univariate_coeffs.push_back(Poly_1(3));
  univariate_coeffs.push_back(Poly_1(0));
  univariate_coeffs.push_back(Poly_1(5));
  Poly_2 F = // 5*y^2 + 3
    construct_polynomial(univariate_coeffs.begin(),univariate_coeffs.end());
  std::cout << "The bivariate polynomial F: " << F << std::endl;
 
 
  // construction from an iterator range over monomials
 
  std::list<std::pair<CGAL::Exponent_vector, Integer> > innermost_coeffs;
  innermost_coeffs.push_back(std::make_pair(CGAL::Exponent_vector(0,0),-2));
  innermost_coeffs.push_back(std::make_pair(CGAL::Exponent_vector(3,5),2));
  Poly_2 G = // (2*x^3)*y^5 + (-2)
    construct_polynomial(innermost_coeffs.begin(),innermost_coeffs.end());
  std::cout << "The bivariate polynomial G: " << G << std::endl;
 
  //construction using shift
  PT_2::Shift shift;
  Poly_2 x = shift(Poly_2(1),1,0); // 'multiply' 1 by x_0^1
  Poly_2 y = shift(Poly_2(1),1,1); // 'multiply' 1 by x_1^1
 
  Poly_2 H = 5 * x * y + 3 * y * y; // = 3*y^2 + (5*x)*y
  std::cout << "The bivariate polynomial H: " << H << std::endl;
}

Coefficient Access

In order to obtain a certain coefficient the traits class provides several functors. Note that the functors do not allow a write access to the coefficients.

• PolynomialTraits_d::GetCoefficient: a model of this concept provides access to a coefficient in the univariate view, that is, it returns elements of \( R[x_0,\dots,x_{d-2}]\).
• PolynomialTraits_d::GetInnermostCoefficient: a model of this concept provides access to a coefficient in the multivariate view, that is, it returns elements of \( R\).
• PolynomialTraits_d::LeadingCoefficient: a model of this concept provides access to the leading coefficient in the univariate view.
• PolynomialTraits_d::InnermostLeadingCoefficient: a model of this concept provides access to the leading coefficient in the multivariate view, that is, it returns the (innermost) coefficient of the leading multivariate monomial. See also PolynomialTraits_d::DegreeVector.

Example

The following example illustrates the application of the functors discussed above:
File Polynomial/coefficient_access.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
  typedef CGAL::Polynomial_type_generator<int,2>::Type Poly_2;
  typedef CGAL::Polynomial_traits_d<Poly_2>            PT_2;
 
  //construction using shift
  Poly_2 x = PT_2::Shift()(Poly_2(1),1,0); // = x^1
  Poly_2 y = PT_2::Shift()(Poly_2(1),1,1); // = y^1
 
  Poly_2 F // = (11*x^2 + 5*x)*y^4 + (7*x^2)*y^3
    = 11 * CGAL::ipower(y,4) * CGAL::ipower(x,2)
    + 5 * CGAL::ipower(y,4)  * CGAL::ipower(x,1)
    + 7 * CGAL::ipower(y,3)  * CGAL::ipower(x,2);
  std::cout << "The bivariate polynomial F: " << F <<"\n"<< std::endl;
 
  PT_2::Get_coefficient get_coefficient;
  std::cout << "Coefficient of y^0: "<< get_coefficient(F,0) << std::endl;
  std::cout << "Coefficient of y^1: "<< get_coefficient(F,1) << std::endl;
  std::cout << "Coefficient of y^2: "<< get_coefficient(F,2) << std::endl;
  std::cout << "Coefficient of y^3: "<< get_coefficient(F,3) << std::endl;
  std::cout << "Coefficient of y^4: "<< get_coefficient(F,4) << std::endl;
  std::cout << "Coefficient of y^5: "<< get_coefficient(F,5) << std::endl;
  std::cout << std::endl;
 
  PT_2::Leading_coefficient lcoeff;
  std::cout << "Leading coefficient with respect to y:           "
            << lcoeff(F)   // = 11*x^2 + 5*x
            << std::endl;
 
  PT_2::Get_innermost_coefficient get_icoeff;
  std::cout << "Innermost coefficient of monomial x^1y^4:        "
            << get_icoeff(F,CGAL::Exponent_vector(1,4)) // = 5
            << std::endl;
 
  PT_2::Innermost_leading_coefficient ilcoeff;
  std::cout << "Innermost leading coefficient with respect to y: "
            << ilcoeff(F) // = 11
            << std::endl;
}

Degree, Total Degree and Degree Vector

There are three functors in PolynomialTraits_d related to the degree of a polynomial.

• PolynomialTraits_d::Degree: a model of this concept returns the degree of the polynomial in the univariate view. By default this is the degree with respect to the outermost variable, but it is also possible to select another variable.
• PolynomialTraits_d::TotalDegree: a model of this concept returns the total degree of a polynomial. The polynomial is considered as a multivariate polynomial. The total degree is the maximum over the sums of the exponents of each multivariate monomial.
• PolynomialTraits_d::DegreeVector: a model of this concept returns the exponent vector of the leading monomial, where the monomial order is lexicographic and starts with the outermost variable. See also PolynomialTraits_d::InnermostLeadingCoefficient.

Example

The following example illustrates the application of the functors discussed above:
File Polynomial/degree.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
  typedef CGAL::Polynomial_type_generator<int,2>::Type Poly_2;
  typedef CGAL::Polynomial_traits_d<Poly_2>            PT_2;
 
  //construction using shift
  Poly_2 x = PT_2::Shift()(Poly_2(1),1,0); // x_0^1
  Poly_2 y = PT_2::Shift()(Poly_2(1),1,1); // x_1^1
 
  Poly_2 F // = (11*x^2 + 5*x)*y^4 + (7*x^2)*y^3
    = 11 * CGAL::ipower(y,4) * CGAL::ipower(x,2)
    + 5 * CGAL::ipower(y,4)  * CGAL::ipower(x,1)
    + 7 * CGAL::ipower(y,3)  * CGAL::ipower(x,2);
  std::cout << "The bivariate polynomial F: " << F <<"\n"<< std::endl;
 
  PT_2::Degree degree;
  PT_2::Total_degree total_degree;
  PT_2::Degree_vector degree_vector;
 
  std::cout << "The degree of F with respect to y: "<< degree(F)       // = 4
            << std::endl;
  std::cout << "The degree of F with respect to x: "<< degree(F,0)     // = 2
            << std::endl;
  std::cout << "The total degree of F            : "<< total_degree(F) // = 6
            << std::endl;
  std::cout << "The degree vector of F           : "<< degree_vector(F)// = (2,4)
            << std::endl;
 
 
 
}

Changing the Order of Variables

Given for instance a bivariate polynomial it is conceivable that one wants to interchange the role of \( x\) and \( y\). That is one wants to interpret the \( x\) as \( y\) and vice versa. For such a case the polynomial traits provides PolynomialTraits_d::Swap:

Given a polynomial \( p\) and to two indices \( i\) and \( j\), the functor returns the polynomial in which \( x_i\) is substituted by \( x_j\) and vice versa, that is, the variables swap their positions. The order of the other variables remains untouched.

Another scenario is, that a particular variable should be moved to another position, for instance, it should become the outermost variable while the relative order of the other variables remains unchanged. For such a case the polynomial traits provides PolynomialTraits_d::Move.

Of course there is also a general method to interchange the order of variables, namely PolynomialTraits_d::Permute.

Example

The following example illustrates the application of the functors discussed above:
File Polynomial/swap_move.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
  typedef CGAL::Polynomial_type_generator<int,3>::Type Poly_3;
  typedef CGAL::Polynomial_traits_d<Poly_3>            PT_3;
 
  //construction using shift
  Poly_3 x = PT_3::Shift()(Poly_3(1),1,0); // x_0^1
  Poly_3 y = PT_3::Shift()(Poly_3(1),1,1); // x_1^1
  Poly_3 z = PT_3::Shift()(Poly_3(1),1,2); // x_2^1
 
 
  Poly_3 F = x*y*y*z*z*z;
  std::cout << "The trivariate polynomial F: " << F << std::endl;
  std::cout << std::endl;
 
  PT_3::Swap swap;
  PT_3::Move move;
  PT_3::Permute permute;
 
  std::cout << "x and z swapped: "<< swap(F,0,2) // = x^3*y^2*z
            << std::endl;
  std::cout << "x and y swapped: "<< swap(F,0,1) // = x^2*y*z^3
            << std::endl << std::endl;
 
  std::cout << "x moved to outermost position           : "
            << move(F,0,2)                       // = x^2*y^3*z
            << std::endl;
  std::cout << "Same as swap(swap(F,0,1),1,2)           : "
            << swap(swap(F,0,1),1,2)             // = x^2*y^3*z
            << std::endl;
 
  std::cout << "Same as the permutation (0,1,2)->(2,0,1): ";
  std::vector<int> perm;
  perm.push_back(2);perm.push_back(0);perm.push_back(1);
  std::cout << permute(F,perm.begin(),perm.end())// = x^2*y^3*z
            << std::endl;
 
}

GCD and More

Since the concept PolynomialTraits_d refines the concept AlgebraicStructureTraits the polynomial traits provides functors for integral division, division with remainder, greatest common divisor, etc. But note that the algebraic structure of a polynomial depends on the algebraic structure of the innermost coefficient, for instance, a gcd is available if and only if the innermost coefficient is a Field or a UniqueFactorizationDomain. Hence, we can not provide a \( gcd\) if the innermost coefficient is just an IntegralDomain since it is simply not well definedAn example for such a number type is the template Sqrt_extension<NT,ROOT> representing an algebraic extension of degree two. This is just an IntegralDomain if NT is not a Field. . However, if we would consider the polynomial over the quotient field of the integral domain the \( gcd\) would be well defined. The only problem is that the result can not be represented over the ring since it contains denominators. Therefore, the PolynomialTraits_d requires functors such as PolynomialTraits_d::GcdUpToConstantFactor. This functor computes the gcd of two polynomials up to a constant factor (utcf). That is, it returns the correct gcd for polynomials over the quotient field, but multiplied by some constant such that the result is representable with coefficients in the ring.

However, note that these 'utcf' functions are usually a bit faster than their strict counterparts. This is due to the fact that the 'utcf' functions are allowed to skip the computation of the correct constant factor. Note that in many cases the constant factor is in fact not needed. In particular if the polynomials are supposed to represent some zero set, that is, an algebraic curve or surface.

The concepts for the related functors are:

• AlgebraicStructureTraits::Gcd

  PolynomialTraits_d::GcdUpToConstantFactor

• AlgebraicStructureTraits::IntegralDivision

  PolynomialTraits_d::IntegralDivisionUpToConstantFactor

• PolynomialTraits_d::UnivariateContent

  PolynomialTraits_d::UnivariateContentUpToConstantFactor

• PolynomialTraits_d::SquareFreeFactorize

  PolynomialTraits_d::SquareFreeFactorizeUpToConstantFactor

Another analog functionality is the pseudo division. The related functors replace the usual division with remainder in case the Polynomial is not a EuclideanRing.

The concepts for the related functors are:

• AlgebraicStructureTraits::Div_mod

  PolynomialTraits_d::PseudoDivision

• AlgebraicStructureTraits::Div

  PolynomialTraits_d::PseudoDivisionQuotient

• AlgebraicStructureTraits::Mod

  PolynomialTraits_d::PseudoDivisionRemainder

Example

The following example illustrates the application of some functors discussed above:
File Polynomial/gcd_up_to_constant_factor.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
  typedef CGAL::Polynomial_type_generator<int,1>::Type Poly_1;
  typedef CGAL::Polynomial_traits_d<Poly_1>            PT_1;
 
  PT_1::Shift                     shift;
  PT_1::Gcd                       gcd;
  PT_1::Gcd_up_to_constant_factor gcd_utcf;
  PT_1::Multivariate_content      mcontent;
  PT_1::Canonicalize              canonicalize;
 
  //construction using shift
  Poly_1 x = shift(Poly_1(1),1,0); // x^1
 
  // common factor 7 * (x^2-2)
  Poly_1 F = 21*(x-5)*(x*x-2); // = 21*x^3 + (-105)*x^2 + (-42)*x + 210
  Poly_1 G = 14*(x-3)*(x*x-2); // = 14*x^3 + (-42)*x^2 + (-28)*x + 84
 
  std::cout << "The univariate polynomial F: " << F << std::endl;
  std::cout << "The univariate polynomial G: " << G << std::endl;
  std::cout << "Common multivariate content:              "
            << CGAL::gcd(mcontent(F),mcontent(G)) // = 7
            << std::endl;
  std::cout << "The gcd of F and G:                       "
            << gcd(F,G)                           // = 7*x^2 + (-14)
            << std::endl;
  std::cout << "The gcd up to constant factor of F and G: "
            << gcd_utcf(F,G)                      // = x^2 + (-2)
            << std::endl;
  std::cout << "Same as canonicalized gcd of F and G:     "
            << canonicalize(gcd_utcf(F,G))        // = x^2 + (-2)
            << std::endl;
 
}

Evaluation and Substitution

Of course, it should also be possible to evaluate a polynomial or substitute its variables. We also require a special functor to determine whether a polynomial is zero at a given point. In case the inner most coefficient is RealEmbeddable the traits also must provide a function to compute the sign at a given point.

The concepts for the related functors are:

• PolynomialTraits_d::Substitute
• PolynomialTraits_d::Evaluate
• PolynomialTraits_d::IsZeroAt
• PolynomialTraits_d::SignAt

The traits is also required to provide variants of these functors that interpret the polynomial as a homogeneous polynomial by adding a virtual homogeneous variable such that each term has the same degree, namely the degree of the polynomial. Of course there is a difference between the univariate and multivariate view. For instance the polynomial

\[ 5x^3 + 7x - 3 \]

has degree 3, hence it is interpreted as the homogeneous polynomial

\[ 5x^3 + 7xw^2 -3w^3 \]

by adding the homogeneous variable \( w\). In case of the multivariate view each term is filled up by the homogeneous variable such that the degree of each term is equal to the total degree of the polynomial. Note that these functors may significantly improve efficiency. For instance, it is possible to determine the sign of a polynomial over integer coefficients at a rational point without changing the coefficient domain of the polynomial. For more details have a look at the following concepts:

• PolynomialTraits_d::SubstituteHomogeneous
• PolynomialTraits_d::EvaluateHomogeneous
• PolynomialTraits_d::IsZeroAtHomogeneous
• PolynomialTraits_d::SignAtHomogeneous

Note that substitute allows the substitution of the variables by any type that is ExplicitInteroperable with the innermost coefficient type. This is a very powerful tool since it allows the substitution of the variables by polynomials. However, for some standard manipulations such as translation or scaling we require special functors since they are expected to be faster than their equivalent implementation using substitution:

• PolynomialTraits_d::Shift
• PolynomialTraits_d::Negate
• PolynomialTraits_d::Invert
• PolynomialTraits_d::Translate
• PolynomialTraits_d::TranslateHomogeneous
• PolynomialTraits_d::Scale
• PolynomialTraits_d::ScaleHomogeneous

Example

The following example illustrates the application of some functors discussed above:
File Polynomial/substitute.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
  typedef CGAL::Polynomial_type_generator<int,2>::Type Poly_2;
  typedef CGAL::Polynomial_traits_d<Poly_2>            PT_2;
 
  //construction using shift
  Poly_2 x = PT_2::Shift()(Poly_2(1),1,0); // x^1
  Poly_2 y = PT_2::Shift()(Poly_2(1),1,1); // y^1
 
 
  Poly_2 F = 2*x*y + 3*CGAL::ipower(y,3);
  std::cout << "The bivariate polynomial F: " << F // = 3*y^3 + (2*x)*y
            << std::endl << std::endl;
 
  PT_2::Evaluate evaluate;
  PT_2::Evaluate_homogeneous hevaluate;
 
  // Evaluation considers a polynomials as univariate:
  std::cout << "F(5): " << evaluate(F,5)      // = 10*x + 375
            << std::endl;
  // Evaluate_homogeneous considers F as a homogeneous polynomial in
  // the outermost variable only, that is, F is interpreted as
  // F(u,v) = 2*x*u*v^2 + 3 * u^3
  std::cout << "F(5,7): " << hevaluate(F,5,7) // = 490*x + 375
            << std::endl << std::endl;
 
 
  PT_2::Substitute substitute;
  PT_2::Substitute_homogeneous hsubstitute;
 
  // Substitute considers a polynomials as multivariate, that is, the
  // new values for the variables are given by an iterator range
  // Note that the value type only has to be interoperable with the innermost
  // coefficient
  std::list<Poly_2> replacements;
  replacements.push_back(x-1); // replace x by x-1
  replacements.push_back(y);   // replace y by y, i.e., do nothing
 
  std::cout << "The bivariate polynomial F: " << F // = 3*y^3 + (2*x)*y
            << std::endl;
  std::cout << "F(x-1,y):   " // = 3*y^3 + (2*x + (-2))*y
            << substitute(F,replacements.begin(),replacements.end())
            << std::endl;
  // Substitute_homogeneous considers F as a homogeneous polynomial in
  // all variable, that is, F is interpreted as
  // F(x,y,w) = 2*x*y*w + 3 * y^3
  replacements.push_back(y);  // replace z by y
 
  std::cout << "F(x-1,y,y): " // = 3*y^3 + (2*x + (-2))*y^2
            << hsubstitute(F,replacements.begin(),replacements.end())
            << std::endl;
}

Resultants, Subresultants and Sturm-Habicht Sequences

The PolynomialTraits_d concept also provides more sophisticated functors for computations with polynomials - computing the resultant of two polynomials, their polynomial subresultant sequence, with or without cofactors, and their principal subresultant coefficients.

• PolynomialTraits_d::Resultant
• PolynomialTraits_d::PolynomialSubresultants
• PolynomialTraits_d::PolynomialSubresultantsWithCofactors
• PolynomialTraits_d::PrincipalSubresultants

Moreover, functors to compute the Sturm-Habicht sequence, with or without cofactors, and for the principal Sturm-Habicht coefficients exist.

• PolynomialTraits_d::SturmHabichtSequence
• PolynomialTraits_d::SturmHabichtSequenceWithCofactors
• PolynomialTraits_d::PrincipalSturmHabichtSequence

For a formal definition of all used terms, we refer to the corresponding reference pages.

The principal Sturm-Habicht sequence allows to count the number of real roots of a polynomial using the function

• number_of_real_roots().

As input, this function requires an iterator range that represents the principal Sturm-Habicht coefficients. This might look complicated at a first sight, as one has to store the principal Sturm-Habicht sequence temporarily. However, we remark an important property of the (principal) Sturm-Habicht sequence. Having a polynomial \( f_t(x)\) that depends on a parameter \( t\), and its (principal) Sturm-Habicht coefficients \( \mathrm{stha}_0(f_t),\ldots,\mathrm{stha}_n(f_t)\), evaluating \( \mathrm{stha}_0(f_t)\) for \( t=t_0\) yields a valid (principal) Sturm-Habicht sequence for \( f_{t_0}\). The same holds for (principal) subresultants. Thus, it is enough in such situations to compute the sequence once for the parameter \( t\), and call number_of_real_roots() for each specialized parameter value.

We finally remark that computing subresultants and Sturm-Habicht sequences introduces an enormous coefficient blow-up. An application of the functors therefore does not make sense for built-in integers except for toy examples. To avoid overflows, one should use arbitrary size integer types in real applications.

Example

The following example illustrates how two compute resultants of two polynomials, and how to count the number of distinct real roots of a polynomial using its principal Sturm-Habicht coefficients.

File Polynomial/subresultants.cpp

#include <CGAL/Polynomial.h>
#include <CGAL/Polynomial_traits_d.h>
#include <CGAL/Polynomial_type_generator.h>
 
#include <CGAL/Exact_integer.h>
 
int main(){
  CGAL::IO::set_pretty_mode(std::cout);
 
  typedef CGAL::Exact_integer Int;
 
  typedef CGAL::Polynomial_type_generator<Int,1>::Type Poly_1;
  typedef CGAL::Polynomial_traits_d<Poly_1>            PT_1;
 
  //construction using shift
  Poly_1 x = PT_1::Shift()(Poly_1(1),1); // x^1
 
  Poly_1 F // = (x+1)^2*(x-1)*(2x-1)=2x^4+x^3-3x^2-x+1
    =   2 * CGAL::ipower(x,4) + 1 * CGAL::ipower(x,3)
      - 3 * CGAL::ipower(x,2) - 1 * CGAL::ipower(x,1)
      + 1 * CGAL::ipower(x,0);
  std::cout << "F=" << F << std::endl;
 
  Poly_1 G // = (x+1)*(x+3)=x^2+4*x+3
    =   1 * CGAL::ipower(x,2) + 4 * CGAL::ipower(x,1) + 3 * CGAL::ipower(x,0);
  std::cout << "G=" << G << std::endl;
 
  // Resultant computation:
  PT_1::Resultant resultant;
 
  std::cout << "The resultant of F and G is: " << resultant(F,G) << std::endl;
  // It is zero, because F and G have a common factor
 
  // Real root counting:
  PT_1::Principal_sturm_habicht_sequence stha;
  std::vector<Int> psc;
 
  stha(F,std::back_inserter(psc));
 
  int roots = CGAL::number_of_real_roots(psc.begin(),psc.end());
 
  std::cout << "The number of real roots of F is: " << roots << std::endl; // 3
 
  roots =  CGAL::number_of_real_roots(G);
 
  std::cout << "The number of real roots of G is: " << roots << std::endl; // 2
 
  return 0;
 
}

Design and Implementation History

This package is the result of the integration process of the NumeriX library of Exacus [1] into CGAL.

The class Polynomial<Coeff> had been started by Michael Seel within CGAL as part of the Nef_2 package. As part of the Exacus project it got significantly improved by Arno Eigenwillig and Michael Hemmer.

However, due to the recursive definition the class was rather restricted to the univariate view. Moreover, it is clear that depending on the context other classes that are symmetric in all variables or dedicated for sparse polynomials may be more efficient. As a consequence this package introduced the Polynomial_traits_d<Polynomial_d> giving also the symmetric view on polynomials and the opportunity to introduce and use other classes representing polynomials within CGAL.