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 >
T	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.