Bessel functions in Python that work with large exponents

Question

I've got some code that uses the modified Bessel functions of both 1st and 2nd order (iv and kv). Annoyingly they seem to have limits, those are iv(0,713) and kv(0,697), add one to each and you get infinity and 0 respectively. This is a problem for me because I need to use values higher than this, often up to 2000 or more. When I try to divide by these I end up dividing by 0 or by infinity which means I either get errors or zeros, neither of which I want.

I'm using the scipy bessel functions, are there any better functions that can cope with much smaller and much larger numbers, or a way of modifying Python to work with these big numbers. I'm unsure what the real issue here is as to why Python can't work these out much beyond 700, is it the function or is it Python?

I don't know if Python is already doing it but I'd only need the first 5-10 digits *10^x for example; that is to say I wouldn't need all 1000 digits, perhaps this is the problem with how Python is working it out compared to how Wolfram Alpha is working it out?

Answer 1

The iv and kv functions in Scipy are more or less as good as you can get if using double precision machine floating point. As noted in the comments above, you are working in the range where the results overflow from the floating point range.

You can use the mpmath library, which does adjustable precision (software) floating point, to get around this. (It's similar to MPFR, but in Python):

In [1]: import mpmath

In [2]: mpmath.besseli(0, 1714)
mpf('2.3156788070459683e+742')

In [3]: mpmath.besselk(0, 1714)
mpf('1.2597398974570405e-746')