[SciPy-user] special.jn(3,x) became inexact in recent versions

Amit Aronovitch 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")

 The solution:
 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)

Possible Explanation:
  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).

       Amit A.

More information about the SciPy-user mailing list