[Bug-gsl] gsl_sf_bessel__J_CF1 bug again

[Koichi Takahashi]

Hi,

Unfortunately, I found some additional sets of numbers that still
crash gsl_sf_bessel_J_CF1() even with the cvs version.   Symptom
is exactly the same as what I reported before.  I tested on
x86_64.


main()
{
//    this used to crash, but now fixed in the current cvs.
//    double a = gsl_sf_bessel_jl( 30,  3875.6138424501978 );

//    at least the following three sets of values still crashes.
    double a = gsl_sf_bessel_jl( 49, 9912.6308956132361 );
//    double a = gsl_sf_bessel_jl( 49, 9950.3478935215589 );
//    double a = gsl_sf_bessel_jl( 52, 9930.5181281798232 );


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


Let me know if there is anything I could do to help you fixing
this issue.

thanks,
Koichi




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