[Bug-gsl] Re: gsl_sf_bessel_jl bug

[Jonathan Taylor]

> 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().

I've also seen this problem (thought I'd raised a bug but can't find  
any reference to it now). The problem seems to be that  
gsl_sf_bessel_J_CF1 checks for overflow, but not for underflow. It  
should be easy to make a local fix to your own gsl build by adding  
the following analogous code after the overflow check:

const double RECUR_SMALL = GSL_SQRT_DBL_MIN;
if(fabs(An) < RECUR_SMALL || fabs(Bn) < RECUR_SMALL) {
      An *= RECUR_BIG;
      Bn *= RECUR_BIG;
      Anm1 *= RECUR_BIG;
      Bnm1 *= RECUR_BIG;
      Anm2 *= RECUR_BIG;
      Bnm2 *= RECUR_BIG;
}