[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [bug #40755] Sporadic nan's from gsl_sf_bessel_Jn an related f
From: |
Peter Barfuss |
Subject: |
[Bug-gsl] [bug #40755] Sporadic nan's from gsl_sf_bessel_Jn an related functions |
Date: |
Sun, 14 Oct 2018 08:27:46 -0400 (EDT) |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/68.0.3440.75 Safari/537.36 |
Follow-up Comment #1, bug #40755 (project gsl):
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...
_______________________________________________________
Reply to this item at:
<https://savannah.gnu.org/bugs/?40755>
_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [bug #40755] Sporadic nan's from gsl_sf_bessel_Jn an related functions,
Peter Barfuss <=