On Mon, Feb 9, 2009 at 4:11 AM, Pauli Virtanen <pav@iki.fi> wrote:
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.
Yes, I agree - I already asked about this a few weeks ago after some problems with other functions. I think cephes and specfun are not reliable - R does not use it, they have their own implementation of core math functions (sometimes inspired from cephes/specfun, but not that often).
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.
Yes, it would be a lot of work - I think we should focus on the tests before rewriting some functions. I would like to have a core scipy.special which is reliable: bessel, gamma/digamma/co, chebychev, basically most functions in R core would be a good start - and already quite heavy from a work POV.
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.
we could use their test-suite, maybe. cheers, David