Hi,
Jonathan's fix suppresses the error.
I hope a new ver with this kind of fix would come
up soon as it seems fairly serious.
Below is the change I made.
Thanks a lot!
Koichi
--- specfunc/bessel.c-bak 2007-10-17 21:29:44.000000000 -0700
+++ specfunc/bessel.c 2007-10-17 21:43:55.000000000 -0700
@@ -487,6 +487,7 @@
double * ratio, double * sgn)
{
const double RECUR_BIG = GSL_SQRT_DBL_MAX;
+ const double RECUR_SMALL = GSL_SQRT_DBL_MIN;
const int maxiter = 10000;
int n = 1;
double Anm2 = 1.0;
@@ -504,6 +505,7 @@
while(n < maxiter) {
double old_fn;
double del;
+ double fabsAn, fabsBn;
n++;
Anm2 = Anm1;
Bnm2 = Bnm1;
@@ -513,7 +515,10 @@
An = Anm1 + an*Anm2;
Bn = Bnm1 + an*Bnm2;
- if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
+ fabsAn = fabs(An);
+ fabsBn = fabs(Bn);
+
+ if(fabsAn > RECUR_BIG || fabsBn > RECUR_BIG) {
An /= RECUR_BIG;
Bn /= RECUR_BIG;
Anm1 /= RECUR_BIG;
@@ -521,6 +526,14 @@
Anm2 /= RECUR_BIG;
Bnm2 /= RECUR_BIG;
}
+ else if(fabsAn < RECUR_SMALL || fabsBn < RECUR_SMALL) {
+ An *= RECUR_BIG;
+ Bn *= RECUR_BIG;
+ Anm1 *= RECUR_BIG;
+ Bnm1 *= RECUR_BIG;
+ Anm2 *= RECUR_BIG;
+ Bnm2 *= RECUR_BIG;
+ }
old_fn = fn;
fn = An/Bn;
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
