Table of Contents
- Author(s) of this documentation:
- Felix Fontein, Jacques-Olivier Lachaud
- Since
- 0.5
LICENSE: LGPL
Part of Mathematical package.
Module Overview
This module provides classes to represent polynomials, like \( x^3-3x-1 \), and multi-variate polynomials, like \( xy^2-yz \), \( x^4-y^2z^2+z^4 \), etc. The number of variables is fixed at compilation time. Essentially, a n-variate polynomial is defined as a standard polynomial whose coefficients are n-1-variate polynomials.
The following examples are related to this module:
Multivariate polynomial
What is a multivariate polynomial ?
A multivariate monomial is a product of powers of zero, one or several variables, like \( 1, x, x^2, xy, y^3, xy^2z, x^4y^4z^4 \). Sometimes the coefficient in front of the monomial is also included in the definition of the monomial. A multivariate polynomial is a sum of products of a coefficient with a multivariate monomial.
A polynomial is a special case of multivariate polynomial where there is only one variable. The terms variable and indeterminate are synonymous. When you define a multivariate polynomial, you must precise its number of variables. The variables are numbered from 0 to n - 1, if n is the number of variables.
The class MPolynomial
The class MPolynomial is the base class for representing arbitrary multivariate polynomials, i.e. an element of \( K[X_0, ..., X_{-1}] \), where K is some ring. The class MPolynomial<n,K> is parameterized by n, the number of indeterminates, and by K, the type of the coefficients, which should be at least a ring.
If you wish to declare a monovariate polynomial P1(x) over the int ring and a 2-variate polynomial P2(x,y) over the double field, you write the following lines:
In fact, the complete polynomial type is MPolynomial <n, K, A>, where A is an allocator for K, for example std::allocator<K>; this is also the default parameter. Usually this parameter does not need to be changed.
This is implemented recursively: MPolynomial <n, K> is a polynomial in one indeterminate with coefficients in MPolynomial <n-1, K>. The last instance, MPolynomial <0, K>, is essentially equal to K.
Generally, a multivariate polynomial is created by summing multivariate monomials and constants. The following snippet creates the 3-variate polynomial "durchblick", \( Pd(x,y,z) = x^3y+xz^3+y^3z+z^3+5z \).
A constant multivariate polynomial is created simply with:
Creation is thus accomplished using the functions mmonomial<K, A>(e_0, ..., e_{n-1}), where A again is std::allocator<K> by default. This creates a polynomial of type MPolynomial <n, K, A> which consists of exactly one monomial \( X_0^{e_0} * ... * X_{n-1}^{e_{n-1}}\) having coefficient 1.
Otherwise, polynomials are created through operations, derivations, etc.
Interface of the class MPolynomial
Let us describe the interface:
- MPolynomial::isZero: Every polynomial has a predicate isZero(), testing whether the polynomial is the zero polynomial. Note that it is not the same to be a zero polynomial in n variables or to be a zero polynomial in n' variables, n != n'.
- Moreover, it has basic arithmetic, i.e. addition, subtraction, multiplication; this is realized via operator overloading. Assignment and comparison operators are also provided.
- MPolynomial::degree, MPolynomial::leading: Polynomials have a degree function
degree(), returning the degree as a polynomial in the first indeterminate, as well as a methodleading(), returning the leading term (which is of typeconstMPolynomial<n-1, K> &). Note that the degree of the zero polynomial is -1. For instance, \( 1 + xy^2 + x^4yz \) is a degree 4 polynomial in x. Its leading term is \( x^4yz \).
Q(x,y,z)=1+xy^2+x^4yz = (1 + 1 X_1^2 X_0 + 1 X_2 X_1 X_0^4) degree = 4 leading= 1 X_1 X_0
- MPolynomial::operator[]: The i-th coefficient can be obtained by writing f[i] if f is of type MPolynomial <n, K>; then f[i] is of type MPolynomial <n - 1, K>. In case i is larger than the degree of f, a zero polynomial will be returned, and in case f is not constant, the internal space for the coefficients will be enlarged to have space for the i-th coefficient. Note that one should call f.normalize() afterwards if one accessed f[i] for non-const f to ensure that degree and leading coefficients will be computed correctly afterwards.
You may notice that Q[i] are n-1-variate polynomials, i.e., x,y stands for the former y,z variables.
Q[0] = 1
Q[1] = 1 X_0^2
Q[2] = 0
Q[3] = 0
Q[4] = 1 X_1 X_0- MPolynomial::operator(): The polynomial can be evaluated using operator(). If f is of type MPolynomial <n, K, A>, then f(x) with x of type S will return an object of type MPolynomialEvaluator <n, K, A, S>. This object can be casted to MPolynomial <n-1, S> to obtain f with the first indeterminate evaluated as x. Note that S must be a type to which elements of K can be casted. One can also continue evaluating with the MPolynomialEvaluator <n, K, A, S> object, yielding MPolynomialEvaluatorImpl <k, K, ..., A, S> objects, k < n. These object hierarchy ensures that polynomial evaluation is efficient, i.e. after (good enough) optimization of the compiler, is essentially a block of code of (n-k+1) nested for-loops.
double v = Q(1)(2)(3); // evaluation at (x,y,z)=(1,2,3)
- Polynomials can be written to std::ostream's using operator<<. There is also a member function MPolynomial::selfDisplay(std::ostream &) const. Note that the first indeterminate in MPolynomial<n, K> is denoted by X_0, and the last by X_{n-1}.
- Finally, there is a method swap() to swap two polynomial's contents. The polynomials have to be of the same type.This is for instance useful in some algorithms of the STL.
How monomials and coefficients are stored ?
The coefficients are stored in a std::vector<K>; while MPolynomial<1, K> uses std::vector<K>, MPolynomial<n, K> for n > 1 uses std::vector <MPolynomial<n-1, K> *>, i.e. pointers to coefficients are stored. This is implemented using the "intelligent" vector IVector<K> template.
Computing partial derivatives
The library also offers to compute partial derivatives. Given a polynomial f of type MPolynomial<n, K>, one can compute the partial derivative with respect to the j-th variable of f by writing derivative<j>(f). Note that K is required to allow multiplication by int's for this to work.
Q(x,y,z)=1+xy^2+x^4yz = (1 + 1 X_1^2 X_0 + 1 X_2 X_1 X_0^4)
dQ/dx = (1 X_1^2 + 4 X_2 X_1 X_0^3) // y^2+4zyx^3
dQ/dy = (2 X_1 X_0 + 1 X_2 X_0^4) // 2yx+zx^4
dQ/dz = 1 X_1 X_0^4 // yx^4
You may have a look at example polynomial-derivative.cpp to see how to get derivatives of simple polynomials:
Along the same lines, the example polynomial2-derivative.cpp computes first- and second-order partial derivatives of a 2-variate polynomial.
Computing euclidean division and greatest commond divisor (gcd)
Finally, there exist functions euclidDiv() (parameters f, g, q, r), requiring parameters f, g of type const MPolynomial <1, K> & and q, r of type MPolynomial <1, K> &. This function computes q, r of f, g such that degree(r) < degree(g) and f = q*g + r, i.e. it does long division of f by g, storing the quotient in q and the remainder in r. Note that euclidDiv() only works if K is a field, or if everything appearing is divisible by g.leading().
There is another function, gcd() (parameters f, g), accepting f, g of type const MPolynomial <1, K> & and returning a MPolynomial <1, K> object. It computes the monic greatest common divisor of f and g using the euclidean algorithm. In every step, the remainder is made monic. This function is ignoring potential round-off errors, and this function needs that K is (kind of) a field. Use with care.
Efficiency considerations.
Best suited polynomials
The class is tuned so as to be efficient for rather dense polynomials. More precisely, this class is not space nor time efficient for a polynomial like \( 1+x+x^{100000} \). This class is rather efficient for low-degree polynomials, even if there are many monomials. In fact, it factors the computation of the first variables.
For instance, \( x^2y^2+x^2z^2+x^4y^2z+x^4y^2z^3 \) is represented as \( x^2(1(z^2)+y^2(1))+x^4(y^2(z+z^3)) \). As one can see, \( x^2 \) is computed once, which is nice, but \( y^2 \) is computed twice. When you design a polynomial, you should choose the variable order well so as to choose the one that factors at best computations.
Repetitive evaluations in space
Consider the following simple code, which integrates some information in cubic part of the space:
This function takes 206ms for a step 0.01, for 8000000 evaluations.
In fact, you can factor evaluations of the first variables at the beginning of the inner loops.
This function takes 69ms for a step 0.01, for 8000000 evaluations. The C version where the function is explicitely compiled takes 22ms.
Input and output for multivariate polynomials
You may simply output a polynomial an output stream with the usual operator<<.
You may also input a polynomial from an input stream with the usual operator>>. However, it is required to include file "MPolynomialReader.h".
You may read successively several polynomials from the input stream: you just have to separate them with newlines. The input stream is read as long as the input is valid for the multivariate polynomial grammar.
If you wish to handle better error recovery when creating a polynomial, you can use the method MPolynomialReader::read directly. You may look at example polynomial-read.cpp for a more concrete example.
In input, you may write the variables either simply x, y, z, t for the first four variables, or more generically X_0, X_1, ..., X_m, where m is an integer number that is smaller than the number of variables of the polynomial.
For instance, you may test polynomial-derivative.cpp and polynomial2-derivative.cpp as follows:
$ ./examples/math/polynomial-derivative "1+x+x^2-3*x^4" P(X_0) = (1 + 1 X_0 + 1 X_0^2 + -3 X_0^4) P'(X_0) = (1 + 2 X_0 + -12 X_0^3) P''(X_0) = (2 + -36 X_0^2) $ ./examples/math/polynomial-derivative "(2 + -36 X_0^2)" P(X_0) = (2 + -36 X_0^2) P'(X_0) = -72 X_0 P''(X_0) = -72
$ ./examples/math/polynomial2-derivative "(x-3)^2 + (y-2)^2 - 4" P(X_0,X_1) = ((9 + -4 X_1 + 1 X_1^2) + -6 X_0 + 1 X_0^2) dP/dX_0(X_0,X_1) = (-6 + 2 X_0) dP/dX_1(X_0,X_1) = (-4 + 2 X_1) d/dX_1 dP/dX_0(X_0,X_1) = 0 d/dX_0 dP/dX_1(X_0,X_1) = 0 d/dX_0 dP/dX_0(X_0,X_1) = 2 d/dX_1 dP/dX_1(X_0,X_1) = 2
Displaying implicit 3-variate polynomials
If you wish to display surfaces defined as 0-locii of polynomials \( P(x,y,z) \), then you may have a look at example trackImplicitPolynomialSurfaceToOFF.cpp. The following program extract a linear approximation of the given polynomial as a triangulated surface in OFF format.
$ ./examples/topology/trackImplicitPolynomialSurfaceToOFF "(x^2+y^2+2*z^2-1)*(z^2x-0.1)" -2 -2 -2 2 2 2 0.02 # Digital surface has 112826 surfels. # output in marching-cube.off (in 1809ms) # You may display it with your favorite OFF displayer (like geomview, etc). $ ctmviewer marching-cube.off
