Root-Finder

Root-Finder is a header-only univariate polynomial solver, which finds/counts all REAL roots of any polynomial inside a given interval. (Please use the up-to-date master branch.)

Feature

0. The solver is a C++11 header-only library, which is highly optimized on the premise of instruction set independence.

1. It only contains one header file "root_finder.hpp".

2. The interface only contains two functions. One is for roots finding while the other one is for roots counting.

3. As for low order polynomials (linear, quadratic, cubic, and quartic polynomials), the solver calculates their closed-form solutions. In this case, the solver only takes less than 0.25 microsecond for Intel i7-8700 CPU.

4. As for high order polynomials (order >= 5), the solver implements 2 different methods to find all roots. The recommended one is Real Roots Isolation Method. The other one is based on Companion Matrix Method.

5. The Real Roots Isolation Method uses geometrical properties of polynomial roots to tighten a given interval. Both Cauchy’s bound as well as Kojima’s bound are adopted. Normally, the latter can be tighter than the former for several magnitude in most cases. Technically, Fujiwara’s bound is always better than Kojima's bound, while Kojima's bound is more numerically friendly and is tight enough.

6. The Real Roots Isolation Method is much faster and much more stable than the Companion Matrix Method. However, due to truncation error of float point number, the former is recommended for at most 32-order polynomials, while the latter is only recommended for at most 20-order polynomials.

7. We perform benchmark example between our Real Roots Isolation Method and TOMS493: Jenkins–Traub Algorithm on two different platforms. The latter one is commonly known as the "RPOLY" algorithm. In the benchmark, all roots of a polynomial are required. For 8-order polynomials, our method is about 5% faster than RPOLY under Intel i7-8700 CPU. In general, our library outperforms the RPOLY algorithm even when all real roots are needed. Moreover, RPOLY is designed to find all roots of a polynomial while ours can handle any given interval. Therefore, when roots in an interval are required, our method performs far better.

8. The library is also capable of efficiently counting the number of roots inside an interval, of which RPOLY is incapable. For 8-order polynomials, our roots counter only takes about 0.2 microsecond under Intel i7-8700 CPU and 0.8 microsecond under Intel i5-5200U CPU. The time consumption for high order polynomial is extremely low as those low order polynomials which have closed-form solutions.

9. Only algorithmic acceleration is explicitly implemented in the code. Furthur instruction-dependent optimization can be done to get better performance.

Interface

Only two functions below are needed.

Roots Finding Function:

std::set<double> RootFinder::solvePolynomial(const Eigen::VectorXd &coeffs, double lbound, double ubound, double tol, bool isolation = true)

Inputs:

coeffs: 
    Eigen VectorXd for coefficients of a polynomial. For example, the polynomial a(n)*x^n+a(n-1)*x^(n-1)+...+a(1)*x+a(0) can be expressed by 
    coeffs=[a(n), a(n-1), ..., a(1), a(0)].

lbound and ubound:
    Open interval (lbound, ubound). If lbound = -INFINITY and ubound = INFINITY, then all real roots can be found by the solver. Note that polynomial cannot be zero at these two boundaries.

tol:
    tolerance for precision, only works when order is higher than 4.

isolation:
    switch for Method used. Default one is Real Roots Isolation Method.

Outputs:

std::set<double> which stores all searched real roots.

Example:

Eigen::VectorXd coeffs(6);
coeffs << 1.0, -2.0, 3.0, -4.0, 5.0, -6.0;
std::set<double> allRoots = RootFinder::solvePolynomial(coeffs, -100.0, 100.0, 0.00001);

Roots Counting Function:

int RootFinder::countRoots(const Eigen::VectorXd &coeffs, double lbound, double ubound)

Inputs:

coeffs: 
    Eigen VectorXd for coefficients of a polynomial. For example, the polynomial a(n)*x^n+a(n-1)*x^(n-1)+...+a(1)*x+a(0) can be expressed by 
    coeffs=[a(n), a(n-1), ..., a(1), a(0)].

lbound and ubound:
    Open interval (lbound, ubound). Note that polynomial cannot be zero at these two boundaries.

Outputs:

The number of distinct roots inside (lbound, ubound).

Example:

Eigen::VectorXd coeffs(6);
coeffs << 1.0, -2.0, 3.0, -4.0, 5.0, -6.0;
int numRoots = RootFinder::countRoots(coeffs, -10, 0.5);

Compile

sudo apt install libeigen3-dev
mkdir build
cd build
cmake ..
make

Note

The "RPOLY" algorithm which we benchmarks againt comes from a modified version by Helen Oleynikova, which originates in http://www.netlib.org/toms/.