Re: [Bug-gsl] incorrect values by bessel_steed_array

[Hu Zhan]

Hi,

I tried make check. The entire library passed the test, but the test for
gsl_sf_bessel_jl_steed_array was hard-coded with LMAX = 50.

On an intel p4 with pre-compiled gsl-1.5 for redhat, I got the correct
results. However, when I compile the test program with the source code
of gsl_sf_bessel_jl_steed_array (gcc -Wall steed.c -lm; see the attached
program) instead of linking to the library, the results are still
incorrect:

      0               1
    0.1               0
    0.2               0
    0.3               0
    0.4               0
    0.5               0
    0.6               0
    0.7               0
    0.8        0.896695

There is one line in gsl_sf_bessel_jl_steed_array:

    W = x_inv / sqrt(FP*FP + F*F);

When LMAX is large, FP and F can be larger than 1e160, and cause W = 0.
The option -mfpmath=387 (90-bit arithmetics internally) can make it work
on the opteron, but not pentium.

The problem is mitigated when I change the line to

    if (FP > 1e100) {
      B = F / FP;
      W = x_inv / (FP * sqrt(1. + B * B));
    }
    else if (F > 1e100) {
      B = FP / F;
      W = x_inv / (F * sqrt(1. + B * B));
    }
    else
      W = x_inv / sqrt(FP*FP + F*F);

Results:

      0               1
    0.1        0.998334
    0.2        0.993347
    0.3        0.985067
    0.4        0.973546
    0.5        0.958851
    0.6        0.941071
    0.7        0.920311
    0.8        0.896695
 


On Tue, 2005-02-01 at 12:30, Brian Gough wrote:
> Hu Zhan writes:
>  > When I use gsl_sf_bessel_jl_steed_array(LMAX, x, wrk) with a large LMAX,
>  > e.g. LMAX = 80, the function returns wrong values at small l and x. My
>  > system is Debian Linux with Opteron CPU.
> 
> Hi,
> 
> Thanks for the data.  I can't reproduce your output, the program works
> fine on my system against gsl-1.6, the largest relative error is
> O(1e-14) i.e. close to machine precision.
> 
> It looks like you have a machine- or compiler- specific problem. I'd
> suggest running "make check" for the whole library to see if it throws
> up any other problems.

steed.c
Description: Text Data