ebertolazzi/quarticrootsflocke

C++ rewrite of rpoly.f

Quartic Roots Solver

Polynomial Roots Solver

Port to C++ of the algorithm of Norbert Flocke for polynomial roots up to degree 4. Also Jenkins-Traub real polynomial root finder is ported but is experimental for the moment.

Usage

#include "PolynomialRoots.hh"

// ....

double coeffs[] = { 8, -8, 16,-16, 8,-8 }; // polynomial coeffs

double zeror[5], zeroi[5];
int    info[5];
int    degree = 5;

int ok = PolynomialRoots::roots( coeffs, degree, zeror, zeroi ); // ok < 0 failed
cout << " ok = " << ok << '\n' ;
for ( int i = 0 ; i < degree ; ++i )
  cout << zeror[i] << " + I* " << zeroi[i] << '\n';

To solve quadratic, cubic or quartic use specialized classes

Quadratic qsolve(a,b,c);
qsolve.info(cout);

Cubic csolve(a,b,c,d);
csolve.info(cout);

Quartic q4(a,b,c,d,e);
q4.info(cout);

look at the class definition to see how to access the computed roots.

Repository and Online Documentation

• https://github.com/ebertolazzi/quarticRootsFlocke
• http://ebertolazzi.github.io/quarticRootsFlocke

References

• N.Flocke

  Algorithm 954: An Accurate and Efficient Cubic and Quartic

  Equation Solver for Physical Applications

  ACM TOMS, vol 41, n.4, 2015
• M.A. Jenkins and J.F.Traub

  A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration

  SIAM Journal on Numerical Analysis

  Vol. 7, No. 4 (Dec., 1970), pp. 545-566

Author

Enrico Bertolazzi

Dipartimento di Ingegneria Industriale

