Univariate and multivariate polynomials

A basic implemententation is proposed for handling univariate and multivariate polynomials (represented by their coefficients). You should use those classes carefully for two reasons :

The implementation is not optimal for multivariate polynomials (loops are not unrolled as for tiny vectors). Maybe this drawback will be removed in future versions of Montjoie
The use of Horner's algorithm may lead to instability when the order of polynomial becomes too large. For example, you should not use those classes for Legendre polynoms, but rather the functions GetJacobiPolynomial and EvaluateJacobiPolynomial, or the class LegendrePolynomial detailed below.

So these classes are useful for developing, testing new methods. Here we briefly describe how to use them.

Basic Use

First, we show examples for 1-D polynomials (univariate) :

// first we create two useful polynoms : 1 and x
UnivariatePolynomial<double> x(1), one(0);
x(0) = 0; x(1) = 1.0;
one(0) = 1.0;

// then you can construct other polynomials with those basic polynoms
// for example chebyshev polynomials are constructed  :
UnivariatePolynomial<double> Pn_m1, Pn, Pn_p1;
Pn_m1 = one;
Pn = x;
int n = 10;
for (int k = 1; k < n; k++)
  {
    Pn_p1 = 2.0*x*Pn - Pn_m1;
    Pn_m1 = Pn;
    Pn = Pn_p1;
  }

// of course you can enter directly coefficients
// P = a0 + a1 X + a2 X^2 + .. + an x^n
UnivariatePolynomial<double> P;
P.SetOrder(n);
P(0) = a0;
P(1) = a1;
// ...
P(n) = an;

// and evaluate y = P(x) for a given x with Horner algorithm :
double y = P.Evaluate(2.53);

// and operators -, +, * can be used 
UnivariatePolynomial<double> Q, R, S;
S = 2.5*P*Q + R;

// polynoms can be derived or integrated (with 0-constant)
UnivariatePolynomial<double> dP_dx, Primitive_P;
DerivatePolynomial(P, dP_dx);
IntegratePolynomial(P, Primitive_P);

For multivariate polynomials, the degree of the polynom is the sum of the degree for each variable.

// for 3-D polynoms, we create 1, x, y and z
// the first argument of the constructor is the number of variables
// by default, we set all the coefficients to 0
MultivariatePolynomial<double> x(3, 1), one(3, 0), y(3, 1), z(3,
1);
one(0, 0, 0) = 1.0
x(1, 0, 0) = 1.0;
y(0, 1, 0) = 1.0;
z(0, 0, 1) = 1.0;

// then you can construct other polynomials with those basic polynoms
MultivariatePolynomial<double> P;
P = x*(x + 2.5*y*z); 

// of course you can enter directly coefficients
P.SetOrder(3, n);
P(0, 0, 0) = a000;
P(1, 0, 0) = a100;
// ...
P(n, 0, 0) = an00;
P(0, n, 0) = a0n0;
P(0, 0, n) = a00n;

// and evaluate y = P(x) for a given x with Horner algorithm :
TinyVector<double, 3> pt(-3.4, 2.31, 0.56);
double value = P.Evaluate(pt);

// and operators -, +, * can be used 
MultivariatePolynomial<double> Q, R, S;
S = 2.5*P*Q + R;

// polynoms can be derived or integrated (with 0-constant)
MultivariatePolynomial<double> dP_dx, dP_dy, Primitive_Px;
DerivatePolynomial(P, dP_dx, 0);
DerivatePolynomial(P, dP_dy, 1);
IntegratePolynomial(P, Primitive_Px, 0);

// operator () has been overloaded for 2-D and 3-D polynoms
// for higher dimensions, you have to use IndexToPowers, PowersToIndex
// for example in 4-D :
MultivariatePolynomial<double> t;
one.SetOrder(4, 1);
x.SetOrder(4, 1);
y.SetOrder(4, 1);
z.SetOrder(4, 1);
t.SetOrder(4, 1);

IVect powers(4);
powers.Fill(0);
int pos = one.PowersToIndex(powers);
one(pos) = 1.0;

