[Bug-gsl] Bug in gsl_poly_complex_solve() producing inaccurate roots of polynomials

[salman]

Hi,

I am using gsl_poly_complex_solve() to determine the roots of some basic 
polynomials. The results I am getting are fairly off from the correct 
roots. Here is the relevant section of the code (deg_L is the degree of 
the polynomial I am considering):

  gsl_poly_complex_workspace *u=gsl_poly_complex_workspace_alloc(deg_L+1);
  gsl_poly_complex_solve(coeffs, deg_L+1, u, roots);
  gsl_poly_complex_workspace_free(u);
  for(i=0; i<deg_L; i++)
    printf("\t %.5f, %.5f\n", roots[i], roots[2*i+1]);

Here is the output from one particular polynomial:

f(x) = 1 -5*x -25*x^2 +125*x^3

         0.20000, 0.00000
         0.00000, -0.00000
         0.20000, 0.00000

The roots of this polynomial should be 0.2, 0.2, -0.2. Here's another 
example without multiple roots:

f(x) = 1 -1*x -5*x^2 +125*x^3

         -0.20000, 0.00000
         0.00000, 0.16000
         0.12000, -0.16000

The roots of this polynomial should be -0.2, 0.12+0.16*i, 0.12-0.16*i.

I am not sure what could be causing the error. Being off by as much as 0.2 
is fairly significant. Perhaps I have a bug in my code, but that is 
doubtful -- it does nothing more than read in these polynomials and then 
tries to determine their roots via gsl_poly_complex_solve(). I checked my 
roots by hand and using PARI/GP, so I'm positive they're accurate (and 
clearly, 0 is not a root of the first polynomial). This leads me to think 
that there is a bug somewhere in the routine, or perhaps the precision 
settings need to be adjusted.

If I am doing something wrong, please let me know.

Thanks,
Salman