[Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again

[Jonny Taylor]

While the symptom is the same, the cause is different. For those  
numbers it seems that 10,000 iterations is simply not enough.  
Interestingly, In fact, in all three cases it requires less than 50  
additional iterations to converge!?

The naive fix is simply to increase the maximum permitted number of  
iterations, but a more sophisticated fix would probably need to  
justify a maximum number of iterations, and propose an alternative  
method of generating the result in the cases where the number of  
iterations is prohibitive...

Jonny


On 2 Dec 2007, at 08:11, Koichi Takahashi wrote:

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