ENH: special: add the analytic continuation of Riemann zeta

scipy/special/_cephes.pxd

@@ -73,6 +73,7 @@ cdef extern from "cephes.h" nogil:
     double poch(double x, double m)
     double psi(double x)
     double rgamma(double x)
+    double riemann_zeta(double x)
     double round(double x)
     int shichi(double x, double *si, double *ci)
     int sici(double x, double *si, double *ci)

scipy/special/add_newdocs.py

@@ -7140,6 +7140,11 @@ def add_newdoc(name, doc):
 
     """)
 
+add_newdoc("_riemann_zeta",
+    """
+    Internal function, use `zeta` instead.
+    """)
+
 add_newdoc("_struve_asymp_large_z",
     """
     _struve_asymp_large_z(v, z, is_h)

scipy/special/basic.py

@@ -13,8 +13,8 @@
                    extract, inexact, nan, zeros, sinc)
 from . import _ufuncs as ufuncs
 from ._ufuncs import (mathieu_a, mathieu_b, iv, jv, gamma,
-                      psi, _zeta, hankel1, hankel2, yv, kv, ndtri,
-                      poch, binom, hyp0f1, wofz, _voigt)
+                      psi, hankel1, hankel2, yv, kv, ndtri,
+                      poch, binom, hyp0f1)
 from . import specfun
 from . import orthogonal
 from ._comb import _comb_int
@@ -2246,7 +2246,7 @@ def voigt(x, sigma=1.0, gamma=1.0, mu=0.0):
     ----------
     .. [1] https://en.wikipedia.org/wiki/Voigt_profile
     """
-    return _voigt(x, sigma, gamma, mu)
+    return ufuncs._voigt(x, sigma, gamma, mu)
 
 
 def zeta(x, q=None, out=None):
@@ -2302,5 +2302,6 @@ def zeta(x, q=None, out=None):
 
     """
     if q is None:
-        q = 1
-    return _zeta(x, q, out)
+        return ufuncs._riemann_zeta(x, out)
+    else:
+        return ufuncs._zeta(x, q, out)

scipy/special/cephes/cephes_names.h

@@ -71,6 +71,7 @@
 #define poch cephes_poch
 #define psi cephes_psi
 #define rgamma cephes_rgamma
+#define riemann_zeta cephes_riemann_zeta
 #define round cephes_round
 #define shichi cephes_shichi
 #define sici cephes_sici

scipy/special/cephes/zetac.c

@@ -29,10 +29,6 @@
  * Extension of the function definition for x < 1 is implemented.
  * Zero is returned for x > log2(NPY_INFINITY).
  *
- * An overflow error may occur for large negative x, due to the
- * Gamma function in the reflection formula, so nan is returned
- * for x < -30.8148.
- *
  * ACCURACY:
  *
  * Tabulated values have full machine accuracy.
@@ -181,7 +177,7 @@ static double TAYLOR0[10] = {
 
 extern double MACHEP;
 
-static double zetac_reflection(double);
+static double zeta_reflection(double);
 static double zetac_smallneg(double);
 static double zetac_positive(double);
 
@@ -201,14 +197,37 @@ double zetac(double x)
 	return zetac_smallneg(x);
     }
     else if (x < 0.0) {
-	return zetac_reflection(-x);
+	return zeta_reflection(-x) - 1;
     }
     else {
 	return zetac_positive(x);
     }
 }
 
 
+/*
+ * Riemann zeta function
+ */
+double riemann_zeta(double x)
+{
+  if (npy_isnan(x)) {
+    return x;
+  }
+  else if (x == -NPY_INFINITY) {
+    return NPY_NAN;
+  }
+  else if (x < 0.0 && x > -0.01) {
+    return 1 + zetac_smallneg(x);
+  }
+  else if (x < 0.0) {
+    return zeta_reflection(-x);
+  }
+  else {
+    return 1 + zetac_positive(x);
+  }
+}
+
+
 /*
  * Compute zetac for positive arguments
  */
@@ -289,14 +308,14 @@ static NPY_INLINE double zetac_smallneg(double x)
  * Compute zetac using the reflection formula (see DLMF 25.4.2) plus
  * the Lanczos approximation for Gamma to avoid overflow.
  */
-static NPY_INLINE double zetac_reflection(double x)
+static NPY_INLINE double zeta_reflection(double x)
 {
     double s, hx, x_shift;
 
     hx = x / 2;
     if (hx == floor(hx)) {
 	/* Hit a zero of the sine factor */
-	return -1.0;
+	return 0;
     }
 
     /* Group large terms together to prevent overflow */
@@ -305,5 +324,5 @@ static NPY_INLINE double zetac_reflection(double x)
     x_shift = fmod(x, 4);
     s *= -SQRT_2_PI * sin(0.5 * NPY_PI * x_shift);
     s *= lanczos_sum_expg_scaled(x + 1) * zeta(x + 1, 1);
-    return s - 1.0;
+    return s;
 }

