GNU 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: |
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:
|
Rhys Ulerich <rhysu> |
No files currently attached
Depends on the following items: None found
Items that depend on this one: None found
Follows 1 latest change.
Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|
2013-11-30 | rhysu | Carbon-Copy | - | Added -email is unavailable- |
Powered by Savane 3.13-f8d8.
Corresponding source code
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...