[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] Bug in polynomial root finder
From: |
Munagala Ramanath |
Subject: |
Re: [Bug-gsl] Bug in polynomial root finder |
Date: |
Sun, 6 Nov 2011 16:54:33 -0800 (PST) |
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