powers(0) = 1;
pos = x.PowersToIndex(powers);
x(pos) = 1.0;

powers(0) = 0; powers(1) = 1;
pos = y.PowersToIndex(powers);
y(pos) = 1.0;

powers(1) = 0; powers(2) = 1;
pos = z.PowersToIndex(powers);
z(pos) = 1.0;

powers(2) = 0; powers(3) = 1;
pos = t.PowersToIndex(powers);
t(pos) = 1.0;

Methods for univariate polynomials

constructors
operators
GetOrder	returns the degree of polynomial
SetOrder	modifies the degree of polynomial and removes previous coefficients
ResizeOrder	modifies the degree of polynomial and keeps previous coefficients
Zero	sets all coefficients to 0
Evaluate	evaluates the polynomial with Horner algorithm

Methods for the class LegendrePolynomial

constructor
EvaluateDerivative	evaluates P_n(cos θ) and its derivatives
EvaluatePn	evaluates Legendre polynomials at a given point
EvaluateDerive	evaluates Legendre polynomials and derivatives

Methods for the class AssociatedLegendrePolynomial

Init	inits the computation of associated Legendre polynomials
EvaluatePnm	evaluates associated Legendre polynomials

Functions for univariate polynomials

Pow	Computes a power of the polynomial
DerivatePolynomial	Computes the derivative of a polynomial
IntegratePolynomial	Computes the primitive of a polynomial
ComputeLegendrePolynome	Computes Legendre polynomials
Add	adds two polynomials
Mlt	computes the product of two polynomials
MltAdd	adds the product of two polynomals
SolvePolynomialEquation	finds the roots of a polynomial

Methods for multivariate polynomials

constructors
operators
GetM	returns the number of coefficients of the polynomial
GetDimension	returns the number of variables (called dimension)
GetOrder	returns the degree of polynomial
SetOrder	modifies the degree and dimension of polynomial and removes previous coefficients
ResizeOrder	modifies the degree and dimension of polynomial and keeps previous coefficients
Zero	sets all coefficients to 0
Evaluate	evaluates the polynomial with Horner algorithm
PowersToIndex	computes scalar index from the tuple (i, j, ...)
IndexToPowers	computes the tuple (i, j, ...) from the scalar index
PrintRational	displays coefficients as fractions
PrintMonomial	displays a single monomial of the polynomial

Functions for multivariate polynomials

Pow	Computes a power of the polynomial
DerivatePolynomial	Computes the derivative of a polynomial
Derivate	Computes the derivative of a polynomial
GetCurlPolynomial	Computes the curl of a 3-D vector of polynomials
GetDivPolynomial	Computes the divergence of a 3-D vector of polynomials
IntegratePolynomial	Computes the primitive of a polynomial
Mlt	computes the product of a polynomial with a scalar
MltAdd	adds the product of two polynomals
IntegrateOnSimplex	computes the integral of the polynomial over the simplex (triangle in 2-D, tetrahedron in 3-D)
Add	adds two polynomials
ScalePolynomVariable	all the variables of the polynomial are multiplied by the same coefficient
CompressPolynom	sets to 0 small coefficients of the polynomial
OppositeX	x-variable is replaced by its opposite
OppositeY	y-variable is replaced by its opposite
PermuteXY	x-variable and y-variable are swapped
GenerateSymPol	generates polynomials obtained by symmetry with respect to diagonals of the square [-1, 1]²
GenerateSym4_Edges	generates polynomials obtained by symmetry with respect to axis x=0 and y=0 of the square [-1, 1]²
GenerateSym8	generates polynomials obtained by symmetry with respect to diagonals and axis x=0 and y=0 of the square [-1, 1]²

Constructors for UnivariatePolynomial

Syntax

  UnivariatePolynomial();
  UnivariatePolynomial(int order);

The default constructor creates a polynomial equal to 0 (of degree 0). The second constructor takes the degree of the polynomial as argument and initializes all coefficients to 0.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

Location :

Share/UnivariatePolynomialInline.cxx

Operators for UnivariatePolynomial

