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.
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
