On Mon, Feb 9, 2009 at 3:26 AM, Pauli Virtanen <pav@iki.fi> wrote:
Mon, 09 Feb 2009 01:20:11 -0800, Ondrej Certik wrote: [clip]
Right. For double precision I think mpmath is not so fast. Fredrik, is it difficult to make mpmath fast even for double precision?
Last time I asked:
http://groups.google.com/group/mpmath/browse_thread/thread/ bca53c3382945c34/<http://groups.google.com/group/mpmath/browse_thread/thread/%0Abca53c3382945c...>
you replied:
" SciPy already provides a truckload of machine precision special functions, with excellent (fast and robust) implementations. It'd be hard to top that. "
But apparently, maybe mpmath can be useful.
I'd say that mpmath faces the same robustness and testing issues as Scipy with regard to special functions. (In addition, since it's written in Python, I'd assume it also faces additional performance issues.)
Also, algorithms that work well in arbitrary precision might not work for limited precision, due to loss of precision or under/overflows in intermediate steps. Looking at the Bessel function implementations in mpmath/functions.py, I'd say that at least besselj and besseli would face overflow issues for large arguments if they were working in double precision. This kind of issues are actually the most difficult to get right.
To clarify: Definitely I think that mpmath is great work, and I'm happy to see people working on it, including improvements to its special function library.
But at the present, I think the path of least resistance for Scipy is to continue using, testing, and improving existing implementations of special function codes, written in C or F77, directly for limited precision.
Here are some references for both testing and implementations: http://math.nist.gov/mcsd/Reports/2001/nesf/ . I found Cody's book(s) a good reference back in the day. I wonder if anything happened with this proposal: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.7298 ? Note the references to Cody which are to works with both implementations and tests. I think Cody was responsible for specfun, have there been problems with that package? Chuck