The operators +, -, * are implemented for univariate polynomials and as well as operators +=, -= and *=. The operator * can mix scalars with polynomials.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2), Q(3), R;
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// Q = x^3 - 7 x
Q(3) = 1.0; Q(1) = -7.0;

// R = 2.0 P Q - 3.0 P
R = 2.0*P*Q - 3.0*P;

// Q = Q + R
Q += R;

Location :

Share/UnivariatePolynomialInline.cxx

GetOrder

Syntax

  int GetOrder() const;

This method returns the degree of the polynomial.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// P.GetOrder() should return 2
int deg = P.GetOrder();

Location :

Share/UnivariatePolynomialInline.cxx

SetOrder

Syntax

  void SetOrder(int k);

This method modifies the degree of the polynomial and sets all coefficients of the polynomial to 0.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// we change the degree of P : all coefficients are set to 0
P.SetOrder(4);
P(0) = 2.5; P(2) = -4.1; P(4) = 6.0;

Location :

Share/UnivariatePolynomial.cxx

ResizeOrder

Syntax

  void ResizeOrder(int k);

This method modifies the degree of the polynomial and sets new coefficients of the polynomial to 0. Old coefficients are kept.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// we change the degree of P and we keep old coefficients
P.SetOrder(4);
// P = 4x^4 + 3 x^3 + 2x^2 + 5x - 8
P(4) = 4.0; P(3) = 3.0;

Location :

Share/UnivariatePolynomial.cxx

Zero

Syntax

  void Zero();

This method sets all coefficients to 0.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// P = 0
// P is not resized, coefficients are just set to 0
P.Zero();

Location :

vector/VectorInline.cxx

Evaluate

Syntax

  T Evaluate(T x);

This method evaluates the polynomial at a given point x with Horner's algorithm.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// we compute y = P(x)
Real_wp x = 5.0;
Real_wp y = P.Evaluate(x);

Location :

Share/UnivariatePolynomial.cxx

Constructor for LegendrePolynomial

Syntax

LegendrePolynomial(int n);

The constructor takes the maximum number of Legendre polynomials to evaluate.

Example :

// coefficients of the recurrence are computed
int n = 10;
LegendrePolynomial<Real_wp> P(n);

// now we can evaluate P_0(x), P_1(x), ..., P_{n-1}(x)
// where P_k is the k-th Legendre polynomial
VectReal_wp EvalP;
Real_wp x = 0.4; 
P.EvaluatePn(n, x, EvalP);

Location :

Share/UnivariatePolynomial.cxx

EvaluateDerivative

Syntax

void EvaluateDerivative(int n, T cos_teta, T sin_teta, Vector& Pn0, Vector& Pn1, Vector& dPn1);

This method evaluates P_n(cos θ) where P_n are Legendre polynomials. It also evaluates sin θ P'_n(cos θ) and the derivative of this expression.

Example :

// coefficients of the recurrence are computed
int n = 10;
LegendrePolynomial<Real_wp> P(n);

// P_k(cos teta) are evaluated for k = 0..n-1
// P_k^1 = sin teta P'_k(cos teta) are also evaluated
// and the derivatives of P_k^1 = cos teta P'_k(cos_teta) - sin^2 teta P''_k(cos_teta)
// where P_k is the k-th Legendre polynomial
VectReal_wp EvalP, EvalP1, dP1;
Real_wp teta = pi_wp/4, cos_teta = cos(teta), sin_teta = sin(teta); 
P.EvaluateDerivative(n, cos_teta, sin_teta, EvalP, EvalP1, dP1);

Location :

Share/UnivariatePolynomial.cxx

EvaluatePn

Syntax

void EvaluatePn(int n, T x, Vector& EvalP);

This method evaluates the Legendre polynomials P₀, P₁, ..., P_n-1 at a given point.

Example :

// coefficients of the recurrence are computed
int n = 10;
LegendrePolynomial<Real_wp> P(n);

// now we can evaluate P_0(x), P_1(x), ..., P_{n-1}(x)
// where P_k is the k-th Legendre polynomial
VectReal_wp EvalP;
Real_wp x = 0.4; 
P.EvaluatePn(n, x, EvalP);

