Solving Polynomials
Polynomials are everywhere in computer graphics and engineering applications. This solution explains how to use the functions and classes in the cyPolynomial code release that provide high-performance solutions to polynomial equations of degrees 2 to 10+.
High-Performance Polynomial Root Finding
Given a polynomial function of x, we are often interested in finding its real roots,
the values of x that would make the polynomial function equal to zero.
Quadratic (degree 2) polynomials have a well-known and highly-efficient solution.
A common misconception is that solving polynomials of degree 3+ would be too expensive.
Indeed, they are significantly more expensive than quadratics, but
a
high-performance numerical solution to polynomials of degree up to 10+ exists,
as described in the paper I presented at the High-Performance Graphics conference in 2022:
-
Cem Yuksel. 2022. High-Performance Polynomial Root Finding for Graphics.
Proc. ACM Comput. Graph. Interact. Tech. 5, 3, Article 7 (July 2022), 15 pages.
If you are interested in understanding how it works, I would highly recommend that you read the paper and/or watch the talks I gave on this topic. You can find them all here: http://www.cemyuksel.com/?x=polynomials
Below, I explain how to use the C++ source code I developed for solving polynomials in my cyPolynomial code release. I do not describe how the polynomial solver works here. I will only discuss what you need to know to use this source code in your projects.
Common C++ Template Parameters
The cyPolynomial code release includes different polynomial root finding functions, sharing the following C++ template parameters:
template <int N, // polynomial's degree
typename ftype, // floating-point representation (i.e. float or double)
bool boundError, // = false
typename RootFinder, // = RootFinderNewton
typename RootCallback // the function to be called when a root is found
>
N is the degree of the polynomial.
ftype determines the floating-point format.
It must be float or double.
boundError is set to false by default, which is recommended.
When it is set to true the root finder performs additional operations
to ensure that the given error threshold bounds the error.
In practice, however, this is not needed, except for special cases.
See the discussion below on minimizing the error for more information.
RootFinder is the class that performs numerical root finding,
which is needed for polynomials of degree 3+.
The default is the RootFinderNewton class.
Unless you have a custom numerical root finding class, you should use this default class.
See the discussion below on custom numerical root finder
for more information on how you could implement your custom numerical root finder.
RootCallback is a function that is called each time a root is found.
This template parameter is only used by the PolynomialForEachRoot function, explained below.
Root Finding Methods
A polynomial of degree N can have up to N roots,
but, depending on your application, you may not need to find all of them.
Therefore, the cyPolynomial code release
includes different root finding functions summarized below.
These functions take an array of polynomial coefficients, coef.
Note that a polynomial of degree N has N+1 coefficients.
All functions expect the coefficients to be in the order of increasing degrees,
starting with the lowest-degree (i.e. constant) term.
For example, a quadratic (degree 2) polynomial of x can be computed as
coef[0] + x*coef[1] + x*x*coef[2].
All of these functions come with two variants: bounded and unbounded.
The bounded variants take parameters xMin and xMax
that determine the interval bounds for the roots you are interested in.
This prevents unnecessarily computing the roots out of this interval.
The unbounded variant finds roots in the unbounded interval from -∞ to ∞, which is computationally more expensive.
These functions also take an error threshold parameter xError.
This should be the maximum error your application can tolerate.
If you are not sure what to use, simply use the default value,
which provides a good balance between performance and accuracy.
For more information, see the Performance vs. Accuracy/Error
section below.
Find All Roots
PolynomialRoots finds all roots and returns the number of roots found.
This function is recommended if you know that you will need all roots.
int PolynomialRoots( ftype roots[N], // roots to be returned
ftype const coef[N+1], // polynomial's coefficients
ftype xMin, // minimum x value for the interval
ftype xMax, // maxinum x value for the interval
ftype xError // the error threshold
);
int PolynomialRoots( ftype roots[N], // roots to be returned
ftype const coef[N+1], // polynomial's coefficients
ftype xError // the error threshold
);
Find the First Root
If you are interested in the very first root only,
you can use PolynomialFirstRoot.
This function returns true
if a root is found and false when there is no root.
bool PolynomialFirstRoot( ftype roots[N], // roots to be returned
ftype const coef[N+1], // polynomial's coefficients
ftype xMin, // minimum x value for the interval
ftype xMax, // maxinum x value for the interval
ftype xError // the error threshold
);
bool PolynomialFirstRoot( ftype roots[N], // roots to be returned
ftype const coef[N+1], // polynomial's coefficients
ftype xError // the error threshold
);
Check If There Is a Root
PolynomialHasRoot only checks if there is a root
without actually computing the root. It returns true
if a root is found and false otherwise.
This function is ideal when you simply need to know if there is a solution,
but do not care about the actual value of the root.
Note that polynomials of odd degrees always have at least one root,
but the root may or may not be within the given interval bounds.
bool PolynomialHasRoot( ftype const coef[N+1], // polynomial's coefficients
ftype xMin, // minimum x value for the interval
ftype xMax, // maxinum x value for the interval
ftype xError // the error threshold
);
bool PolynomialHasRoot( ftype const coef[N+1], // polynomial's coefficients
ftype xError // the error threshold
);
Find Roots One-by-One
Finally, we have PolynomialForEachRoot functions
that take a callback function instead of an array of roots to be returned:
bool PolynomialForEachRoot( RootCallback callback, // called when a root is found
ftype const coef[N+1], // polynomial's coefficients
ftype xMin, // minimum x value for the interval
ftype xMax, // maxinum x value for the interval
ftype xError // the error threshold
);
bool PolynomialForEachRoot( RootCallback callback, // called when a root is found
ftype const coef[N+1], // polynomial's coefficients
ftype xError // the error threshold
);
These functions find the roots one by one, starting with the smallest root,
and call the given callback function.
The callback function must be in the form
bool RootCallback( ftype root );
PolynomialForEachRoot calls this function for each root.
When this callback function returns true, it terminates root finding.
If it returns false, root finding continues.
Therefore, PolynomialForEachRoot is recommended if you may sometimes
but not always need all roots.
When you find the root you need, you can simply return true from your callback function.
Then, root finding for the remaining roots can be skipped.
This is particularly effective for polynomials of degrees 4+.
Specialized Functions for Quadratic and Cubic Polynomials
There are special variants of these functions for quadratic (degree 2) and cubic (degree 3) polynomials. You can call them directly when you are working with quadratic and cubic polynomials. Alternatively, you can simply use the functions above, which automatically call the quadratic and cubic variants if your polynomial is degree 2 or 3.
The Polynomial Class and Examples
You may need to build your polynomial (i.e. compute its coefficients) on the fly.
For that, you can use the Polynomial class.
It provides methods for manipulating polynomials, such as adding or multiplying two polynomials.
It also has member functions for root finding, which directly call the functions above.
Overall this class provides a simpler implementation.
The following example builds a degree 5 polynomial and computes its roots.
cy::Polynomial<double,5> poly = { -0.01, 1.0, -2.0, 3.0, -4.0, 0.5 };
double roots[5];
double intervalMin = 0.0;
double intervalMax = 1.0;
double errorThreshold = 0.0001;
int numRoots = poly.Roots( roots, intervalMin, intervalMax, errorThreshold );
The following example uses the same polynomial to compute up to 2 roots using a lambda function as the callback.
int numRootsFound = 0;
double rootsFound[2];
bool foundTwoRoots = poly.ForEachRoot(
[&]( double root )
{
rootsFound[ numRootsFound++ ] = root;
return numRootsFound == 2;
}
, intervalMin, intervalMax, errorThreshold );
Performance vs. Accuracy/Error
Only quadratic (degree 2) polynomials have an efficient analytical solution for finding their roots. With polynomials of degree 3+, we use numerical iterations to find the roots.
Numerical iterations are common for computing various mathematical functions. The number of iterations determines the accuracy as well as the computation cost. Ideally, you would like to have the fewest number of iterations that would satisfy your accuracy needs to achieve the best performance.
We automatically determine the number of iterations by using an error threshold parameter xError.
If the solver determines that the error is below this threshold, it returns its current best guess
of the root without performing more iterations.
The ideal error threshold you should pick depends on your application and the range of values you are working with.
The polynomial solver typically produces results with significantly lower error than this threshold.
For example, in my experiments with randomly-generated polynomials,
I observed average error values that are a thousand times smaller with single precision
(i.e. float) and a million times smaller with double precision (i.e. double).
Therefore, this error threshold should not be the desired average error,
but it should be the maximum error you can tolerate.
If you are not sure how to pick your error threshold xError value, you can simply use the default values.
Minimizing the Error
If you like, you can set the error threshold xError to zero.
This would significantly increase the computation time (to 2× or more),
but the accuracy would be close to the precision of the floating-point representation.
There is also the template parameter boundError, which is false by default.
If you set boundError to true, some additional computation is performed ensure that the error
bounded by the given error threshold.
This is seldom useful in practice, because the error is almost always bounded by the error threshold,
without this additional computation.
Its main use is when you set xError to zero. This, together with boundError set to true,
ensures that you get the best representation of the root possible by the floating-point representation you use
(i.e. float or double).
Floating-Point Representation
Note that evaluating polynomials and finding their roots can be highly inaccurate for relatively high-order polynomials
if your floating-point representation does not have sufficient precision.
Single precision (i.e. float) can be sufficient for polynomials of degree 3 to 5, if you are OK with low accuracy.
For higher-order polynomials, I would recommend using double precision (i.e. double).
In fact, even double precision can be insufficient for polynomials of degree 20+.
Custom Numerical Root Finder
The default numerical root finder,
RootFinderNewton,
works quite well in general.
I have experimented with multiple different methods for numerical root finding,
and in my experiments RootFinderNewton
almost always outperformed them.
There are, however, special cases. For example, when there are three roots at the same position, the default numerical root finder performs poorly. A simple bisection method can be faster in this case, though it is almost always much slower in every other case.
Nonetheless, your application may benefit from a specialized numerical root finder or you might be interested in testing a new numerical root finding method. Therefore, you might want to implement your own custom numerical root finder.
Whatever your reason is, the root finding methods in the cyPolynomial code release
are implemented in a way that makes it easy to replace the default numerical root finder
with your custom numerical root finder.
All you need to do is to implement the relevant methods of
the RootFinderNewton class.
Below is a complete implementation of the pure bisection method
for numerical root finding within a finite interval. In this case, we only need to implement
the FindClosed method, as shown below.
class RootFinderBisection
{
public:
template <int N, typename ftype, bool boundError>
static inline ftype FindClosed (
ftype const coef[N+1], // polynomial's coefficients
ftype const deriv[N], // its derivative's coefficients (not used here)
ftype x0, // the minimum limit of the interval
ftype x1, // the maximum limit of the interval
ftype y0, // polynomial's value at x0
ftype y1, // polynomial's value at x1
ftype xError // the error threshold
)
{
ftype twoEpsilon = 2*xError;
ftype xbounds[2] = { x0, x1 };
while ( xbounds[1] - xbounds[0] > twoEpsilon ) {
ftype xr = (xbounds[0] + xbounds[1]) * ftype(0.5);
ftype yr = PolynomialEval<N,ftype>( poly, xr );
if ( yr == 0 ) return xr;
int i = IsDifferentSign( y0, yr );
xbounds[i] = xr;
}
return (xbounds[0] + xbounds[1]) * ftype(0.5);
}
};
Conclusion
I have tried to summarize how you could take advantage of the functions and classes in the cyPolynomial code release and use it for your projects. If you are interested in finding out more about how this polynomial solver actually works, I would recommend that you check out my paper cited above and/or watch my videos on this topic. You can find them all here: http://www.cemyuksel.com/?x=polynomials