scipy/special/functions.json

@@ -1387,5 +1387,10 @@
         "cephes.h": {
             "zetac": "d->d"
         }
+    },
+    "_riemann_zeta": {
+	"cephes.h": {
+	    "riemann_zeta": "d->d"
+	}
     }
 }

scipy/special/tests/test_mpmath.py

@@ -1849,12 +1849,21 @@ def test_wrightomega(self):
                             lambda z: _mpmath_wrightomega(z, 25),
                             [ComplexArg()], rtol=1e-14, nan_ok=False)
 
-    def test_zeta(self):
+    def test_hurwitz_zeta(self):
         assert_mpmath_equal(sc.zeta,
                             exception_to_nan(mpmath.zeta),
                             [Arg(a=1, b=1e10, inclusive_a=False),
                              Arg(a=0, inclusive_a=False)])
 
+    def test_riemann_zeta(self):
+        assert_mpmath_equal(
+            sc.zeta,
+            mpmath.zeta,
+            [Arg(-100, 100)],
+            nan_ok=False,
+            rtol=1e-13,
+        )
+
     def test_zetac(self):
         assert_mpmath_equal(sc.zetac,
                             lambda x: mpmath.zeta(x) - 1,

scipy/special/tests/test_zeta.py

@@ -9,31 +9,37 @@ def test_zeta():
     assert_allclose(sc.zeta(2,2), np.pi**2/6 - 1, rtol=1e-12)
 
 
-def test_zeta_1arg():
-    assert_allclose(sc.zeta(2), np.pi**2/6, rtol=1e-12)
-    assert_allclose(sc.zeta(4), np.pi**4/90, rtol=1e-12)
+def test_zetac():
+    # Expected values in the following were computed using Wolfram
+    # Alpha's `Zeta[x] - 1`
+    x = [-2.1, 0.8, 0.9999, 9, 50, 75]
+    desired = [
+        -0.9972705002153750,
+        -5.437538415895550,
+        -10000.42279161673,
+        0.002008392826082214,
+        8.881784210930816e-16,
+        2.646977960169853e-23,
+    ]
+    assert_allclose(sc.zetac(x), desired, rtol=1e-12)
 
 
-def test_zetac():
-    assert_equal(sc.zetac(0), -1.5)
-    assert_equal(sc.zetac(1.0), np.inf)
-    # Expected values in the following were computed using
-    # Wolfram Alpha `Zeta[x] - 1`:
-    rtol = 1e-12
-    assert_allclose(sc.zetac(-2.1), -0.9972705002153750, rtol=rtol)
-    assert_allclose(sc.zetac(0.8), -5.437538415895550, rtol=rtol)
-    assert_allclose(sc.zetac(0.9999), -10000.42279161673, rtol=rtol)
-    assert_allclose(sc.zetac(9), 0.002008392826082214, rtol=rtol)
-    assert_allclose(sc.zetac(50), 8.881784210930816e-16, rtol=rtol)
-    assert_allclose(sc.zetac(75), 2.646977960169853e-23, rtol=rtol)
-
-
-def test_zetac_negative_even():
-    pts = [-2, -50, -100]
-    for p in pts:
-        assert_equal(sc.zetac(p), -1)
-
-
-def test_zetac_inf():
-    assert_equal(sc.zetac(np.inf), 0.0)
-    assert_(np.isnan(sc.zetac(-np.inf)))
+def test_zetac_special_cases():
+    assert sc.zetac(np.inf) == 0
+    assert np.isnan(sc.zetac(-np.inf))
+    assert sc.zetac(0) == -1.5
+    assert sc.zetac(1.0) == np.inf
+
+    assert_equal(sc.zetac([-2, -50, -100]), -1)
+
+
+def test_riemann_zeta_special_cases():
+    assert np.isnan(sc.zeta(np.nan))
+    assert sc.zeta(np.inf) == 1
+    assert sc.zeta(0) == -0.5
+
+    # Riemann zeta is zero add negative even integers.
+    assert_equal(sc.zeta([-2, -4, -6, -8, -10]), 0)
+
+    assert_allclose(sc.zeta(2), np.pi**2/6, rtol=1e-12)
+    assert_allclose(sc.zeta(4), np.pi**4/90, rtol=1e-12)