Re: [Bug-gsl] gsl_sf_bessel_jl bug

[Brian Gough]

At Tue, 16 Oct 2007 17:58:33 +0100,
Jonathan Taylor 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().
> 
> 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;
> }

Thanks, I've applied this patch.

-- 
Brian Gough