[Help-gsl] Polynomial root finding

[Will Leckie]

Hi,

Great work on the GSL, I love using it in my research.

I'm using the GSL in my research to solve arbitrary polynomials of order up to 
4.  I've written a little wrapper around the GSL to use 
gsl_poly_complex_solve_quadratic() for quadratics and 
gsl_poly_complex_solve_cubic() for cubics, and I use gsl_poly_complex_solve() 
to solve quartics.  

I also implemented my own analytical quartic solver based on Ferrari's method 
(http://mathworld.wolfram.com/QuarticEquation.html) using the GSL complex math 
operations such as gsl_complex_mul etc. to manipulate the complex numbers 
involved in Ferrari's solution.

When comparing the speed of the GSL gsl_poly_complex_solve() solution of a 
quartic to my analytical version, I find that the analytic routine runs just a 
touch faster:

*********** timing GSL version...
GSL version took 1.078 seconds for 100000 runs; averaging 1.078e-05 seconds per 
run.

*********** Analytic version...
Analytic version took 0.797 seconds for 100000 runs; averaging 7.97e-06 seconds 
per run.


(I'm running a 3.2 GHz Pentium IV with 1GB RAM, Cygwin on Windows XP for the 
timing tests, using the cstdlib rand() function to generate "random" 
coefficients between 0.0 and 10.0 for the quartics for each run)

When comparing the "precision" of the two quartic solvers, I find a very small 
(machine-level?) difference between the roots that they find, typically 10^-30, 
using gsl_complex_sub() to calculate the difference between roots.

I also find that sometimes the GSL gsl_poly_complex_solve() routine returns a 
non-zero status int meaning that it didn't find roots or it found 
unstable/overflow/underflow conditions, I presume.  

So, my question is, what is my advantage in using the GSL 
gsl_poly_complex_solve() routine to solve my quartics?  I have a feeling the 
GSL gsl_poly_complex_solve() routine is more robust to large and small 
coefficients in the polynomial, because I haven't taken too much care in 
looking out for those in my analytic code.  I'm sure the GSL developers looked 
at Ferrari's method when considering whether or not to implement a 
gsl_poly_complex_solve_quartic() routine, since an analytic solution is 
acheivable, and I'm wondering why they didn't implement one?  Is there some 
obvious flaw when implementing Ferrari's method computationally?  What are the 
other advantages of solving quartics using the iterative method in 
gsl_poly_complex_solve() instead?  

If you'd like to see any of the code I've used, either for Ferrari's method, or 
for testing, I'd be happy to share it.

Thanks again for the great library.

Cheers,
Will