[SciPy-dev] Bessel functions from Boost
Pauli Virtanen
pav@iki...
Sun Feb 8 13:11:24 CST 2009
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
More information about the Scipy-dev
mailing list