[SciPy-Dev] Time efficient code to calculate range of Bessel Functions

[Vlad Blazhko]

Hello,

I would like to contribute time efficient and numerically stable cython
code for calculating a range of Bessel's functions of first and second
order (jv and yv) to the subpackage special in scipy.

It is useful since often you have a sum that uses jv of the one argument,
but of several orders that increase by one. They can be efficiently
computed through the recurrence relation. While recurrence relation works
fine for the yv function, it quickly diverges for the jv function, so for
the jv function I have implemented Miller's recurrence algorithm, which is
stable.

In my experience it is usually faster 5-6 times then using just jv and yv,
and such code usually lies in the heart of simulations (in hot places).

So I propose the following:

Functions names: jv_range, yv_range
Example of code:
1) jv_range(v_from=0.3, n=5, z=np.ndarray) -> computes values for orders
0.3, 1.3, 2.3, 3.3, 4.3 and puts them in last axis, i.e. the shape of
result is z.shape + (n,)
2) jv_range(v_from=-2.3, n=6, z=np.ndarray) -> computes values for orders
-2.3, -1.3, -0.3, 0.7, 1.7, 2.7

If the functions don't have enough accuracy for given orders, then they
will fallback to jv and yv functions.

Do you think such functions will be a good contribution to the special
subpackage?
Feel free to give a suggestion on the functions interface
And please let me know if you would be so kind to review my PR.

Best regards,
Vlad Blazhko
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/scipy-dev/attachments/20210116/17474435/attachment.html>