Re: [Bug-gsl] bug in polynome function

[eknecronzontas]

No. It's not a bug. You're returning a pointer to locally allocated storage. 
If you want the function, roots(), to return all the roots, then you'll have
to use malloc to allocate storage on the heap, or make the roots pointer
arguments to your function. 

(BTW, My advice would be to upgrade to a newer version of GSL as well.)

----- Original Message ----
From: edward hage <address@hidden>
To: address@hidden
Sent: Thursday, April 12, 2007 8:37:32 AM
Subject: [Bug-gsl] bug in polynome function

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




_______________________________________________
Bug-gsl mailing list
address@hidden
http://lists.gnu.org/mailman/listinfo/bug-gsl