Repository: scipy/scipy
Pull request: ENH: special: Extend Riemann Zeta function to complex inputs
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
https://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes #14073

#### What does this implement/fix?
<!--Please explain your changes.-->
This PR extends the Riemann Zeta function to complex inputs. It follows a Python reference implementation written by @mdhaber and posted in a Colab notebook [here](https://colab.research.google.com/drive/1zMDSAJlXCLRqMMtJ0e9nDGjQ8iZCnmn5?usp=sharing). There are a couple of minor tweaks I'm aware of which could improve accuracy slightly, but I will leave those for a follow-up to make this easier to review. As can be seen in the linked notebook, results are sufficiently accurate for this to get in as is.

The primary ways the implementation here differs from the reference implementation are that

1. Coefficients for expansions are stored in lookup tables.
2. Some factors (well logs of factors) in the Euler-Maclaurin sum are updated incrementally in each iteration, rather than recomputed from scratch.
3. A slight algebraic manipulation as been made to the Borwein series. It is now able to handle z with arbitrarily large real part, so there is no need to have a special case for z with large real part.
4. The logsumexp calculation for the log of sinpi must be done manually, since logsumexp is not something that's available currently in `xsf`. I added simple (private) implementation of complex `exppi` ($exp(\pi z)$) for use here.

I've added test cases comparing with the `mpmath` reference implementation.

#### Additional information
<!--Any additional information you think is important.-->
Some things that may need attention:

The Riemann zeta function and the Hurwitz Zeta function are currently combined into a single function `zeta` which dispatches to the ufunc `_riemann_zeta` when an argument `q=None`, and otherwise dispatches to a ufunc `_zeta` for the Hurwitz zeta function. I've only added complex support for `_riemann_zeta`, and `zeta` will error out when called on complex input when `q` is not `None`. Having accepted input types depend on a keyword argument isn't ideal, but I think it's acceptable here, and complex support for Hurwitz Zeta is on the way (@mdhaber has written a reference implementation for this too). We could However change this so that complex inputs are still accepted when `q` is not `None`, but `NaN` is returned when off the real line. This is similar to what currently happens for inputs `x < 1`. They are only supported for `q = None`, but not for any other `q`.

Is there a better place for the Python reference implementation than a colab notebook? Perhaps we could add reference implementations to the `xsf` repo? I think it's fine not to worry about this now though.





Problem statement:
scipy.special.zeta and zetac don't take complex variables
The scipy.special.zeta and scipy.special.zetac don't take complex inputs.

A previously closed issue on this argued that the complex use case is not supported. 

However, there isn't really much of a use case for zeta() without complex variables.

This  limitation is like having a euclidian distance function that only supported one-dimensional inputs.

Hints from issue discussion:
We use the real-domain implementation of `zeta` for the [zeta/Zipf probability distribution](https://en.wikipedia.org/wiki/Zeta_distribution), which has its own use cases. Adding complex-domain support would, of course, be valuable, but it's not like the current implementation is useless.
[Here is the previous discussion.](https://github.com/scipy/scipy/issues/9182#issuecomment-416663307)

Can you clarify what you want to compute in the complex domain? Is floating-point accuracy suitable for that computation? In the previous discussion, it was not, which is why the issue was closed.
AFAICT, neither Boost nor GSL computes `zeta` over the complex domain either. The only hits I see for libraries that compute `zeta` over the complex domain are arbitrary-precision libraries. If that's what you're interested in, please see [`mpmath`](https://mpmath.org/).
To provide an example where one needs to evaluate the zeta function with complex arguments. In quantum physics, the correlation function of an Ohmic environment is given by
```math
C(t) = c \Bigl[ \zeta\Bigl( s+1, \frac{1 + \beta\omega_c - \mathrm i \omega_c t}{\beta\omega_c} \Bigr) + \zeta\Bigl( s+1, \frac{1 + \mathrm i \omega_c t}{\beta\omega_c} \Bigr) \Bigr]
```
where $c$ is some prefactor, $\beta$ and $\omega_c$ the inverse temperature and the cutoff frequency of the environment and $s \geq 0$ a parameter determining its low-frequency behavior (see e.g. [these lecture notes](https://www1.itp.tu-berlin.de/brandes/public_html/publications/notes.pdf)).

Not being able to evaluate this expression with scipy is unfortunate. Floating point accuracy would be sufficient.
I have this as one of the planned deliverables for a small development grant I’m working on. eta for getting something out is a few months. It’s probably going to require using extended (not arbitrary) precision internally.

Commit message: Update scipy/special/xsf/zeta.h [skip ci]

Co-authored-by: Matt Haberland <mhaberla@calpoly.edu>

Patch to review:
diff --git a/scipy/special/_basic.py b/scipy/special/_basic.py
index 7a59257d3c1d..a188ea677489 100644
--- a/scipy/special/_basic.py
+++ b/scipy/special/_basic.py
@@ -3393,10 +3393,14 @@ def zeta(x, q=None, out=None):
 
     Parameters
     ----------
-    x : array_like of float
-        Input data, must be real
+    x : array_like of float or complex.
+        Input data
     q : array_like of float, optional
-        Input data, must be real.  Defaults to Riemann zeta.
+        Input data, must be real.  Defaults to Riemann zeta. When `q` is
+        ``None``, complex inputs `x` are supported. If `q` is not ``None``,
+        then currently only real inputs `x` with ``x >= 1`` are supported,
+        even when ``q = 1.0`` (corresponding to the Riemann zeta function).
+        
     out : ndarray, optional
         Output array for the computed values.
 
@@ -3420,6 +3424,10 @@ def zeta(x, q=None, out=None):
     see [dlmf]_ for details. The Riemann zeta function corresponds to
     the case when ``q = 1``.
 
+    For complex inputs with ``q = None``, points with
+    ``abs(z.imag) > 1e9`` and ``0 <= abs(z.real) < 2.5`` are currently not
+    supported due to slow convergence causing excessive runtime.
+
     References
     ----------
     .. [dlmf] NIST, Digital Library of Mathematical Functions,
@@ -3438,6 +3446,11 @@ def zeta(x, q=None, out=None):
     >>> zeta(4), np.pi**4/90
     (1.0823232337111381, 1.082323233711138)
 
+    First nontrivial zero:
+
+    >>> zeta(0.5 + 14.134725141734695j)
+    0 + 0j
+
     Relation to the `polygamma` function:
 
     >>> m = 3
diff --git a/scipy/special/_special_ufuncs.cpp b/scipy/special/_special_ufuncs.cpp
index edb55b0c694c..88f33c505129 100644
--- a/scipy/special/_special_ufuncs.cpp
+++ b/scipy/special/_special_ufuncs.cpp
@@ -216,7 +216,8 @@ PyMODINIT_FUNC PyInit__special_ufuncs() {
     PyModule_AddObjectRef(_special_ufuncs, "_sinpi", _sinpi);
 
     PyObject *_zeta = xsf::numpy::ufunc(
-        {static_cast<xsf::numpy::ff_f>(xsf::zeta), static_cast<xsf::numpy::dd_d>(xsf::zeta)}, "_zeta", _zeta_doc);
+	{static_cast<xsf::numpy::ff_f>(xsf::zeta), static_cast<xsf::numpy::Ff_F>(xsf::zeta),
+	 static_cast<xsf::numpy::dd_d>(xsf::zeta), static_cast<xsf::numpy::Dd_D>(xsf::zeta)}, "_zeta", _zeta_doc);
     PyModule_AddObjectRef(_special_ufuncs, "_zeta", _zeta);
 
     PyObject *airy =
@@ -847,7 +848,9 @@ PyMODINIT_FUNC PyInit__special_ufuncs() {
     PyModule_AddObjectRef(_special_ufuncs, "rgamma", rgamma);
 
     PyObject *_riemann_zeta = xsf::numpy::ufunc(
-        {static_cast<xsf::numpy::d_d>(xsf::riemann_zeta), static_cast<xsf::numpy::f_f>(xsf::riemann_zeta)},
+	 {static_cast<xsf::numpy::d_d>(xsf::riemann_zeta), static_cast<xsf::numpy::D_D>(xsf::riemann_zeta),
+	 static_cast<xsf::numpy::f_f>(xsf::riemann_zeta), static_cast<xsf::numpy::F_F>(xsf::riemann_zeta)
+	},
         "_riemann_zeta", _riemann_zeta_doc);
     PyModule_AddObjectRef(_special_ufuncs, "_riemann_zeta", _riemann_zeta);
 
diff --git a/scipy/special/tests/test_zeta.py b/scipy/special/tests/test_zeta.py
index 82b3245cac53..d05ff6c99a0d 100644
--- a/scipy/special/tests/test_zeta.py
+++ b/scipy/special/tests/test_zeta.py
@@ -1,5 +1,8 @@
+import scipy
 import scipy.special as sc
 import numpy as np
+import pytest
+
 from numpy.testing import assert_equal, assert_allclose
 
 
@@ -47,3 +50,255 @@ def test_riemann_zeta_avoid_overflow():
     s = -260.00000000001
     desired = -5.6966307844402683127e+297  # Computed with Mpmath
     assert_allclose(sc.zeta(s), desired, atol=0, rtol=5e-14)
+
+# Some of the test cases below fail for intel compilers
+cpp_compiler = scipy.__config__.CONFIG["Compilers"]["c++"]["name"]
+using_intel_compiler = cpp_compiler == "intel-cl"
+
+@pytest.mark.parametrize(
+    "z, desired, rtol",
+    [
+        ## Test cases taken from mpmath with the script:
+
+        # import numpy as np
+        # import scipy.stats as stats
+
+        # from mpmath import mp
+
+        # # seed = np.random.SeedSequence().entropy
+        # seed = 154689806791763421822480125722191067828
+        # rng = np.random.default_rng(seed)
+        # default_rtol = 1e-13
+
+        # # A small point in each quadrant outside of the critical strip
+        # cases = []
+        # for x_sign, y_sign in [1, 1], [1, -1], [-1, 1], [-1, -1]:
+        #     x = x_sign * rng.uniform(2, 8)
+        #     y = y_sign * rng.uniform(2, 8)
+        #     z = x + y*1j
+        #     reference = complex(mp.zeta(z))
+        #     cases.append((z, reference, default_rtol))
+
+        # # Moderately large imaginary part in each quadrant outside of critical strip
+        # for x_sign, y_sign in [1, 1], [1, -1], [-1, 1], [-1, -1]:
+        #     x = x_sign * rng.uniform(2, 8)
+        #     y = y_sign * rng.uniform(50, 80)
+        #     z = x + y*1j
+        #     reference = complex(mp.zeta(z))
+        #     cases.append((z, reference, default_rtol))
+
+        # # points in critical strip
+        # x = rng.uniform(0.0, 1.0, size=5)
+        # y = np.exp(rng.uniform(0, 5, size=5))
+        # z = x + y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     cases.append((complex(t), reference, default_rtol))
+        # z = x - y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     cases.append((complex(t), reference, default_rtol))
+
+        # # Near small trivial zeros
+        # x = np.array([-2, -4, -6, -8])
+        # y = np.array([1e-15, -1e-15])
+        # x, y = np.meshgrid(x, y)
+        # x, y = x.ravel(), y.ravel()
+        # z = x + y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     cases.append((complex(t), reference, 1e-7))
+
+        # # Some other points near real axis
+        # x = np.array([-0.5, 0, 0.2, 0.75])
+        # y = np.array([1e-15, -1e-15])
+        # x, y = np.meshgrid(x, y)
+        # x, y = x.ravel(), y.ravel()
+        # z = x + y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     cases.append((complex(t), reference, 1e-7))
+
+        # # Moderately large real part
+        # x = np.array([49.33915930750887, 50.55805244181687])
+        # y = rng.uniform(20, 100, size=3)
+        # x, y = np.meshgrid(x, y)
+        # x, y = x.ravel(), y.ravel()
+        # z = x + y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     cases.append((complex(t), reference, default_rtol))
+
+        # # Very large imaginary part
+        # x = np.array([0.5, 34.812847097948854, 50.55805244181687])
+        # y = np.array([1e6, -1e6])
+        # x, y = np.meshgrid(x, y)
+        # x, y = x.ravel(), y.ravel()
+        # z = x + y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     rtol = 1e-7 if t.real == 0.5 else default_rtol
+        #     cases.append((complex(t), reference, rtol))
+        #
+        # # Naive implementation of reflection formula suffers internal overflow
+        # x = -rng.uniform(200, 300, 3)
+        # y = np.array([rng.uniform(10, 30), -rng.uniform(10, 30)])
+        # x, y = np.meshgrid(x, y)
+        # x, y = x.ravel(), y.ravel()
+        # z = x + y*1j
+        # for t in z:
+        #     reference = complex(mp.zeta(t))
+        #     cases.append((complex(t), reference, default_rtol))
+        #
+        # A small point in each quadrant outside of the critical strip
+        ((3.12838509346655+7.111085974836645j),
+         (1.0192654793474945+0.08795174413289127j),
+         1e-13),
+        ((7.06791362314716-7.219497492626728j),
+         (1.0020740683598117-0.006752725913243711j),
+         1e-13),
+        ((-6.806227077655519+2.724411451005281j),
+         (0.06312488213559667-0.061641496333765956j),
+         1e-13),
+        ((-3.0170751511621026-6.3686522550665945j),
+         (-0.10330747857150148-1.214541994832571j),
+         1e-13),
+        # Moderately large imaginary part in each quadrant outside of critical strip
+        ((6.133994402212294+60.03091448000761j),
+         (0.9885701843417336+0.009636925981078128j),
+         1e-13),
+        ((6.17268142822657-64.74883149743795j),
+         (1.0080474225840865+0.012032804974965354j),
+         1e-13),
+        ((-3.462191939791879+76.16258975567534j),
+         (18672.072070850158+2908.5104826247184j),
+         1e-13),
+        ((-6.955735216531752-74.75791554155748j),
+         (-77672258.72276545+71625206.0401107j),
+         1e-13),
+        # Points in critical strip
+        ((0.4088038289823922+1.4596830498094384j),
+         (0.3032837969400845-0.47272237994110344j),
+         1e-13),
+        ((0.9673493951209633+4.918968547259143j),
+         (0.7488756907431944+0.17281553371482428j),
+         1e-13),
+        ((0.8692482679977754+66.6142398421354j),
+         (0.5831942469552066-0.26848904799062334j),
+         1e-13),
+        ((0.42771847720003764+21.783747851715468j),
+         (0.4767032638444329+0.6898148744603123j),
+         1e-13),
+        ((0.20479494678428956+33.17656449538932j),
+         (-0.6983038977487848+0.18060923618150224j),
+         1e-13),
+        ((0.4088038289823922-1.4596830498094384j),
+         (0.3032837969400845+0.47272237994110344j),
+         1e-13),
+        ((0.9673493951209633-4.918968547259143j),
+         (0.7488756907431944-0.17281553371482428j),
+         1e-13),
+        ((0.8692482679977754-66.6142398421354j),
+         (0.5831942469552066+0.26848904799062334j),
+         1e-13),
+        ((0.42771847720003764-21.783747851715468j),
+         (0.4767032638444329-0.6898148744603123j),
+         1e-13),
+        ((0.20479494678428956-33.17656449538932j),
+         (-0.6983038977487848-0.18060923618150224j),
+         1e-13),
+        # Near small trivial zeros
+        ((-2+1e-15j), (3.288175809370978e-32-3.0448457058393275e-17j), 1e-07),
+        ((-4+1e-15j), (-2.868707923051182e-33+7.983811450268625e-18j), 1e-07),
+        ((-6+1e-15j), (-1.7064292323640116e-34-5.8997591435159376e-18j), 1e-07),
+        ((-8+1e-15j), (2.5060859548261706e-33+8.316161985602247e-18j), 1e-07),
+        ((-2-1e-15j), (3.288175809371319e-32+3.0448457058393275e-17j), 1e-07),
+        ((-4-1e-15j), (-2.8687079230520114e-33-7.983811450268625e-18j), 1e-07),
+        ((-6-1e-15j), (-1.70642923235801e-34+5.8997591435159376e-18j), 1e-07),
+        ((-8-1e-15j), (2.5060859548253293e-33-8.316161985602247e-18j), 1e-07),
+        # Some other points near real axis
+        ((-0.5+1e-15j), (-0.20788622497735457-3.608543395999408e-16j), 1e-07),
+        (1e-15j, (-0.5-9.189384239689193e-16j), 1e-07),
+        ((0.2+1e-15j), (-0.7339209248963406-1.4828001150329085e-15j), 1e-07),
+        ((0.75+1e-15j), (-3.4412853869452227-1.5924832114302393e-14j), 1e-13),
+        ((-0.5-1e-15j), (-0.20788622497735457+3.608543395999408e-16j), 1e-07),
+        (-1e-15j, (-0.5+9.189387416062746e-16j), 1e-07),
+        ((0.2-1e-15j), (-0.7339209248963406+1.4828004007675122e-15j), 1e-07),
+        ((0.75-1e-15j), (-3.4412853869452227+1.5924831974403957e-14j), 1e-13),
+        # Moderately large real part
+        ((49.33915930750887+53.213478698903955j),
+         (1.0000000000000009+1.0212452494616078e-15j),
+         1e-13),
+        ((50.55805244181687+53.213478698903955j),
+         (1.0000000000000004+4.387394180390787e-16j),
+         1e-13),
+        ((49.33915930750887+40.6366015728302j),
+         (0.9999999999999986-1.502268709924849e-16j),
+         1e-13),
+        ((50.55805244181687+40.6366015728302j),
+         (0.9999999999999994-6.453929613571651e-17j),
+         1e-13),
+        ((49.33915930750887+85.83555435273925j),
+         (0.9999999999999987-2.7014400611995846e-16j),
+         1e-13),
+        ((50.55805244181687+85.83555435273925j),
+         (0.9999999999999994-1.160571605555322e-16j),
+         1e-13),
+        # Very large imaginary part
+        ((0.5+1e6j), (0.0760890697382271+2.805102101019299j), 1e-07),
+        ((34.812847097948854+1e6j),
+         (1.0000000000102545+3.150848654056419e-11j),
+         1e-13),
+        ((50.55805244181687+1e6j),
+         (1.0000000000000002+5.736517078070873e-16j),
+         1e-13),
+        ((0.5-1e6j), (0.0760890697382271-2.805102101019299j), 1e-07),
+        ((34.812847097948854-1e6j),
+         (1.0000000000102545-3.150848654056419e-11j),
+         1e-13),
+        ((50.55805244181687-1e6j),
+         (1.0000000000000002-5.736517078070873e-16j),
+         1e-13),
+        # Naive implementation of reflection formula suffers internal overflow
+        pytest.param(
+            (-217.40285743524163+13.992648136816397j),
+            (-6.012818500554211e+249-1.926943776932387e+250j),
+            5e-13,
+            marks=pytest.mark.xfail(
+                condition=using_intel_compiler,
+                reason="Fails with intel compiler.",
+            ),
+        ),
+        pytest.param(
+            (-237.71710702931668+13.992648136816397j),
+            (-8.823803086106129e+281-5.009074181335139e+281j),
+            1e-13,
+            marks=pytest.mark.xfail(
+                condition=using_intel_compiler,
+                reason="Fails with intel compiler.",
+            ),
+        ),
+        ((-294.86605461349745+13.992648136816397j), (-np.inf+np.inf*1j), 1e-13),
+        pytest.param(
+            (-217.40285743524163-16.147667799398363j),
+            (-5.111612904844256e+251-4.907132127666742e+250j),
+            5e-13,
+            marks=pytest.mark.xfail(
+                condition=using_intel_compiler,
+                reason="Fails with intel compiler.",
+            ),
+        ),
+        pytest.param(
+            (-237.71710702931668-16.147667799398363j),
+            (-1.3256112779883167e+283-2.253002003455494e+283j),
+            5e-13,
+            marks=pytest.mark.xfail(
+                condition=using_intel_compiler,
+                reason="Fails with intel compiler.",
+            ),
+        ),
+        ((-294.86605461349745-16.147667799398363j), (np.inf-np.inf*1j), 1e-13),
+    ]
+)
+def test_riemann_zeta_complex(z, desired, rtol):
+    assert_allclose(sc.zeta(z), desired, rtol=rtol)
diff --git a/scipy/special/xsf/numpy.h b/scipy/special/xsf/numpy.h
index 0e713d4c4f94..57f9b29392c2 100644
--- a/scipy/special/xsf/numpy.h
+++ b/scipy/special/xsf/numpy.h
@@ -176,6 +176,8 @@ namespace numpy {
     using ld_d = double (*)(long int, double);
     using lF_F = cfloat (*)(long int, cfloat);
     using lD_D = cdouble (*)(long int, cdouble);
+    using Dd_D = cdouble (*) (cdouble, double);
+    using Ff_F = cfloat (*) (cfloat, float);
 
     // autodiff, 2 inputs, 1 output
     using autodiff0_if_f = autodiff0_float (*)(int, autodiff0_float);
diff --git a/scipy/special/xsf/zeta.h b/scipy/special/xsf/zeta.h
index 443d887b0185..3434159637af 100644
--- a/scipy/special/xsf/zeta.h
+++ b/scipy/special/xsf/zeta.h
@@ -1,10 +1,371 @@
+/* Complex riemann-zeta function implementation based on Python implementation
+ * written by Matt Haberland (@mdhaber) in:
+ * https://colab.research.google.com/drive/1zMDSAJlXCLRqMMtJ0e9nDGjQ8iZCnmn5?usp=sharing
+ */
+
 #pragma once
 
+#include "config.h"
+#include "error.h"
+#include "gamma.h"
+#include "trig.h"
+
+#include "cephes/const.h"
 #include "cephes/zeta.h"
 #include "cephes/zetac.h"
 
 namespace xsf {
 
+namespace detail {
+
+    /* Log of absolute value of expansion coefficients for Euler-Maclaurin
+     * summation formula. log(|B2k / (2k)!|)
+     *
+     * See https://en.wikipedia.org/wiki/Riemann_zeta_function#Numerical_algorithms
+     *
+     * Generated with the script
+     *
+     * import numpy as np
+     * from mpmath import mp
+     *
+     * mp.dps = 10000
+     *
+     * results = []
+     * for k in range(51):
+     *     results.append(
+     *         float(
+     *             mp.log(abs(mp.bernoulli(2*k)/(mp.factorial(2*k))))
+     *         )
+     *     )
+     */
+    constexpr double zeta_em_log_abs_coeff_lookup[] = {
+	0.0,
+	-2.4849066497880004,
+	-6.579251212010101,
+	-10.31692083029347,
+	-14.005800284407405,
+	-17.68462940266784,
+	-21.361131560073222,
+	-25.037070502911423,
+	-28.712870599846948,
+	-32.388636197522295,
+	-36.06439319366539,
+	-39.740148041995184,
+	-43.41590235365616,
+	-47.091656531181485,
+	-50.76741067517639,
+	-54.44316481078909,
+	-58.11891894430628,
+	-61.79467307729959,
+	-65.47042721016194,
+	-69.14618134299154,
+	-72.82193547581296,
+	-76.49768960863234,
+	-80.1734437414512,
+	-83.84919787426993,
+	-87.52495200708863,
+	-91.20070613990733,
+	-94.87646027272602,
+	-98.55221440554472,
+	-102.2279685383634,
+	-105.9037226711821,
+	-109.57947680400078,
+	-113.25523093681947,
+	-116.93098506963817,
+	-120.60673920245685,
+	-124.28249333527555,
+	-127.95824746809424,
+	-131.63400160091294,
+	-135.30975573373163,
+	-138.9855098665503,
+	-142.661263999369,
+	-146.3370181321877,
+	-150.0127722650064,
+	-153.6885263978251,
+	-157.36428053064375,
+	-161.04003466346245,
+	-164.71578879628115,
+	-168.39154292909984,
+	-172.06729706191854,
+	-175.74305119473723,
+	-179.4188053275559,
+	-183.0945594603746
+    };
+
+    // Complex log of expansion coefficients for Euler-Maclaurin summation formula.
+    XSF_HOST_DEVICE inline std::complex<double> zeta_em_log_coeff(std::size_t n) {
+	std::complex<double> J(0.0, 1.0);
+	std::complex<double> result;
+	if (n < 50) {
+	    result = zeta_em_log_abs_coeff_lookup[n];
+	} else {
+	    /* Asymptotic formula
+	     * Uses https://dlmf.nist.gov/24.11#E1 to approximate B_{2n} and
+	     * Stirling's approximation for (2n)!.
+	     */
+	    result = std::log(2.0) - 2.0*n*std::log(2*M_PI);
+	}
+	if (n % 2 == 0) {
+	    /* B_{2n}/(2n)! is negative for even n. This contributes a term
+	     * pi*i when taking the log. */
+	    result += M_PI * J;
+	}
+	return result;
+    }
+
+    /* Compute riemann_zeta for complex input z using the Euler-Maclaurin formula.
+     * Computation of individual terms in expansion are logarithmized to avoid
+     * overflow. TODO: only logarithmize when necessary. */
+    XSF_HOST_DEVICE inline std::complex<double> zeta_euler_maclaurin(std::complex<double> z) {
+	if (z == 1.0) {
+	    /* Return NaN at pole since value depends on how z approaches 1.0. */
+	    return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
+	}
+	std::size_t n = static_cast<std::size_t>(std::max(std::abs(z.imag()) / 4.0, 50.0));
+	std::size_t m = n;
+	std::complex<double> result = 0.0;
+	for (std::size_t i = 1; i < n; i++) {
+	    std::complex<double> term = std::pow(static_cast<double>(i), -z);
+	    result += term;
+	    // When z.real() > 1, series converges and we can consider early termination
+	    if (z.real() > 1 && std::abs(term) / std::abs(result) <= std::numeric_limits<double>::epsilon()) {
+	     	return result;
+	    }
+	}
+	double N = static_cast<double>(n);
+	std::complex<double> b = std::pow(n, -z);
+	result += b * (0.5 + N / (z - 1.0));
+	std::complex<double> log_poch = std::log(z);
+	std::complex<double> log_factor = std::log(b) - std::log(N);
+	for (std::size_t i = 1; i <= m; i++) {
+	    std::complex<double> term = std::exp(zeta_em_log_coeff(i) + log_factor + log_poch);
+	    result += term;
+	    if (std::abs(term)/std::abs(result) <= std::numeric_limits<double>::epsilon()) {
+		return result;
+	    }
+	    log_poch += std::log(z + static_cast<double>(2*i - 1)) + std::log(z + static_cast<double>(2*i));
+	    log_factor -= 2*std::log(N);
+	}
+	// Euler-maclaurin absolute error estimate.
+	double error;
+	error = std::abs(std::exp(zeta_em_log_coeff(m + 1) + log_factor + log_poch));
+	error *= std::abs((z + 2.0*(m+1) + 1.0)/(z.real() + 2.0*(m+1) + 1.0));
+	// convert to relative error estimate
+	error /= std::abs(result);
+	if (error > 1e-8) {
+	    if (error > 1e-1) {
+		/* If error estimate predicts we don't even get 1 digit of precision, return NaN
+		 * and signal no result */
+		set_error("zeta", SF_ERROR_NO_RESULT, NULL);
+		return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
+	    }
+	    // Signal reduced precision.
+	    set_error("zeta", SF_ERROR_LOSS, NULL);
+	}
+	return result;
+    }
+
+    /* Lookup table of coefficients for Algorithm 2 from Borwein 1995
+     * Borwein, Peter B.. “An efficient algorithm for the Riemann zeta function.” (1995).
+     *
+     * Stores coefficients as dk / dn, where dn is the final coefficient.
+     *
+     * Generated with the Python script:
+     *
+     * import numpy as np
+     * import math
+     * from mpmath import mp
+     *
+     * mp.dps = 1000
+     * n = 50
+     *
+     * coeffs = []
+     * S = mp.zero
+     * for i in range(n + 1):
+     *     num = math.factorial(n + i - 1) * 4**i
+     *     den = math.factorial(n - i) * math.factorial(2*i)
+     *     S += mp.mpf(num) / mp.mpf(den)
+     *     coeffs.append(S*n)
+     *
+     * dn = coeffs[-1]
+     * coeffs = [float(dk/dn) for dk in coeffs[:-1]]
+     * coeffs = np.asarray(coeffs)
+     */
+    constexpr double zeta_borwein_coeff[] = {
+	1.0555078361382878e-38,
+	5.278594688527578e-35,
+	4.4014687322044963e-32,
+	1.467453546497519e-29,
+	2.617862196688831e-27,
+	2.900097799958025e-25,
+	2.184440361492933e-23,
+	1.1890977312913296e-21,
+	4.8871396166872276e-20,
+	1.5672253698802734e-18,
+	4.022931234264572e-17,
+	8.435973533351745e-16,
+	1.469296379914116e-14,
+	2.1548747054571902e-13,
+	2.6919530537535124e-12,
+	2.8925409162768484e-11,
+	2.6957505693699856e-10,
+	2.194772239130839e-09,
+	1.57078229370057e-08,
+	9.936187220749133e-08,
+	5.581721578217702e-07,
+	2.796271112037765e-06,
+	1.253886254275813e-05,
+	5.049261002939051e-05,
+	0.00018312884459703666,
+	0.0005997690328552426,
+	0.0017780501082460968,
+	0.004781802283665968,
+	0.011690432287131671,
+	0.026034302929535964,
+	0.052922982472754856,
+	0.09842471411648858,
+	0.16789610796540344,
+	0.26350429194768626,
+	0.3819442810600665,
+	0.5137731385068898,
+	0.645292538541866,
+	0.7625449322050362,
+	0.8556063057019102,
+	0.921056062886525,
+	0.9616100580028147,
+	0.9835905431606953,
+	0.9939187229336753,
+	0.9980782525166648,
+	0.9994930141459889,
+	0.999891478844585,
+	0.9999819091990203,
+	0.9999977981288319,
+	0.9999998260580315,
+	0.9999999933099243
+    };
+
+    /* Compute riemann_zeta for complex input z using Algorithm 2 from Borwein 1995. */
+    XSF_HOST_DEVICE inline std::complex<double> zeta_borwein(std::complex<double> z) {
+	std::complex<double> result = 0.0;
+	// Sum in reverse order because smaller terms come later.
+	for (int k = 49; k >= 0; k--) {
+	    double sign = std::pow(-1.0, k);
+	    std::complex<double> den = std::pow(k + 1, z);
+	    std::complex<double> term = sign * (zeta_borwein_coeff[k] - 1.0) / den;
+	    result += term;
+	}
+	return result * -1.0/(1.0 - std::pow(2.0, 1.0 - z));
+    }
+
+    /* Compute riemann zeta for complex z and real part >= 0 */
+    XSF_HOST_DEVICE inline std::complex<double> zeta_right_halfplane(std::complex<double> z) {
+	if (z == 1.0) {
+	    return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
+	}
+	/* Cutoff for using Euler-MacLaurin chosen based on cursory empirical search.
+	 * TODO: Choose cutoffs in a more principled way. */
+	if (z.real() < 50.0 && std::abs(z.imag()) > 50.0) {
+	    if (z.real() >= 0.0 && z.real() < 2.5 && std::abs(z.imag()) > 1e9) {
+		/* Euler-MacLaurin summation starts to take an unreasonable amount of time in this
+		 * region, so just give up and return NaN instead. */
+		set_error("zeta", SF_ERROR_NO_RESULT, NULL);
+		return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
+	    }
+	    return zeta_euler_maclaurin(z);
+	}
+	return zeta_borwein(z);
+    }
+
+    XSF_HOST_DEVICE inline std::complex<double> exppi(std::complex<double> z) {
+	    // exp(pi*z) for complex z.
+	    double x = z.real();
+	    double y = z.imag();
+	    std::complex<double> factor1(xsf::cospi(y), xsf::sinpi(y));
+	    double factor2 = std::exp(M_PI*x);
+	    return factor1 * factor2;
+	}
+
+    XSF_HOST_DEVICE inline std::complex<double> logsinpi(std::complex<double> z) {
+	/* log(sinpi(z)) using sin(z) = (exp(i*pi*z) - exp(-i*pi*z)) / 2i
+	 *
+	 * No attempt is made to choose any particular branch of the logarithm.
+	 * This is an internal function and the intent is that this that the
+	 * result of log(sinpi(z)) will be added to other terms, and the sum
+	 * will then be exponentiated, making the choice of a specific branch
+	 * unnecessary.
+	 */
+	std::complex<double> result = std::log(xsf::sinpi(z));
+	// If it doesn't overflow, just do the regular calculation.
+	if (std::isfinite(result.real()) && !std::isfinite(result.imag())) {
+	    return result;
+	}
+	/* Otherwise factor before taking log. This is where we may end up
+	 * taking a branch other than the principal branch. */
+	std::complex<double> J(0.0, 1.0);
+	/* Calculating log((exp(i*pi*z) - exp(-i*pi*z)) / 2i). Factor out term
+	 * with larger magnitude before taking log. */
+	if (z.imag() > 0 ) {
+	    /* if z.imag() > 0 then, exp(-i*pi*z) has greatest magnitude. Factor it
+	     * out to get:
+	     * log(exp(-i*pi*z)*((exp(2*i*pi*z) - 1.0)/(2i)) =
+	     * log(exp(-i*pi*z)) + log((exp(2*i*pi*z) - 1.0)/(2i)) =
+	     * -i*pi*z + log((exp(2*i*pi*z) - 1.0)/(2i)) */
+	    return -J * M_PI * z + std::log((exppi(2.0 * z * J) - 1.0) / (2.0*J));
+	}
+	/* if z.imag() < 0 then, exp(i*pi*z) has greatest magnitude. Factor similarly
+	 * to above */
+	return J * M_PI * z + std::log((1.0 - exppi(-2.0 * z * J)) / (2.0*J));
+    }
+
+    /* Leading factor in reflection formula for zeta function.
+     * zeta(z) = 2^z * pi^(z-1) * sin(pi*z/2) * gamma(1 - z) * zeta(1 - z)
+     * This computes 2^z * pi^(z - 1) * sin(pi*z/2) * gamma(1 - z)
+     *
+     * Computation is logarithimized to prevent overflow.
+     * TODO: Complexify the cephes zeta_reflection implementation, which uses
+     * the lanczos approximation for the gamma function. */
+    XSF_HOST_DEVICE inline std::complex<double> zeta_reflection_factor_with_logs(std::complex<double> z) {
+ 	std::complex<double> t1 = z * M_LN2;
+	std::complex<double> t2 = (z - 1.0) * xsf::cephes::detail::LOGPI;
+	std::complex<double> t3 = logsinpi(z / 2.0);
+	std::complex<double> t4 = xsf::loggamma(1.0 - z);
+	std::complex<double> factor = std::exp(t1 + t2 + t3 + t4);
+	return factor;
+    }
+
+    XSF_HOST_DEVICE inline std::complex<double> zeta_reflection(std::complex<double> z) {
+	std::complex<double> factor = 2.0 * std::pow(2*M_PI, z - 1.0) * xsf::sinpi(z/2.0) * xsf::gamma(1.0 - z);
+	if (!std::isfinite(factor.real()) || !std::isfinite(factor.imag())) {
+	    // Try again with logs if standard calculation had overflow.
+	    factor = zeta_reflection_factor_with_logs(z);
+	}
+	std::complex<double> result = zeta_right_halfplane(1.0 - z);
+	/* zeta tends to 1.0 as real part tends to +inf. In cases where
+	 * the real part of zeta tends to -inf, then zeta(1 - z) in the
+	 * reflection formula will tend to 1.0. Factor overflows then,
+	 * factor * result below will become NaN. In this case, we just
+	 * return factor to preserve complex infinity. Only zeta(1 - z) == 1.0
+	 * is handled because this is the only practical case where we should
+	 *  expect zeta(1 - z) == x for a real number x when z is not on the
+	 * real line. */
+	return (result == 1.0) ? factor : factor * result;
+    }
+}
+
+XSF_HOST_DEVICE inline std::complex<double> riemann_zeta(std::complex<double> z) {
+    if (z.imag() == 0.0) {
+	return cephes::riemann_zeta(z.real());
+    }
+    if (z.real() >= 0.5) {
+	return detail::zeta_right_halfplane(z);
+    }
+    return detail::zeta_reflection(z);
+}
+
+XSF_HOST_DEVICE inline std::complex<float> riemann_zeta(std::complex<float> z) {
+    return static_cast<std::complex<float>>(riemann_zeta(static_cast<std::complex<double>>(z)));
+}
+
 XSF_HOST_DEVICE inline double riemann_zeta(double x) { return cephes::riemann_zeta(x); }
 
 XSF_HOST_DEVICE inline float riemann_zeta(float x) { return riemann_zeta(static_cast<double>(x)); }
@@ -13,6 +374,19 @@ XSF_HOST_DEVICE inline double zeta(double x, double q) { return cephes::zeta(x,
 
 XSF_HOST_DEVICE inline float zeta(float x, float q) { return zeta(static_cast<double>(x), static_cast<double>(q)); }
 
+XSF_HOST_DEVICE inline std::complex<double> zeta(std::complex<double> z, double q) {
+    if (z.imag() == 0.0) {
+	return zeta(z.real(), q);
+    }
+    // Complex input for Hurwitz Zeta is not currently supported.
+    set_error("zeta", SF_ERROR_DOMAIN, NULL);
+    return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
+}
+
+XSF_HOST_DEVICE inline std::complex<float> zeta(std::complex<float> z, float q) {
+    return static_cast<std::complex<float>>(zeta(static_cast<std::complex<double>>(z), static_cast<float>(q)));
+}
+
 XSF_HOST_DEVICE inline double zetac(double x) { return cephes::zetac(x); }
 
 XSF_HOST_DEVICE inline float zetac(float x) { return zetac(static_cast<double>(x)); }


Code review comments:
diff_hunk: @@ -1,10 +1,371 @@
+/* Complex riemann-zeta function implementation based on Python implementation
+ * written by Matt Haberland (@mdhaber) in:
+ * https://colab.research.google.com/drive/1zMDSAJlXCLRqMMtJ0e9nDGjQ8iZCnmn5?usp=sharing
+ */
+
 #pragma once
 
+#include "config.h"
+#include "error.h"
+#include "gamma.h"
+#include "trig.h"
+
+#include "cephes/const.h"
 #include "cephes/zeta.h"
 #include "cephes/zetac.h"
 
 namespace xsf {
 
+namespace detail {
+
+    /* Log of absolute value of expansion coefficients for Euler-Maclaurin
+     * summation formula. log(|B2k / (2k)!|)
+     *
+     * See https://en.wikipedia.org/wiki/Riemann_zeta_function#Numerical_algorithms
+     *
+     * Generated with the script
+     *
+     * import numpy as np
+     * from mpmath import mp
+     *
+     * mp.dps = 10000
+     *
+     * results = []
+     * for k in range(51):
+     *     results.append(
+     *         float(
+     *             mp.log(abs(mp.bernoulli(2*k)/(mp.factorial(2*k))))
+     *         )
+     *     )
+     */
+    constexpr double zeta_em_log_abs_coeff_lookup[] = {
+	0.0,
+	-2.4849066497880004,
+	-6.579251212010101,
+	-10.31692083029347,
+	-14.005800284407405,
+	-17.68462940266784,
+	-21.361131560073222,
+	-25.037070502911423,
+	-28.712870599846948,
+	-32.388636197522295,
+	-36.06439319366539,
+	-39.740148041995184,
+	-43.41590235365616,
+	-47.091656531181485,
+	-50.76741067517639,
+	-54.44316481078909,
+	-58.11891894430628,
+	-61.79467307729959,
+	-65.47042721016194,
+	-69.14618134299154,
+	-72.82193547581296,
+	-76.49768960863234,
+	-80.1734437414512,
+	-83.84919787426993,
+	-87.52495200708863,
+	-91.20070613990733,
+	-94.87646027272602,
+	-98.55221440554472,
+	-102.2279685383634,
+	-105.9037226711821,
+	-109.57947680400078,
+	-113.25523093681947,
+	-116.93098506963817,
+	-120.60673920245685,
+	-124.28249333527555,
+	-127.95824746809424,
+	-131.63400160091294,
+	-135.30975573373163,
+	-138.9855098665503,
+	-142.661263999369,
+	-146.3370181321877,
+	-150.0127722650064,
+	-153.6885263978251,
+	-157.36428053064375,
+	-161.04003466346245,
+	-164.71578879628115,
+	-168.39154292909984,
+	-172.06729706191854,
+	-175.74305119473723,
+	-179.4188053275559,
+	-183.0945594603746
+    };
+
+    // Complex log of expansion coefficients for Euler-Maclaurin summation formula.
+    XSF_HOST_DEVICE inline std::complex<double> zeta_em_log_coeff(std::size_t n) {
+	std::complex<double> J(0.0, 1.0);
+	std::complex<double> result;
+	if (n < 50) {
+	    result = zeta_em_log_abs_coeff_lookup[n];
+	} else {
+	    /* Asymptotic formula
+	     * Uses https://dlmf.nist.gov/24.11#E1 to approximate B_{2n} and
+	     * Stirling's approximation for (2n)!.
+	     */
+	    result = std::log(2.0) - 2.0*n*std::log(2*M_PI);
+	}
+	if (n % 2 == 0) {
+	    /* B_{2n}/(2n)! is negative for even n. This contributes a term
+	     * pi*i when taking the log. */
+	    result += M_PI * J;
+	}
+	return result;
+    }
+
+    /* Compute riemann_zeta for complex input z using the Euler-Maclaurin formula.
+     * Computation of individual terms in expansion are logarithmized to avoid
+     * overflow. TODO: only logarithmize when necessary. */
+    XSF_HOST_DEVICE inline std::complex<double> zeta_euler_maclaurin(std::complex<double> z) {
+	if (z == 1.0) {
+	    /* Return NaN at pole since value depends on how z approaches 1.0. */
+	    return {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
+	}
+	std::size_t n = static_cast<std::size_t>(std::max(std::abs(z.imag()) / 4.0, 50.0));
+	std::size_t m = n;
+	std::complex<double> result = 0.0;
+	for (std::size_t i = 1; i < n; i++) {
+	    std::complex<double> term = std::pow(static_cast<double>(i), -z);
+	    result += term;
+	    // When z.real() > 1, series converges and we can consider early termination
+	    if (z.real() > 1 && std::abs(term) / std::abs(result) <= std::numeric_limits<double>::epsilon()) {
+	     	return result;
+	    }
+	}
+	double N = static_cast<double>(n);
+	std::complex<double> b = std::pow(n, -z);
+	result += b * (0.5 + N / (z - 1.0));
+	std::complex<double> log_poch = std::log(z);
+	std::complex<double> log_factor = std::log(b) - std::log(N);
+	for (std::size_t i = 1; i <= m; i++) {
+	    std::complex<double> term = std::exp(zeta_em_log_coeff(i) + log_factor + log_poch);
+	    result += term;
+	    if (std::abs(term)/std::abs(result) <= std::numeric_limits<double>::epsilon()) {
+		return result;
+	    }
+	    log_poch += std::log(z + static_cast<double>(2*i - 1)) + std::log(z + static_cast<double>(2*i));
+	    log_factor -= 2*std::log(N);
+	}
+	// Euler-maclaurin absolute error estimate.
+	double error;
+	error = std::abs(std::exp(zeta_em_log_coeff(m + 1) + log_factor + log_poch));
+	error *= std::abs((z + 2.0*(m+1) + 1.0)/(z.real() + 2.0*(m+1) + 1.0));
path: scipy/special/xsf/zeta.h
text: @user1:
The sum is
```math
\sum_{k=1}^m T_{k, n}(s)
```
and your loop is `i = 1; i <= m`, so the last execution is with `i = m`. 

The formula for the error is:
![{6A46A812-47FB-42DA-AAB5-FECEEED7D986}](https://github.com/user-attachments/assets/a4a8a9c0-af9f-4818-8a7e-e2570c39cd71)

So shouldn't this just be:
```suggestion
	error *= std::abs((z + 2.0*m + 1.0)/(z.real() + 2.0*m + 1.0));
```
IIUC, the indices should "just match" the formula; we shouldn't need to correct for any one offsets anymore.

@author:
Yes, good catch. Committing and restarting CI since I just restarted in anyway.