[Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again

[Koichi Takahashi]

At this range, x is still not large enough to use the asymptotic form?
In gsl_sf_bessel_jl_e,

  else if(x > 1000.0 && x > 100.0*l*l)
  {
    //asymptotic
  }
  else
  {
    //CF1
  }

so, for example, 100 * 50 * 50 = 250,000.
For l=50, the iteration starts to exceed 10,000 around x=9900.
If we want to stick to 10,000 max iteration, we have to switch to
the asymptotic version with something like x > 3*l*l.  Maybe this is
too small?   If so I'd consider increasing the max iteration.
Maybe we should do both.  Or there can be a better method.

koichi




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