Location :

Share/UnivariatePolynomial.cxx

EvaluateDerive

Syntax

void EvaluateDerive(int n, T x, Vector& EvalP, Vector& dP);

This method evaluates the Legendre polynomials P₀, P₁, ..., P_n-1 at a given point and the derivatives.

Example :

// coefficients of the recurrence are computed
int n = 10;
LegendrePolynomial<Real_wp> P(n);

// now we can evaluate P_0(x), P_1(x), ..., P_{n-1}(x)
// where P_k is the k-th Legendre polynomial
// and the derivatives P'_0(x), P'_1(x), ..., P'_{n-1}(x)
VectReal_wp EvalP, derP;
Real_wp x = 0.4; 
P.EvaluateDerive(n, x, EvalP, derP);

Location :

Share/UnivariatePolynomial.cxx

Init

Syntax

void Init(int Lmax)

This method comptes coefficients involved in reccurence relation used to evaluate normalized associated Legendre polynomials. They are defined as:

where are the usual associated Legendre polynomials satisfying the recurrence:

The normalization factor is designed such that the L₂ norm of polynomials over the sphere is equal to one.

Example :

// coefficients of the recurrence are computed
// we need to know the maximal value of l
int Lmax = 100;
AssociatedLegendrePolynomial<Real_wp> P(Lmax);

// now we can evaluate \bar{P}_l^m(cos teta)
// these coefficients are stored in EvalP(l)(m)
// and computed for l = 0..Lmax
Vector<VectReal_wp> EvalP;
Real_wp teta = pi_wp/4;
P.EvaluatePnm(Lmax, Lmax, teta, EvalP);

Location :

Share/UnivariatePolynomial.cxx

EvaluatePnm

Syntax

void EvaluatePnm(int Lmax, int Mmax, T teta, Vector<Vector>& EvalP);

This methods evaluates normalized associated Legendre polynomials. They are defined as:

where are the usual associated Legendre polynomials satisfying the recurrence:

The normalization factor is designed such that the L₂ norm of polynomials over the sphere is equal to one. The polynomials are evaluated for l=0..Lmax and for m=0..min(l, Mmax).

Example :

// coefficients of the recurrence are computed
// we need to know the maximal value of l
int Lmax = 100;
AssociatedLegendrePolynomial<Real_wp> P(Lmax);

// now we can evaluate \bar{P}_l^m(cos teta)
// these coefficients are stored in EvalP(l)(m)
// and computed for l = 0..Lmax
Vector<VectReal_wp> EvalP;
Real_wp teta = pi_wp/4;
// The second argument can be lower than Lmax, if you do not want to compute all polynomials m=0..l
P.EvaluatePnm(Lmax, Lmax, teta, EvalP);

Location :

Share/UnivariatePolynomial.cxx

Pow

Syntax

  UnivariatePolynomial Pow(const UnivariatePolynomial& P, int n);

This functions returns the power of a polynomial.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P, Q;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// if you want to compute Q = P^4
Q = Pow(P, 4);

Location :

Share/UnivariatePolynomial.cxx

DerivatePolynomial

Syntax

  void DerivatePolynomial(const UnivariatePolynomial& P, UnivariatePolynomial& Q);

This functions computes the derivative of a polynomial.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P, Q;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// if you want to compute Q = P' = 4x + 5 here
DerivatePolynomial(P, Q);

Location :

Share/UnivariatePolynomial.cxx

IntegratePolynomial

Syntax

  void IntegratePolynomial(const UnivariatePolynomial& P, UnivariatePolynomial& Q);

This functions computes the primitive of a polynomial that cancels for x=0.

Example :

// P is equal to 0
UnivariatePolynomial<Real_wp> P, Q;

// P = 2x^2 + 5x - 8
UnivariatePolynomial<Real_wp> P(2);
P(0) = -8.0; P(1) = 5.0; P(2) = 2.0;

// if you want to compute Q = \int P = 2/3 x^3 + 5/2 x^2 - 8x here
IntegratePolynomial(P, Q);

