[bug-gsl] Re: gsl_sf_bessel_jl bug (2007-10-25)
|
From: |
Brian Gough <bjg@gnu.org> |
|
Subject: |
Re: gsl_sf_bessel_jl bug |
|
Date: |
Thu, 25 Oct 2007 12:34:54 +0100 |
|
To: |
Koichi Takahashi <ktakahashi@..........> |
|
Cc: |
bug-gsl@gnu.org, Jonathan Taylor <j.m.taylor@............> |
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
