Nasty bug in jn(3)

[Steve Kargl]

There's a nasty little bug lurking in jn(3).

#include <stdio.h>
#include <math.h>

int
main(void)
{
	double z;
	int i, n;

	z = 2.4048255576957729;
	for (n = 2; n < 10; n++)
		printf("%d %e\n", n, jn(n,z));

	return (0);
}

troutmask:kargl[446] cc -o z testjn.c -lm
troutmask:kargl[447] ./z
2 4.317548e-01
3 -inf
4 4.069027e-02
5 -inf
6 3.247664e-03
7 -inf
8 7.495602e-05
9 -inf

I can assure you that -inf is not a valid value for
an integer order Bessel function at z = 2.40482555...
A quick inspection of e_jn.c suggest a similar
problem maybe found at other zeros of j0(x).

-- 
Steve