bug in special function iv for large argumen

[mbakker7]

The modified Bessel function iv returns incorrect values for large argument.

iv(0,100) gives 5.72185663838e+041

while the equivalent is

jv(0,complex(0,100)) (1.07375170713e+042+0j )

This latter results is the correct answer, as verified in Abramowitz and Stegun, Table 9.11, Page 428. Using the jv is a work-around, but this should probably be fixed. The two results start to deviate for arguments over about 50, with jv giving the correct answer. The same goes wrong for iv of other integer order (I did not check non-integer orders).

On a related note, there are implementations for Bessel functions of integer order (jn, kn) but not for the modified Bessel function In. I guess this is because the function would be called 'in', but it would be nice to have a special function for integer order, and I am pretty sure they are around.

Thanks, Mark

[BlackShift]

The problem is not the large argument, but the fact that the argument is real:

>>> scipy.special.iv(0,100+0j).real
1.0737517071310738e+42
>>> scipy.special.iv(0,100)
5.7218566383804826e+41

This bug has been there for some time (at least a year) and I think it has been mentioned before, however I cannot find a ticket. The work around is to always include an imaginary part, even when you are only interested in real numbers.

[pv]

Cephes's implementation of hyp1f1 was buggy (didn't properly switch to the asymptotic expansion), and this implementation was used to implement Iv.

Fixed in r5477 and r5478.

[pv]

Not completely fixed yet: simultaneous large v and large x yield NaNs?:

>>> import scipy.special as sc
>>> sc.iv(1200, 1300)
nan

Moreover, for v = -1/2, -3/2, ..., Cephes tries to evaluate hyperg at a singularity:

>>> sc.iv(-0.5, 1)
0.0
>>> sc.iv(-0.5 + 1e-13, 1)
1.2312002145930347

[pv]

Fixed in r5641, r5640, r5639.

Reviewed at http://codereview.appspot.com/20078