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

[Ralph Silva]

Hi Salman,

you have a small error in your code.

You have this     printf("\t %.5f, %.5f\n", roots[i], roots[2*i+1]);

But it should be  printf("\t %.5f, %.5f\n", roots[2*i], roots[2*i+1]);

Now you get the exact results.

Please, see the manual for further information.
(
http://www.gnu.org/software/gsl/manual/html_node/Roots-of-Polynomials-Examples.html)

Best regards,
Ralph.


On Dec 12, 2007 5:23 AM, salman <address@hidden> wrote:

> 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.
> >
>
>
> _______________________________________________
> Bug-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/bug-gsl
>