Repository: mpmath/mpmath

Issue:
Bug: issue with spherical bessel functions: z is not at the denominator.
Hello, I have noticed that the spherical bessel functions [spherical_jn](https://github.com/mpmath/mpmath/blob/97c45f21e37d53070c813138bab83d6120f47df7/mpmath/functions/bessel.py#L197)  and [spherical_yn](https://github.com/mpmath/mpmath/blob/97c45f21e37d53070c813138bab83d6120f47df7/mpmath/functions/bessel.py#L224) are not providing correct results.

Here a gist with a comparison against the analogous functions of scipy: https://gist.github.com/Ferrangelo/67c5065b6594b6ec54c1038e11a5d02a

The problem lies in the formulas defining the functions. In these two lines, z is not at the denominator as it should be: [L221_bessel.py](https://github.com/mpmath/mpmath/blob/97c45f21e37d53070c813138bab83d6120f47df7/mpmath/functions/bessel.py#L221), [L243_bessel.py](https://github.com/mpmath/mpmath/blob/97c45f21e37d53070c813138bab83d6120f47df7/mpmath/functions/bessel.py#L243)

A quick fix would be to modify them as follows:

`return ctx.besselj(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/(2*z))`

`return ctx.bessely(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/(2*z))`

But then the division by zero should be handled somehow.

The functions were still passing the tests because they are tested only for `z = 1`, and it doesn't matter whether z is a the numerator or denominator if it is 1.

Discussion:

Fix (patch):
diff --git a/mpmath/functions/bessel.py b/mpmath/functions/bessel.py
index 42b7f20..d8fe37c 100644
--- a/mpmath/functions/bessel.py
+++ b/mpmath/functions/bessel.py
@@ -218,7 +218,7 @@ def spherical_jn(ctx, n, z):
     mpf('0.84147098480789639')
 
     """
-    return ctx.besselj(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/2*z)
+    return ctx.besselj(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/(2*z))
 
 @defun
 def spherical_yn(ctx, n, z):
@@ -240,7 +240,7 @@ def spherical_yn(ctx, n, z):
     mpf('-0.54030230586813965')
 
     """
-    return ctx.bessely(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/2*z)
+    return ctx.bessely(n + ctx.one/2, z) * ctx.sqrt(ctx.pi/(2*z))
 
 @defun_wrapped
 def whitm(ctx,k,m,z,**kwargs):


Test changes:
diff --git a/mpmath/tests/test_functions2.py b/mpmath/tests/test_functions2.py
index 0392ecb..cd925c8 100644
--- a/mpmath/tests/test_functions2.py
+++ b/mpmath/tests/test_functions2.py
@@ -2454,3 +2454,7 @@ def test_issue_908():
 def test_issue_637():
     assert hankel1(1, 1 + 30j).ae(-7.25495e-15 - 1.17346e-14j)
     assert hankel2(1, 1 - 30j).ae(-7.25495e-15 + 1.17346e-14j)
+
+def test_issue_991():
+    assert spherical_jn(0, 1.3).ae(0.74119860416707)
+    assert spherical_yn(0, 1.3).ae(-0.20576832971122)