[SciPy-user] special.jn(3,x) became inexact in recent versions
aronovitch at gmail.com
Sat Jan 20 18:36:48 CST 2007
I've been holding back from installing recent release of scipy for a
long time now - the reason being some numerical errors I did not have
time to track down.
I just found out the problem (by tracing the calculations down to a
binary lib function...), so here's a description - in case it would help
some other users:
in the old version:
special.jn(3,62) gives 0.101359688424111 (correct)
in new version:
special.jn(3,62) gives -0.0027352572609119909 (wrong!)
(plotting a graph of jn in this region shows several "jumps")
quoting from the C source (jn.c taken from netlib's source)
* Not suitable for large n or x. Use jv() instead.
indeed special.jv(3.62) gives the correct result (0.101)
Older releases of cephes used the original names (j0,jn,...), while
the newer release contains macros converting these names to "cephes_jn" etc.
It is possible that the older cephes module did not really use the
cephes functions, but the ones from glibc instead (glibc contains many
mathematical functions, including a "jn", which does work for large x).
I wonder how do cephes functions compare to glibc ones
performance-wise. Might be a good idea to generate ufunc wrappers to
those (on platforms where glibc is available).
More information about the SciPy-user