incorrect bessel functions of large order #4366

Closed
stevengj opened this Issue Sep 25, 2013 · 9 comments

Projects

None yet

5 participants

@stevengj
Member

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
Member

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
klkuhlm commented Sep 26, 2013

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
Member

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
Member

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
Member

I believe this is fixed now?

@jiahao
Member
jiahao commented Dec 29, 2013

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
Member
jiahao commented Feb 22, 2014

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

@jiahao jiahao closed this Feb 22, 2014
@JeffBezanson
Member

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

julia> besselj(2^31, 0.1)
0.0012489586587999257
@JeffBezanson JeffBezanson reopened this Feb 22, 2014
@jiahao jiahao added a commit that referenced this issue Feb 22, 2014
@jiahao jiahao Fix Bessel functions for very large order (#4366)
Our version of openlibm uses a recurrence relation to build higher order
Bessel functions from j0 and j1. For orders greater than int32_t, the
recurrence loop never runs thanks to integer overflow in the loop index.
The fix here is to simply throw a DomainError for such enormous values
of the order of the Bessel function.
b2fea25
@jiahao
Member
jiahao commented Feb 22, 2014

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

@jiahao jiahao closed this Feb 22, 2014
@jiahao jiahao added a commit that referenced this issue Feb 27, 2014
@jiahao jiahao Better fix for #4366
- Instead of throwing DomainError for besselj(n, x) when
  n>typemax(Int32), call besselj(float64(n), x) instead.
- Fixes also wraparound for n<typemin(Int32).

The behavior for n::Int is now symmetric and complements the runtime
dispatch behavior for n::FloatingPoint.
4b8526e
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment