Implement ufuncs for CSPHJY, SPHJ, SPHY, CSPHIK, SPHI, SPHIK (allow fast vectorization)

[certik]

Currently, the only way to use vectorized spherical Bessel functions is using a variation of:

In [3]: import numpy as np

In [4]: jn_vect = np.vectorize(sph_jn)


In [9]: jn_vect(0, [0.1, 0.2, 0.3, 0.5])
Out[9]:
(array([ 0.99833417,  0.99334665,  0.98506736,  0.95885108]),
 array([-0.03330001, -0.06640038, -0.09910289, -0.16253703]))

In [10]: jn_vect([0] * 4, [0.1, 0.2, 0.3, 0.5])
Out[10]:
(array([ 0.99833417,  0.99334665,  0.98506736,  0.95885108]),
 array([-0.03330001, -0.06640038, -0.09910289, -0.16253703]))

But this is slow, per the docstring of vectorize:

     The `vectorize` function is provided primarily for convenience, not for
    performance. The implementation is essentially a for loop.

As such, the solution is to implement proper ufuncs support for these special functions on the C level. Per @pv's suggestion, nowadays, this is fairly simple to do, take a look at:
generate_ufuncs.py
specfun_wrappers.h
specfun_wrappers.c

[cairijun]

I think it is a bit difficult because the original versions of these functions return variable-length arrays. It seems that generate_ufuncs.py cannot deal with this situations. It may also violate the definition of the ufunc in http://docs.scipy.org/doc/numpy/reference/ufuncs.html

That is, a ufunc is a “vectorized” wrapper for a function that takes a fixed number of scalar inputs and produces a fixed number of scalar outputs.

[pv]

The point is more to just throw away unnecessary results and not return
them to Python-space. Moreover, I believe these functions are not very
challenging to calculate, so that a reimplementation can also be written
from scratch.

[tpudlik]

I have started working on this and hope to make a PR within a couple weeks. Because the current algorithm for SPHJ suffers from stability issues (#2165, #1641), I'm preparing a comparison of the algorithms described in the literature to find one robust for a broad range of orders and argument values. For real values of the argument, there is thorough prior work (Jablonski 1995), but for complex values the evidence it is less extensive (the best I've found is Liang-Wu Cai 2011).

[ev-br]

Have been done in gh-5623, so closing.