Re: [Bug-gsl] gsl_sf_bessel_jl bug

[Koichi Takahashi]

Hi Brian,

I tried this on some different configurations,
with gsl-1.8, 1.9, 1.10
compiled with gcc-4.1.2, 4.1.3
on x86 and x86_64.

I constantly get the error on all x86_64 environments,
and the correct result on all x86 machines.

Now I'm inclined to think that this is caused by
the difference between 387/SSE FP units.

Here is the test program I used.

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

main()
{

//    double a = gsl_sf_bessel_jl( 30,  3875.6 ); // OK
//    double a = gsl_sf_bessel_jl( 30,  3875.6138424501978 ); // Error
//    double a = gsl_sf_bessel_jl( 30,  3875.62 ); // OK
      double a = gsl_sf_bessel_jl( 30,  3878.62 ); // Error
//    double a = gsl_sf_bessel_jl( 30,  4000 ); // OK

      printf("%g\n",a);
}

Here I have only l=30 but I get the same error for different
orders, for example, l=21, x=3810.7005297824849.


This is a show-stopper for the simulator I'm developing, so
let me know if there is anything else I could do for you.

thanks!
Koichi




> At Tue, 16 Oct 2007 01:45:59 -0700,
> Koichi Takahashi wrote:
> > gsl_sf_bessel_jl() crashes with at least one specific set of argument
> > values.
> > double a = gsl_sf_bessel_jl( 30,   3875.6138424501978 );
> > fails in gsl_sf_bessel_J_CF1 () and invokes gsl_error().
> > Slightly different values, like
> > gsl_sf_bessel_jl( 30,   3875.610000000000 );
> > or
> > gsl_sf_bessel_jl( 29,   3875.6138424501978 );
> > seem to return correct results.
>
> Hello
>
> Thanks for the email.  Please could you send an example program which
> shows the error and details of the version of GSL and your compiler
> and operating system (the problem does not occur with gsl-1.10 on x86
> with gcc).