DGtal
0.9.4beta

LICENSE: LGPL
Part of Mathematical package.
Module Overview
This module provides classes to represent polynomials, like \( x^33x1 \), and multivariate polynomials, like \( xy^2yz \), \( x^4y^2z^2+z^4 \), etc. The number of variables is fixed at compilation time. Essentially, a nvariate polynomial is defined as a standard polynomial whose coefficients are n1variate polynomials.
The following examples are related to this module:
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 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 2variate 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 <n1, 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 3variate 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_{n1}), 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_{n1}^{e_{n1}}\) having coefficient 1.
Otherwise, polynomials are created through operations, derivations, etc.
Let us describe the interface:
degree()
, returning the degree as a polynomial in the first indeterminate, as well as a method leading()
, returning the leading term (which is of type const
MPolynomial<n1, 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
You may notice that Q[i] are n1variate 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
double v = Q(1)(2)(3); // evaluation at (x,y,z)=(1,2,3)
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<n1, K> *>, i.e. pointers to coefficients are stored. This is implemented using the "intelligent" vector IVector<K> template.
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 jth 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 polynomialderivative.cpp to see how to get derivatives of simple polynomials:
Along the same lines, the example polynomial2derivative.cpp computes first and secondorder partial derivatives of a 2variate polynomial.
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 roundoff errors, and this function needs that K is (kind of) a field. Use with care.
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 lowdegree 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.
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.
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 polynomialread.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 polynomialderivative.cpp and polynomial2derivative.cpp as follows:
$ ./examples/math/polynomialderivative "1+x+x^23*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/polynomialderivative "(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/polynomial2derivative "(x3)^2 + (y2)^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
If you wish to display surfaces defined as 0locii 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^21)*(z^2x0.1)" 2 2 2 2 2 2 0.02 # Digital surface has 112826 surfels. # output in marchingcube.off (in 1809ms) # You may display it with your favorite OFF displayer (like geomview, etc). $ ctmviewer marchingcube.off