[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again

From: Koichi Takahashi
Subject: [Bug-gsl] Re: gsl_sf_bessel__J_CF1 bug again
Date: Sun, 02 Dec 2007 06:06:11 -0800
User-agent: Thunderbird (X11/20071022)

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)

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.


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


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


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

//    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 );


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


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

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:

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.

reply via email to

[Prev in Thread] Current Thread [Next in Thread]