Module zeta 

The Riemann zeta function, as $\zeta(s) - 1$ with $\zeta$ built on top.

Which one is the primitive

$\zeta(s) \to 1$ fast: $\zeta(40) - 1$ is about $9.1\times10^{-13}$, already far below the mantissa of $\zeta$ itself, and by $s = 80$ the complement is 8.3e-25. So a caller who wants the complement cannot get it by subtracting: measured at s = 80, forming zeta(s) and taking away 1 is 100% wrong, and by s = 200 it returns a flat zero where the true value is 1e-61.

The Euler-Maclaurin sum below opens with the $n = 1$ term, which is that 1, so the complement comes from omitting it rather than cancelling it (exact, with no subtraction anywhere), and still carries digits at s = 700 where the value is around 1e-211. That makes zetac the primitive here and $\zeta = 1 + \zeta_c$ the derived form, the same relationship exp_m1 has to exp.

Algorithm

Euler-Maclaurin, truncated at N direct terms with bernoulli_terms correction terms:

\zeta(s) = \sum_{n=1}^{N-1} n^{-s} + \frac{N^{1-s}}{s-1} + \frac{N^{-s}}{2}
         + \sum_{k\ge1} \frac{B_{2k}}{(2k)!}\,(s)_{2k-1}\,N^{-(s+2k-1)}

The usual alternative is a table of minimax rationals over five or six intervals in s (this is what Boost does, in about a thousand lines). That is excellent scalar code and the wrong shape for a vector unit, where selecting a coefficient table per lane means either a gather or evaluating every interval and discarding all but one. Euler-Maclaurin is one straight-line expression for the whole positive axis instead.

Borwein’s accelerated eta series was the other candidate, and is the more famous one because it converges in the critical strip and over the complex plane. Measured, it needs 22 terms to match this at 10, more than twice the transcendental calls for the same answer, so it lost on cost. If complex s is ever wanted, it becomes interesting again.

Four exponentials, not nine

Every direct term is $n^{-s} = 2^{-s\log_2 n}$ with $\log_2 n$ a compile-time constant, which reads as one exp2 per term. But the Dirichlet terms factor over the primes: with $p_n = n^{-s}$ evaluated for n in 2, 3, 5, 7, the rest are products: $p_4 = p_2^2$, $p_6 = p_2p_3$, $p_8 = p_2^3$, $p_9 = p_3^2$, and $N^{-s} = p_2p_5$ needs no call of its own. Four transcendentals and five multiplies cover all of n = 2..10, at measured accuracy indistinguishable from nine separate calls (4.89e-16 against 4.41e-16).

The count is $\pi(N)$, the prime-counting function, not N, so raising N to tighten the critical strip is cheaper than it looks: N = 16 costs six, N = 20 costs eight.

The correction sum needs no transcendentals at all. $N^{-(s+2k-1)}$ is $N^{-s}$ times a constant, and $(s)_{2k-1}/(2k)!$ advances by a two-factor recurrence whose denominator is a compile-time integer, so the whole tail is one multiply-accumulate ladder over the shipped Bernoulli table.

Accuracy

Against mpmath at 40 digits, worst relative error with N = 10 and 8 correction terms: 4.4e-16 for s in [1.5, 5], 4.3e-16 for [5, 40], 2.3e-15 through the critical strip [0.1, 0.9], and 4.6e-16 approaching the pole. The strip is the weak region, and N is the lever if it ever matters.

Negative s is not reachable by adding terms. The expansion is asymptotic, and its error there gets worse with larger N (measured 3.7e-9 at N = 10, 5.4e-8 at N = 16). It takes the functional equation instead, which lands at $1 - s > 1$, back in the region where the series is at its best.

Constants

N

The truncation point of the direct sum. Ten is the knee: eight direct terms and six correction terms leave 2.4e-14, ten and eight reach 4.4e-16, and more of either buys nothing (12 and 8 measured 2.7e-16). Because the terms factor over the primes, the transcendental cost is pi(10) = 4 rather than 9.

Traits

ZetaConsts

Per-element constants: the base-2 logarithms of the primes under N. log2(pi) for the functional equation comes from FloatConsts. Declared for f32/f64. Add more as needed.

Functions

bernoulli_terms

Correction terms by precision tier. The dropped term bounds the error directly, and the series is convergent-then-asymptotic in this range, so the tiers are: 8 terms is 4.4e-16, 4 is around 1e-11, and 2 is around 1e-7, which straddles f32’s floor, where the whole tail is nearly free anyway.

zeta_core

zeta(s) (ZETAC = false) or zeta(s) - 1 (ZETAC = true), and its derivative when DERIV is set. zeta and zetac differ by a constant, so the one derivative serves both.

zeta_impl

zeta(s) (ZETAC = false) or zeta(s) - 1 (ZETAC = true), for real s.