[SciPy-dev] Bessel functions from Boost

[Pauli Virtanen]

Sun, 08 Feb 2009 11:39:55 -0700, Charles R Harris wrote:
[clip]
> I think it's a good idea. It would also be nice if we picked best of
> breed from several libraries to make up our own special functions
> collection. For instance, there are several log gamma functions. 

I think we need to do this eventually, even if it means lots of work. At 
points the Cephes and Specfun codes seem like the author has not wanted 
to bother with the best possible algorithm, which leads to problems in 
corner cases.

> That's
> probably a big job though and we would need extensive tests for the
> functions before trying it. Do you know of any other project that has
> put together such a test suite?

For Bessel functions we can easily test against the AMOS library, which 
appears to be reliable --- unless the order is close to negative integers 
in which case there can be cancellation errors of the order of 1e-6 in 
the reflection formulas.

Boost itself has tests for its special functions, these are spot tests at 
precomputed points, from 1e-2..1e2 magnitudes in both parameters. [As an 
aside, I noticed Boost's I(v,x) overflows to infty somewhat earlier than 
necessary for very large orders, though.]

GSL has similar point tests. (But it's GPLed.)

Netlib/Specfun has a test suite: http://netlib.org/specfun/ ; in F77.

Anyway, point tests across some magnitudes of parameters should be easy 
to generate. What takes more work is checking the behavior of the 
functions in transition regions where the method of computation changes, 
and asymptotic behavior (overflows, etc.) at large or small parameters 
and near singularities.

> Does boost have versions for log1p? We need a better implementation in
> numpy itself.

It has. 

http://svn.boost.org/svn/boost/trunk/boost/math/special_functions/
log1p.hpp

It's cluttered by C++ templates, but the algorithm looks like some 
serious effort has been put into it.

-- 
Pauli Virtanen