Re: [Bug-gsl] gsl_sf_bessel_jl bug

[Koichi Takahashi]

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;