ebertolazzi / quarticrootsflocke Goto Github PK
View Code? Open in Web Editor NEWC++ rewrite of rpoly.f
License: Other
C++ rewrite of rpoly.f
License: Other
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###
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:
Setting a coeffient to 0 the solution is correct. It seems an ill-conditioned problem.
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
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.
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
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; }
// ...
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
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.
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.
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
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.