Fix jn() precision problems around zero points of j0()

Fix jn() precision problems around zero points of j0()

There appears to be a really nasty bug in jn() from fdlibm, which is
the foundation for most libm implementations (including glibc libm).
The zeroth-order j0() and first-order j1() cylindrical Bessel functions
are used to recursively generate the jn() value, but only the zeroth-order
Bessel function is used to normalize it; however, each of the functions
gets highly imprecise (approaching "bogus") near its zero point, making
the jn() value itself bogus.

But in fact, the zero points of j0() and j1() never coincide, thus j1()
should be used in case it is more precise than j0(). (That is, simply
when its value is further from zero.)

As an example, 2.4048255576957729_8 is the first zero of j0().
The proper value as calculated by Mathematica is 0.19899990535769...
However, jn() returns -inf on 64-bit arch, or 0.185007 on 32-bit arch.
With the proposed patch below, the returned value is 0.199000.

The fix is based on work by Steve Kargl and Tobias Burnus.