bugGNU Scientific Library - Bugs: bug #40755, Sporadic nan's from...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #40755: Sporadic nan's from gsl_sf_bessel_Jn an related functions

Submitter:  Rhys Ulerich <rhysu>
Submitted:  Sat 30 Nov 2013 04:27:29 PM UTC
   
 
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  None
Assigned to:  None Open/Closed:  Open
Release: 

Sun 14 Oct 2018 12:27:45 PM UTC, comment #1: 

This is a silly bug. The issue is as follows: the Hankel asymptotic is fundamentally numerically unstable when invoked with (n*n+1) > (2*x). You will see that the codepath taken when you start getting NaNs is via gsl_sf_bessel_Jnu_asympx_e & if you print the values of t, P & Q in the loop you will see the iterates wildly diverge & you get either P being -Inf or Q being +Inf after ~40 to ~50 iterations. Shortly thereafter you get a NaN somewhere & it proceeds to gremlin its way to the output.

Now, why is this codepath being taken? Well, it specifically is coming from the following branch in specfunc/bessel_Jn.c:
} else if(GSL_ROOT4_DBL_EPSILON * x > (n*n+1.0)) {

Now, GSL_ROOT4_DBL_EPSILON is 2^-13 on an IEEE754-compliant system (i.e. every sane one), and 51471.451913 2^-13 = ~6.28314 which shouldn't be larger than 1+46340^2 or 1+46341^2. And indeed, it is not for the former value. Which is less than 2^31. But it *is larger than the second value, because the n*n+1 is computed as a signed int32_t, meaning the code compares 6.28314 against -4633 & concludes that the former is indeed greater than the latter.

Solutions: cast n to either a uint64_t or a double in the comparison, replace the check with if (x > ((n/GSL_ROOT2_DBL_EPSILON)*(n/GSL_ROOT2_DBL_EPSILON))+1.0) or some similar workaround for the overflow, or... just take out this branch entirely? Every case it hits can just as well be handled via CF1 + backward recurrence... in fact, this is exactly what fdlibm jn(n, x) does here, & it works nicely...

Peter Barfuss <bofh453>
Sat 30 Nov 2013 04:27:29 PM UTC, original submission:  

From http://lists.gnu.org/archive/html/bug-gsl/2013-11/msg00008.html and behavior confirmed on GSL 1.16:

I think this is distinct from other bugs in the gsl bug database. No error is reported by gsl_sf_bessel_Jn_e.

$ cat testgsl.c
/*
** Bug in gsl_sf_bessel_Jn (and related functions)
*/

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf_bessel.h>

int main (int argc, char **argv) {
  double x = M_PI 128 127.999;
  int nmin = 46335, nmax = 46345, i;
  for (i = nmin; i <= nmax; ++i) {
    double gsl_res = gsl_sf_bessel_Jn (i, x);
    double sys_res = jn (i, x);
    printf ("For (%d, %f), gsl gave %f, jn() gave %f\n", i, x, gsl_res,
        sys_res);
  }
  return 0;
}
$ gcc -O testgsl.c -o testgsl -lgsl -lgslcblas -lm
$ ./testgsl
For (46335, 51471.451913), gsl gave -0.004580, jn() gave -0.004580
For (46336, 51471.451913), gsl gave -0.005310, jn() gave -0.005310
For (46337, 51471.451913), gsl gave -0.004980, jn() gave -0.004980
For (46338, 51471.451913), gsl gave -0.003657, jn() gave -0.003657
For (46339, 51471.451913), gsl gave -0.001604, jn() gave -0.001604
For (46340, 51471.451913), gsl gave 0.000769, jn() gave 0.000769
For (46341, 51471.451913), gsl gave -nan, jn() gave 0.002988
For (46342, 51471.451913), gsl gave -nan, jn() gave 0.004612
For (46343, 51471.451913), gsl gave -nan, jn() gave 0.005316
For (46344, 51471.451913), gsl gave -nan, jn() gave 0.004961
For (46345, 51471.451913), gsl gave -nan, jn() gave 0.003618


Thanks
Paul

Rhys Ulerich <rhysu>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by bofh453 (Posted a comment)
  • -email is unavailable- added by rhysu (Submitted the item)
  • -email is unavailable- added by rhysu
  •  

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-11-30 rhysu Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code