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

[salman]

Hi Brian,

Here is a program that reproduces the results:

#include <stdio.h>
#include <gsl_poly.h>

int main(int argc, char **argv) {
    int i;
    double coeffs[4]={1, -5, -25, 125};
    double roots[6];
    gsl_poly_complex_workspace *u = gsl_poly_complex_workspace_alloc(4);
    gsl_poly_complex_solve(coeffs, 4, u, roots);
    gsl_poly_complex_workspace_free(u);
    for(i=0; i<3; i++)
        printf("\t %.5f, %.5f\n", roots[i], roots[2*i+1]);
    printf("\t------------\n");
    coeffs[0]=1; coeffs[1]=-1; coeffs[2]=-5; coeffs[3]=125;
    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc(4);
    gsl_poly_complex_solve(coeffs, 4, w, roots);
    gsl_poly_complex_workspace_free(w);
    for(i=0; i<3; i++)
        printf("\t %.5f, %.5f\n", roots[i], roots[2*i+1]);
}

The output of the program is:

         0.20000, 0.00000
         0.00000, -0.00000
         0.20000, 0.00000
        ------------
         -0.20000, 0.00000
         0.00000, 0.16000
         0.12000, -0.16000

But the first polynomial's roots are 0.2, 0.2, and -0.2. The second 
polynomial has roots 0.2, 0.12+0.16*i, and 0.12-0.16*i.

I'm running Ubuntu with kernel 2.6.20-16-386 Version 2 and GSL 1.8-3. I 
compiled using g++ 4.1.2-1. Is there any other info you need?

Brian Gough wrote:
> At Mon, 10 Dec 2007 16:41:31 -0600 (CST),
> salman wrote:
>
> > 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]);
>
> Thanks for your email.  Please can you send a complete example program
> so we can reproduce the problem.  Also, details of which version of
> GSL you are using and operating system etc.