incorrect bessel functions of large order

[stevengj]

e.g. besselj(n, 0.1) should underflow to zero for n > 110 (as checked via WolframAlpha), but our implementation underflows for n > 100.

Worse, for sufficiently large orders it "un-underflows". For example, besselj(2.0^70, 0.1) gives 2.1386280643e-314 instead of 0.0.

Worse, the behavior for large orders is not pure: it depends upon previous function calls! Consider the following two examples (both after restarting Julia):

julia> besselj(2.0^70, 0.1)
2.1386280643e-314

julia> besselj(2^70, 0.1)
0.9975015620660401

(where 2^70 gives 0, so the the latter is the correct result for besselj(0, 0.1)), but:

julia> besselj(2^70, 0.1)
0.9975015620660401

julia> besselj(2.0^70, 0.1)
0.9975015620660401

where somehow the earlier result for 2^70 has "corrupted" the later result for 2.0^70!!!

[stevengj]

By the way,

julia> ccall(:jn, Float64, (Int32, Float64), 110, 0.1)
4.84e-322

(the correct result, more or less), so the jn function in openlibm doesn't have this problem (although it only supports Int32 arguments). The jn function is also significantly faster than besselj, so we should definitely have besselj for (sufficiently small) integer order call jn via the wonders of multiple dispatch (and similarly for bessely).

[klkuhlm]

Possibly this issue is related to the version of the Amos library used by Julia? The netlib version is an older version, which was later updated (but the updated version has a different license).

I came across a similar issue a few years ago in the octave implementation of bessel functions, which I think is the same as Julia's.

http://comments.gmane.org/gmane.comp.gnu.octave.bugs/7982

[stevengj]

We should really use the SLATEC version (patched with bugfixes as needed) for licensing reasons, since only the SLATEC version comes with a clear statement that it is in the public domain. (In contrast, anything in ACM TOMS is by default unusable.)

[ViralBShah]

I had tried slatec at one point, but just dropping it in did not even make our tests pass. It is probably in the history of openlibm-extras and can be reinstated from the git history. I never dug deep into what happened.

[JeffBezanson]

I believe this is fixed now?

[jiahao]

The purity issue seems to be fixed.

However, besselj still seems to have problems.

julia> besselj(2^29,0.1) #Takes ~8s to run
0.0

julia> besselj(float(2^29),0.1)
0.0

julia> besselj(2^31, 0.1)
0.049937526036242

julia> besselj(float(2^31), 0.1)
6.9318596292179e-310

[jiahao]

The underflowing cases reported in this issue now throw AmosExceptions, so I think we can consider this issue resolved. (ref. #4967)

[JeffBezanson]

The integer order version is still a bit wonky. I get

julia> besselj(2^31, 0.1)
0.0012489586587999257

[jiahao]

Thanks for the reminder for that case. This should now be fixed as well.