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 IEEE754compliant 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...

From http://lists.gnu.org/archive/html/buggsl/201311/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
