[Bug-gsl] bug in polynome function

[edward hage]

The following program delivers different output dependant if in function 
roots the results are shown with printf.
I can not get the results from z to main without the printf statement. 
Clearly this is a bug.


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

double *roots (double p3, double p2, double p1, double p0)
{
 double a[4] = { p0, p1, p2, p3 };
 double z[6];
 double *zp;
 zp = z;

 gsl_poly_complex_workspace * w
     = gsl_poly_complex_workspace_alloc (4);

 gsl_poly_complex_solve (a, 4, w, z);

printf("%+.18f\n", zp[5]); // toggle this on and off

 gsl_poly_complex_workspace_free (w);

 return zp;
}

int main (void)
{
 int i;
 double *z;

 /* p3 p2 p1 p0 for polynome */
 z = roots(1, 2, -0.23,1);

 for (i = 0; i < 3; i++)
   {
     printf ("z%d = %+.18f %+.18f*i\n",
             i, z[2*i], z[2*i+1]);
   }

 return 0;
}


Returns when printf toggled of in roots, this is FALSE:
z0 = +0.145465665064093808 +0.644471441876203199*i
z1 = +0.145465665064093808 +2.151574274584890922*i
z2 = +0.000000000000000000 +2.086480617521523584*i

Returns when printf toggled on in roots, this is GOOD:
+0.000000000000000000
z0 = +0.145465665064093808 +0.644471441876203199*i
z1 = +0.145465665064093808 -0.644471441876203199*i
z2 = -2.290931330128187060 +0.000000000000000000*i


I use:  GNU Scientific Library  version 1.4 -35 as delivered with SuSE 9.1

I would appreciate a repsonse to this problem.
Thank you in advance,
Edward Hage