For the failing test cases, I think it's fine to punt on the Intel oneAPI failures, and I've attempted to xfail those. I also loosened the rtol for a couple test cases.

I tried checking the ∑ T k , n part but it is still challenging to cross the i's and dot the t's. If there's anything else you're willing to do to make it look more like the reference, I'd appreciate it. Every little bit would help, like using -z * log(n) in log_factor instead of log(b) and changing the name of the loop variable.

Noted. Things should be more clear now.

Do you also want someone more familiar with the special/C++-specific stuff to review from that perspective, or are you comfortable with that as long as I've checked the algorithmic stuff?

Feel free to merge based on your algorithmic review. This is all pretty straightforward from a C++ perspective.

Playing around with this, I think it's pretty good. I think you already knew that 0 + 0j is a problem spot.

(The positive axis labels are something like 16 + np.log10(diff), where diff is the distance from the center point.)

z = (1e-15 + 1e-15j)
res = complex(scipy_zeta(z))
ref = complex(mpmath_zeta(z))
print(res, ref, abs(res-ref)/abs(ref))
# (-0.5097353193036964+0.010130110463593661j) (-0.5000000000000009-9.189385332046748e-16j) 0.02809950746540121

Similar story at 1 + 0j.

z = (1+1e-15 + 1e-15j)
res = complex(scipy_zeta(z))
ref = complex(mpmath_zeta(z))
print(res, ref, abs(res-ref)/abs(ref))
# (496745969635765.9-443048778321942.06j) (497279149540589.06-447909238514022.3j) 0.007306002718971802

At negative even integers (trivial zeros), there is only a problem on the real axis (at least imaginary offsets greater than 1e-16).

z = (-2+1e-15)
res = complex(scipy_zeta(z))
ref = complex(mpmath_zeta(z))
print(res, ref, abs(res-ref)/abs(ref))
# (-3.6806848447762256e-17+0j) (-3.3804578090538615e-17+0j) 0.08881253743746413

But the absolute error is OK, and this is a problem with the existing real implementation, not this PR. It might be worth a second look at whether we should just use the complex implementation even when the input is real. IIRC, there were problems at the trivial zeros in my implementation, so I replaced it on the entire real axis, but maybe we should just special-case the trivial zeros.

When we zoom out, though, the error is typically quite good. Consider white to mean esentially zero error: on the left, both mpmath and scipy return an infinite result, and on the right, the result of both is essentially just 1.0.

Looking near the critical strip, we get some inaccuracy near the transition between algorithms.

But typically the error near the critical strip is not bad.

It seems fine very close to the real axis.

There is a transition near 50 + 50j that looks fine.

The transitions at 50 - 50j looks about the same, and the transitions at the negative reals are invisible, probably because the error is all due to use of the log-reflection (but it's on the order of 1e-13, so nothing to worry about).

In the critical strip with imaginary part close to 1e9, the error is around 1e-7. It takes quite a long time where EM is used here, but probably not so long that people would think it hangs. You've added some protection against long execution times for real part less than 2.5, but it can start taking a long time and get inaccurate when the real and imaginary parts are greater, e.g.

z = 2.51 + 1e13*1j
res = complex(scipy_zeta(z))
ref = complex(mpmath_zeta(z))
print(res, ref, abs(res-ref)/abs(ref))
# (0.9381747827300801-0.21060417060191572j) (0.938250825152661-0.2106505908517201j) 9.26485094426602e-05

But I think this PR is about the main algorithms, not the fine tuning or special cases. We know where some additional work may be needed, but let's merge the main implementation. Thanks @steppi! (I'll merge later tonight if you're happy, too.)

Was just updating the release notes for the factorial changes, and noticed that this doesn't have an entry yet.

Was just updating the release notes for the factorial changes, and noticed that this doesn't have an entry yet.

Thanks @h-vetinari. There's a few other things I need to add release notes for too. I'll try to get them all done by the end of the month.