Re: [Bug-gsl] Bug in polynomial root finder

[Munagala Ramanath]

The array of coefficients needs to be reversed. Here is the output and modified 
code:

uname -a
Linux clay 2.6.38-8-generic #42-Ubuntu SMP Mon Apr 11 03:31:24 UTC 2011 x86_64 
x86_64 x86_64 GNU/Linux
[This is Linux Mint 11 (katya) 64 bit]

Output [compiled with: gcc -o gsl_poly -O2 -std=c99 gsl_poly.c -lgsl -lgslcblas 
]:
 ./gsl_poly
gsl: zsolve.c:78: ERROR: root solving qr method failed to converge
Default GSL error handler invoked.
Aborted

Modified code:
#include <stdio.h>
#include <gsl/gsl_poly.h>

int main ()
{
    double b[16] =
          { 1, -4, 2, 10, -17, 10, 6, -16, 12, -16, 16, -8, 28, -8, -48, 32 },
        a[16], z[16*2];
    int i;

    for ( int j = 0; j < 16; ++j ) {
        a[j] = b[ 16 - j - 1 ];
    }

    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (16);

    int status = gsl_poly_complex_solve (a, 16, w, z);

    for (i = 0; i<8; i++)
      printf("%.18e %.18e\n", z[2*i], z[2*i+1]);

    gsl_poly_complex_workspace_free (w);
}

Ram


________________________________
From: Brian Gough <address@hidden>
To: Munagala Ramanath <address@hidden>
Cc: "address@hidden" <address@hidden>
Sent: Sunday, November 6, 2011 12:05 PM
Subject: Re: [Bug-gsl] Bug in polynomial root finder

At Sun, 6 Nov 2011 08:38:22 -0800 (PST),
Munagala Ramanath wrote:
> 
> Thanks for the response.
> 
> I have run it on more than 1.6 million polynomials of order 15 and this is
> the only one on which it has problems.
> 
> I've also run it on some polynomials of order 18 and again it had no problems.
> 
> So it seems the order is not the problem, just some corner case with this
> particular polynomial.
> 

Ah, thanks. That is useful information.

Can you confirm what platform you are one and try the program as a
standalone test case in C -- it works fine for me.

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

int main ()
{
    double a[16] = { 1, -4, 2, 10, -17, 10, 6, -16, 12, -16, 16, -8, 28, -8, 
-48, 32 };
    double z[16*2];
    int i;

    gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc (16);

    int status = gsl_poly_complex_solve (a, 16, w, z);

    for (i = 0; i<8; i++)
      printf("%.18e %.18e\n", z[2*i], z[2*i+1]);

    gsl_poly_complex_workspace_free (w);
}



$ ./a.out 
-1.000000000000002442e+00 0.000000000000000000e+00
-5.827430765394314705e-01 7.380899912666477602e-01
-5.827430765394314705e-01 -7.380899912666477602e-01
-7.652876256274119271e-01 0.000000000000000000e+00
-6.219637508703932394e-01 0.000000000000000000e+00
-5.472734432737202948e-02 8.769244425642294116e-01
-5.472734432737202948e-02 -8.769244425642294116e-01
3.245172121152936628e-01 6.791075284395361455e-01
$ $ uname -a
Linux nc 2.6.35-30-generic #61-Ubuntu SMP Tue Oct 11 15:29:15 UTC 2011 i686 
GNU/Linux