Repository: mpmath/mpmath (python)

Issue:
Spherical Bessel functions
Spherical Bessel functions of the first kind are defined as,

```
from mpmath import mpf, sqrt, pi, besselj

def sph_jn(n, z):
    return besselj(n + mpf(1)/2, z)*sqrt(pi/(2*z))
```

For positive integer `n`, these functions have closed-form expressions in terms of trigonometric functions.  For example, `sph_jn(0, x) = sin(x)/x` (see http://dlmf.nist.gov/10.49.E3).

For negative values of `x`, the `mpmath` implementation of these functions has the wrong sign.  For example,

```
>>> sph_jn(0, -1)
mpc(real='-0.84147098480789639', imag='0.0')
```

The magnitude is correct, but the sign is not!

EDIT: The same problem plagues spherical Bessel functions of the second kind, defined analogously in terms of `bessely`.


Pull request:
Support spherical Bessel functions (jn/yn)
Naming scheme follows to SciPy.  Closes #319

Patch:
diff --git a/docs/functions/bessel.rst b/docs/functions/bessel.rst
index 04f9bdb..6510c9a 100644
--- a/docs/functions/bessel.rst
+++ b/docs/functions/bessel.rst
@@ -33,6 +33,13 @@ Hankel functions
 .. autofunction:: mpmath.hankel2
 
 
+Spherical Bessel functions
+..........................
+
+.. autofunction:: mpmath.spherical_jn
+.. autofunction:: mpmath.spherical_yn
+
+
 Kelvin functions
 ................
 
diff --git a/mpmath/__init__.py b/mpmath/__init__.py
index c036c13..b546432 100644
--- a/mpmath/__init__.py
+++ b/mpmath/__init__.py
@@ -392,6 +392,8 @@
 besselk = mp.besselk
 besseljzero = mp.besseljzero
 besselyzero = mp.besselyzero
+spherical_jn = mp.spherical_jn
+spherical_yn = mp.spherical_yn
 hankel1 = mp.hankel1
 hankel2 = mp.hankel2
 struveh = mp.struveh
diff --git a/mpmath/functions/bessel.py b/mpmath/functions/bessel.py
index 7c69afa..371987e 100644
--- a/mpmath/functions/bessel.py
+++ b/mpmath/functions/bessel.py
@@ -192,6 +192,55 @@ def terms():
                 -ctx.j*ctx.bessely(n,x,**kwargs)]
     return ctx.sum_accurately(terms)
 
+@defun
+def spherical_jn(ctx, n, z):
+    r"""
+    Spherical Bessel function of the first kind.
+
+    This function is a solution to the spherical Bessel equation
+
+    .. math ::
+        z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2}
+          + 2z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 - \nu(\nu + 1)) w = 0.
+
+    It can be defined as
+
+    .. math ::
+        j_\nu(z) = \sqrt{\frac{\pi}{2z}} J_{\nu + \frac{1}{2}}(z),
+
+    where `J_\nu(z)` is the Bessel function of the first kind.
+
+    **Examples**
+
+    >>> from mpmath import spherical_jn
+    >>> spherical_jn(0, 1)
+    mpf('0.84147098480789639')
+
+    """
+    return ctx.besselj(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/2*z)
+
+@defun
+def spherical_yn(ctx, n, z):
+    r"""
+    Spherical Bessel function of the second kind.
+
+    This function is another solution to the spherical Bessel equation, and
+    linearly independent from `j_n`. It can be defined as
+
+    .. math ::
+        j_\nu(z) = \sqrt{\frac{\pi}{2z}} Y_{\nu + \frac{1}{2}}(z),
+
+    where `Y_\nu(z)` is the Bessel function of the second kind.
+
+    **Examples**
+
+    >>> from mpmath import spherical_yn
+    >>> spherical_yn(0, 1)
+    mpf('-0.54030230586813965')
+
+    """
+    return ctx.bessely(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/2*z)
+
 @defun_wrapped
 def whitm(ctx,k,m,z,**kwargs):
     if z == 0:


Test patch:
diff --git a/mpmath/tests/test_functions2.py b/mpmath/tests/test_functions2.py
index 610f7ca..05edd39 100644
--- a/mpmath/tests/test_functions2.py
+++ b/mpmath/tests/test_functions2.py
@@ -11,8 +11,8 @@
                     kei, ker, laguerre, lambertw, ldexp, legendre, legenp,
                     legenq, lerchphi, li, log, lower_gamma, meijerg, mp, mpc,
                     mpf, nan, ncdf, npdf, nthroot, pi, qp, quadts, shi, si,
-                    spherharm, sqrt, struveh, struvel, upper_gamma, whitm,
-                    whitw, zeta)
+                    spherharm, spherical_jn, spherical_yn, sqrt, struveh,
+                    struvel, upper_gamma, whitm, whitw, zeta)
 from mpmath.libmp import BACKEND, NoConvergence
 
 
@@ -86,6 +86,8 @@ def test_bessel():
     assert besselk(0,j).ae(-0.13863371520405399968-1.20196971531720649914j)
     assert (besselk(3, 10**10) * mpf(10)**4342944824).ae(1.1628981033356187851)
     assert besselk(1,inf) == 0
+    assert spherical_jn(0, 1).ae(0.841470984807896)
+    assert spherical_yn(0, 1).ae(-0.54030230586814)
     # test for issue 331, bug reported by Michael Hartmann
     for n in range(10,100,10):
         mp.dps = n