scipy/scipy — pull request #509 by pv (2013-04-13 18:58:01+00:00)
Title: BUG: special: use AMOS for computing all Bessel functions

The Bessel function implementations from Cephes are somewhat bad at dealing correctly with large order functions; yv and kn are missing an uniform large-order expansions, and the jv code seems partially buggy in this regime. These issues are not straightforward to fix; while the correct algorithms to do the evaluation are well-known, there's no point in writing new code for evaluating Bessel functions in 2013.

The AMOS implementation, which is already in Scipy, is written considerably more carefully and deals with these things correctly. It however is complex-valued implementation. It also has some limitations and does not compute results for |z| ~> 1e8.

Another probably good implementation is in the Boost library. Boost is unfortunately too large to bundle with Scipy (~ 20 MB of code) just for a couple of functions. Porting the YV/JV implementation to C is in principle doable, but I'm not going to spend time on that.

What this PR then does is:

- Evaluate kn, jv, yv using AMOS even for real-valued parameters.
- If AMOS refuses to compute, use the Cephes codes (for jv, yv they do contain apparently working large-z expansions), or known asymptotic behavior (for kv/kn).
- Fix some corner case bugs, such as signs of infinities. The sign infinity for jv at negative non-integer orders is still not correct, but probably doesn't matter in practice.

The drawback with using AMOS is that it is somewhat slower. However, I don't expect that this will matter in practice.

Also, the Cephes Bessel I code was thoroughly fixed some time ago, and doesn't seem to need any further fixes.

Closes [Trac #1740](http://projects.scipy.org/scipy/ticket/1740), [Trac #851](http://projects.scipy.org/scipy/ticket/851) and [Trac #1589](http://projects.scipy.org/scipy/ticket/1589)