Files Class List
cyPolynomial.h File Reference

Detailed Description

Polynomial operations and solver class and functions.

Author
Cem Yuksel

This file includes the implementation of the numerical polynomial solver described in

Cem Yuksel. 2022. High-Performance Polynomial Root Finding for Graphics. Proc. ACM Comput. Graph. Interact. Tech. 5, 3, Article 7 (July 2022), 15 pages.

http://www.cemyuksel.com/?x=polynomials

The polynomial functions take an array of polynomial coefficients in the order of increasing degrees.

#include "cyCore.h"

Classes

class  Polynomial< ftype, N >
 A general-purpose polynomial class. More...
 
class  RootFinderNewton
 Numerical root finder using safe Newton iterations. More...
 

Functions

General Polynomial Functions

These functions can be used for evaluating polynomials and their derivatives, computing their derivatives, deflating them, and inflating them.

template<int N, typename ftype >
ftype PolynomialEval (ftype const coef[N+1], ftype x)
 Evaluates the given polynomial of degree N at x.
 
template<int N, typename ftype >
ftype PolynomialEvalWithDeriv (ftype &derivativeValue, ftype const coef[N+1], ftype x)
 Evaluates the given polynomial and its derivative at x.
 
template<int N, typename ftype >
void PolynomialDerivative (ftype deriv[N], ftype const coef[N+1])
 Computes the polynomial's derivative and sets the coefficients of the deriv array.
 
template<int N, typename ftype >
void PolynomialDeflate (ftype defPoly[N], ftype const coef[N+1], ftype root)
 Deflates the given polynomial using one of its known roots.
 
template<int N, typename ftype >
void PolynomialInflate (ftype infPoly[N+2], ftype const coef[N+1], ftype root)
 Inflates the given polynomial using the given root.
 
Polynomial Root Finding Functions

These functions find polynomial roots. They can be limited to a finite bounds between xMin and xMax. Functions for degree 3 and higher polynomials use numerical root finding. By default, they use Newton iterations defined in RootFinderNewton as their numerical root finding method, but they can also be used with a custom class that provides the same interface.

The given xError parameter is passed on to the numerical root finder. The default value is 0, which aims to find the root up to numerical precision. This might be too slow. Therefore, for a high-performance implementation, it is recommended to use a non-zero error threshold for xError.

If boundError is false (the default value), RootFinderNewton does not guarantee satisfying the error bound xError, but it almost always does. Keeping boundError false is recommended for improved performance.

template<typename ftype >
ftype PolynomialDefaultError ()
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<>
float PolynomialDefaultError< float > ()
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<>
double PolynomialDefaultError< double > ()
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
int PolynomialRoots (ftype roots[N], ftype const coef[N+1], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
int CubicRoots (ftype roots[3], ftype const coef[4], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
int QuadraticRoots (ftype roots[2], ftype const coef[3], ftype xMin, ftype xMax)
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
int LinearRoot (ftype &root, ftype const coef[2], ftype xMin, ftype xMax)
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
int PolynomialRoots (ftype roots[N], ftype const coef[N+1], ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial and returns the number of roots found.
 
template<typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
int CubicRoots (ftype roots[3], ftype const coef[4], ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
int QuadraticRoots (ftype roots[2], ftype const coef[3])
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
int LinearRoot (ftype &root, ftype const coef[2])
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
bool PolynomialFirstRoot (ftype &root, ftype const coef[N+1], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Finds the first root of the given polynomial between xMin and xMax and returns true if a root is found.
 
template<typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
bool CubicFirstRoot (ftype &root, ftype const coef[4], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
bool QuadraticFirstRoot (ftype &root, ftype const coef[3], ftype xMin, ftype xMax)
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
bool PolynomialFirstRoot (ftype &root, ftype const coef[N+1], ftype xError=PolynomialDefaultError< ftype >())
 Finds the first root of the given polynomial and returns true if a root is found.
 
template<typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
bool CubicFirstRoot (ftype &root, ftype const coef[4], ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
bool QuadraticFirstRoot (ftype &root, ftype const coef[3])
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
bool PolynomialHasRoot (ftype const coef[N+1], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Returns true if the given polynomial has a root between xMin and xMax.
 
template<typename ftype >
bool CubicHasRoot (ftype const coef[4], ftype xMin, ftype xMax)
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
bool QuadraticHasRoot (ftype const coef[3], ftype xMin, ftype xMax)
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton>
bool PolynomialHasRoot (ftype const coef[N+1], ftype xError=PolynomialDefaultError< ftype >())
 Returns true if the given polynomial has a root.
 
template<typename ftype >
bool CubicHasRoot (ftype const coef[4])
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype >
bool QuadraticHasRoot (ftype const coef[3])
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton, typename RootCallback >
bool PolynomialForEachRoot (RootCallback callback, ftype const coef[N+1], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Calls the given callback function for each root of the given polynomial between xMin and xMax in increasing order. The callback function should have the following form:
 
template<typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton, typename RootCallback >
bool CubicForEachRoot (RootCallback callback, ftype const coef[4], ftype xMin, ftype xMax, ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype , typename RootCallback >
bool QuadraticForEachRoot (RootCallback callback, ftype const coef[3], ftype xMin, ftype xMax)
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton, typename RootCallback >
bool PolynomialForEachRoot (RootCallback callback, ftype const coef[N+1], ftype xError=PolynomialDefaultError< ftype >())
 Calls the given callback function for each root of the given polynomial in increasing order. The callback function should have the following form:
 
template<typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton, typename RootCallback >
bool CubicForEachRoot (RootCallback callback, ftype const coef[4], ftype xError=PolynomialDefaultError< ftype >())
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
template<typename ftype , typename RootCallback >
bool QuadraticForEachRoot (RootCallback callback, ftype const coef[3])
 Finds the roots of the given polynomial between xMin and xMax and returns the number of roots found.
 
Support Functions (Internal)
template<typename T , typename S >
MultSign (T v, S sign)
 Multiplies the given value with the given sign.
 
template<typename T , typename S >
bool IsDifferentSign (T a, S b)
 Returns true if the sign bits are different.
 

Function Documentation

◆ PolynomialEval()

template<int N, typename ftype >
ftype PolynomialEval ( ftype const coef[N+1],
ftype x )

Evaluates the given polynomial of degree N at x.

The coefficients in the coef array must be in the order of increasing degrees.

◆ PolynomialEvalWithDeriv()

template<int N, typename ftype >
ftype PolynomialEvalWithDeriv ( ftype & derivativeValue,
ftype const coef[N+1],
ftype x )

Evaluates the given polynomial and its derivative at x.

This function does not require computing the derivative of the polynomial, but it is slower than evaluating the polynomial and its precomputed derivative separately. Therefore, it is not recommended when the polynomial and it derivative are computed repeatedly. The coefficients in the coef array must be in the order of increasing degrees.

◆ PolynomialDerivative()

template<int N, typename ftype >
void PolynomialDerivative ( ftype deriv[N],
ftype const coef[N+1] )

Computes the polynomial's derivative and sets the coefficients of the deriv array.

Note that a degree N polynomial's derivative is degree N-1 and has N coefficients. The coefficients are in the order of increasing degrees.

◆ PolynomialDeflate()

template<int N, typename ftype >
void PolynomialDeflate ( ftype defPoly[N],
ftype const coef[N+1],
ftype root )

Deflates the given polynomial using one of its known roots.

Stores the coefficients of the deflated polynomial in the defPoly array. Let f(x) represent the given polynomial. This computes the deflated polynomial g(x) of a lower degree such that

f(x) = (x - root) * g(x).

The given root must be a valid root of the given polynomial. Note that the deflated polynomial has degree N-1 with N coefficients. The coefficients are in the order of increasing degrees.

◆ PolynomialInflate()

template<int N, typename ftype >
void PolynomialInflate ( ftype infPoly[N+2],
ftype const coef[N+1],
ftype root )

Inflates the given polynomial using the given root.

Stores the coefficients of the inflated polynomial in the infPoly array. Let f(x) represent the given polynomial. This computes the inflated polynomial g(x) of a higher degree such that

g(x) = (x - root) * f(x).

Note that the inflated polynomial has degree N+1 with N+2 coefficients. The coefficients are in the order of increasing degrees.

◆ PolynomialForEachRoot() [1/2]

template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton, typename RootCallback >
bool PolynomialForEachRoot ( RootCallback callback,
ftype const coef[N+1],
ftype xMin,
ftype xMax,
ftype xError = PolynomialDefaultError<ftype>() )

Calls the given callback function for each root of the given polynomial between xMin and xMax in increasing order. The callback function should have the following form:

bool RootCallback( ftype root );

If the callback function returns true, root finding is terminated without finding any additional roots. If the callback function returns false, root finding continues. Returns true if root finding is terminated by the callback function. Otherwise, returns false.

◆ PolynomialForEachRoot() [2/2]

template<int N, typename ftype , bool boundError = false, typename RootFinder = RootFinderNewton, typename RootCallback >
bool PolynomialForEachRoot ( RootCallback callback,
ftype const coef[N+1],
ftype xError = PolynomialDefaultError<ftype>() )

Calls the given callback function for each root of the given polynomial in increasing order. The callback function should have the following form:

bool RootCallback( ftype root );

If the callback function returns true, root finding is terminated without finding any additional roots. If the callback function returns false, root finding continues. Returns true if root finding is terminated by the callback function. Otherwise, returns false.