Location :

Share/UnivariatePolynomial.cxx

ComputeLegendrePolynome

Syntax

  void ComputeLegendrePolynome(int n0, int n, Vector<UnivariatePolynomial>& P)

This functions computes the n first Legendre polynomials. This function should not be used since coefficients of Legendre polynomials blow up. The class LegendrePolynomial should be used instead to evaluate these polynomials.

Example :

// if you want to know the coefficients of Legendre polynomials
// P_0, P_1, ... P_{n-1}
Vector<UnivariatePolynomial<Real_wp> > P;

int n = 10;
ComputeLegendrePolynome(0, n, P);

Location :

Share/UnivariatePolynomial.cxx

Add

Syntax

  void Add(T alpha, const UnivariatePolynomial& P, UnivariatePolynomial& Q);

This functions adds a polynomial to another one.

Example :

UnivariatePolynomial<Real_wp> P, Q;

// if you want to update Q
// Q = Q + 2.5*P
Add(Real_wp(2.5), P, Q);

Location :

Share/UnivariatePolynomial.cxx

Mlt

Syntax

  void Mlt(T alpha, UnivariatePolynomial& Q);

This functions multiplies a polynomial by a coefficient.

Example :

UnivariatePolynomial<Real_wp> P;

// if you want to multiply P by a scalar
// P = P*4.2
Mlt(Real_wp(2.5), P);

Location :

Share/UnivariatePolynomial.cxx

MltAdd

Syntax

  void MltAdd(T alpha, const UnivariatePolynomial& P, UnivariatePolynomial& Q, T beta, UnivariatePolynomial& R);

This functions computes the product of two polynomials P and Q and adds the result to R. R is replaced by beta R + alpha P Q.

Example :

UnivariatePolynomial<Real_wp> P, Q, R;

// if you want to update R
// R = beta R + alpha P Q 
Real_wp alpha(2), beta(3);
MltAdd(alpha, P, Q, beta, R);

Location :

Share/UnivariatePolynomial.cxx

SolvePolynomialEquation

Syntax

   void SolvePolynomialEquation(const UnivariatePolynomial& P, Vector& R, Vector& Rimag);
   void SolvePolynomialEquation(const UnivariatePolynomial<complex>& P, Vector<complex>& R);

This functions computes the roots of the polynomials P, the roots are R + i Rimag. If the polynomial is complex, only R should be given. This functions drops zero-valued roots.

Example :

UnivariatePolynomial<Real_wp> P(3);

// P = 2 x^3 + 3 x^2 + 5 x
P(3) = 2.0; P(2) = 3.0; P(1) = 5.0;

// if you want to find all the roots of P
// here P has 0 as simple root, this root will not be computed
VectReal_wp R, Rimag;
SolvePolynomialEquation(P, R, Rimag);

// R stores real part of roots and Rimag imaginary parts

Location :

Share/UnivariatePolynomial.cxx

Constructors for MultivariatePolynomial

Syntax

  MultivariatePolynomial();
  MultivariatePolynomial(int dimension, int n);

The default constructor initializes a null polynom of dimension 0 (no variable) and degree 0. The second constructor takes as arguments the dimension (i.e. the number of variables) and the degree. The degree is the total degree of the polynomial, i.e. the sum of degrees for each variable. By default all the coefficients are set to 0.

Example :

MultivariatePolynomial<Real_wp> Q;

// here P = 0, you do not need to call Zero()
MultivariatePolynomial<Real_wp> P(2, 3);

// P = 2 x^3 + 4 x^2 y - x y^2 + 5 y^3 + 3 x^2  - 4 x y - 3.2 y^2 + 5 x - 2.5 y + 8
P(3, 0) = 2.0; P(2, 1) = 4.0; P(1, 2) = -1.0; P(0, 3) = 5.0;
P(2, 0) = 3.0; P(1, 1) = -4.0; P(0, 2) = -3.2;
P(1, 0) = 5.0; P(0, 1) = -2.5;
P(0, 0) = 8.0;