Universita` degli Studi di Trento

email: [email protected]

quarticrootsflocke's Issues

alternative root solver

Hi,
as a faster alternative to Flocke's method, which can be used also for finding roots of complex polynomials,
I would like to draw your attention to my algorithm https://github.com/cridemichel/quarticpp ,
whichi is based on:

[1] A.Orellana and C. De Michele ACM Transactions on Mathematical Software, Vol. 46, No. 2, 20 (2020), https://doi.org/10.1145/3386241
[3] C. De Michele, ACM Transactions on Mathematical Software, Vol. 48 Issue 4 Article No. 46 pp 1-3, https://doi.org/10.1145/3564270

best Cristiano

Invalid solutions found for quartic poly

Solving the following poly, it gives two invalid solutions:
p(x) = ax^4 + bx^3 + cx^2 + dx + e
a = 714285.71428571432
b = -100000
c = 13479.924671428573
d = -2.2737367544323206e-13
e = -22.526485131747265

Matlab solutions with 'roots' are:
r1 = 0.064423196657504 + 0.122215870576901i
r2 = 0.064423196657504 - 0.122215870576901i
r3 = 0.046605762356559
r4 = -0.035452155671566

[Feature] Specify absolute and relative tolerances

In many cases it's not required an high precision on roots solutions. For example if the (real) roots represents a time, having an absolute tolerance of milliseconds or microseconds could be enough in real applications.
Actually Newton iterations and bisection stop using a relative tolerance equal to machine precision: std::numeric_limits<valueType>::epsilon()
It's useful to setup the user tolerance (both absolute and relative) in class constructor or with a proper method. In this way the algorithm should be faster.

Create a new release with master's HEAD

I've been attempting to create a Conda package to publish on conda-forge, but as I'm using the last release tag (1.0.0) from this repo (and using releases is a good practice in conda packaging), it led to some issues (see this conda-forge CI log). These issues are fixed by the last changes on CMake files that are on master after the last release commit.

Making a new release would fix these issues.

Invalid solution found for quartic poly

Solving the following poly, it gives invalid solutions for complex roots and one real root:
p(x) = ax^4 + bx^3 + cx^2 + dx + e
a = -5.898376305898030e-09
b = -6.000000000001007e+06
c = 3.664603343001005e-27
d = 0
e = 9.000000000003021e+02

Matlab solutions with 'roots' are:
r1 = -1.017229096421223e+15
r2 = -0.026566464600067 + 0.046014465225927i
r3 = -0.026566464600067 - 0.046014465225927i
r4 = 0.053132929200150

For example the root residual for r4 is:

• -3.765546080103377e-05 with Matlab
• 9.000000000003021e+02 with quarticRootsFlocke

Setting a coeffient to 0 the solution is correct. It seems an ill-conditioned problem.

Wrong solution for sixth order poly

Consider the sixth order polynomial:

double c6 = 20833333.333333332;
double c5 = 0.0000000000000000;
double c4 =  -1083333.3333333372;
double c3 = -7650000.0000000000;
double c2 = 281666.66666666680;
double c1 = 0.0000000000000000;
double c0 = 658458.33333333337;

The following method call return an incorrect solution:

double coeffs[] = {c6, c5, c4, c3, c2, c1, c0};
double re[6], im[6];
int degree = 6;
int ok = PolynomialRoots::roots( coeffs,  degree,  re,  im );

More precisely, the method provided the following roots:

roots[0] = 0.51539085009371632 + 0.0000000000000000i;
roots[1] = -0.19792110566656965 + 0.44658964697441494i;
roots[2] = -0.19792110566656965 + -0.44658964697441494i;
roots[3] = -0.37675566813624506 + 0.51326747398874339i;
roots[4] = -0.37675566813624506 + -0.51326747398874339i;
roots[5] = -0.37675566813624506 + -0.51326747398874339i;

Instead the roots found with matlab for the same poly are:

roots[0] = 0.633962697511912 + 0.00000000000000i
roots[1] = 0.515390850093718 + 0.00000000000000i
roots[2] = -0.376755668136244 + 0.513267473988743i
roots[3] = -0.376755668136244 - 0.513267473988743i
roots[4] = -0.197921105666570 + 0.446589646974415i
roots[5] = -0.197921105666570 - 0.446589646974415i

As it is possible noticed, there is a root that does not coincide with the one found by matlab.
The issue has been debugged with VS2019 using the commit 8e8fbe3###

Neverending loop for quartic poly

Consider the quartic poly:

double c4 = 5.0 * ( -1462.2326183798939 );
double c3 = 4.0 * ( 1868.4944990527565 );
double c2 = 3.0 * ( -901.83171888631921 );
double c1 = 2.0 * ( 206.27445843560130 );
double c0 = 1.0 * ( -23.818259047289210 );

The following method call goes into a neverending loop:

double re[4], im[4];
int n = QuarticGetRoots( c4, c3, c2, c1, c0, re, im );

Debugging the code, there's a division by 0 inside zeroQuarticByNewtonBisection method:
x -= p/dp
with dp = 0.
Changing a little the poly coeffs (i.e. using c4 = -5.0 * 1462.232618379893) it gives the correct results.
The issue has been debugged with VS2013 and VS2019 using the commit bc377b5

alternative quartic solver

Hi,
as a faster alternative to Flocke's method, which can be used also for finding roots of complex polynomials,
I would like to draw your attention to my algorithm https://github.com/cridemichel/quarticpp ,
whichi is based on:

[1] A.Orellana and C. De Michele ACM Transactions on Mathematical Software, Vol. 46, No. 2, 20 (2020), https://doi.org/10.1145/3386241
[3] C. De Michele, ACM Transactions on Mathematical Software, Vol. 48 Issue 4 Article No. 46 pp 1-3, https://doi.org/10.1145/3564270

best Cristiano

Different found roots for Quartic poly after update

Updating from commit fa763d2 to commit c0be672 I get different solutions for the following Quartic poly:

p = [86594.643180459301, -509.64355468067714, 0.74986322275510664, 0.0000000000000000, -1.7386537088520670e-16]

Old roots (all real):
r0 = -1.5226962261856962e-08
r1 = 1.5227119847090433e-08
r2 = 0.0029426967763176593
r3 = 0.0029426967763176593

New roots (2 real, 2 complex):
r0 = 0.0029426967762312241 + i * 2.2554513436451675e-08
r1 = 0.0029426967762312241 - i * 2.2554513436451675e-08
r2 = -1.5226962261856962e-08
r3 = 1.5227119847090433e-08

Matlab solutions match with the old version.

[Hint] Increase performance for low order poly in Jenkins-Traub algorithm

In roots method of Jenkins-Traub the algorithm iterates until poly degree is 3 then it uses Flocke algorithm. Flocke algorithm is available for 4th order poly too so the iterations could stop when N <= 4.
With this modification my benchmarks show up to 12% running time decrease (tested with sixtic poly).

while ( N > 0 ) {
      // Main loop
      // Start the algorithm for one zero
      if ( N < 4 ) { roots3( p, N, zeror+Degree-N, zeroi+Degree-N ); break; }
      // ...