zeta.h – support for the zeta function

This module implements various algorithms for evaluating the Riemann zeta function and related functions. The functions provided here are mainly intended for internal use, though they may be useful to call directly in some applications where the default algorithm choices are suboptimal. Most applications should use the user-friendly functions in the fmprb and fmpcb modules (or for power series, the functions in the fmprb_poly and fmpcb_poly modules).

Integer zeta values

void zeta_apery_bsplit(fmprb_t x, long prec)

Sets x to Apery’s constant \(\zeta(3)\), computed by applying binary splitting to a hypergeometric series.

void zeta_ui_asymp(fmprb_t z, ulong s, long prec)

Assuming \(s \ge 2\), approximates \(\zeta(s)\) by \(1 + 2^{-s}\) along with a correct error bound. We use the following bounds: for \(s > b\), \(\zeta(s) - 1 < 2^{-b}\), and generally, \(\zeta(s) - (1 + 2^{-s}) < 2^{2-\lfloor 3 s/2 \rfloor}\).

void zeta_ui_euler_product(fmprb_t z, ulong s, long prec)

Computes \(\zeta(s)\) using the Euler product. This is fast only if s is large compared to the precision.

Writing \(P(a,b) = \prod_{a \le p \le b} (1 - p^{-s})\), we have \(1/\zeta(s) = P(a,M) P(M+1,\infty)\).

To bound the error caused by truncating the product at \(M\), we write \(P(M+1,\infty) = 1 - \epsilon(s,M)\). Since \(0 < P(a,M) \le 1\), the absolute error for \(\zeta(s)\) is bounded by \(\epsilon(s,M)\).

According to the analysis in [Fil1992], it holds for all \(s \ge 6\) and \(M \ge 1\) that \(1/P(M+1,\infty) - 1 \le f(s,M) \equiv 2 M^{1-s} / (s/2 - 1)\). Thus, we have \(1/(1-\epsilon(s,M)) - 1 \le f(s,M)\), and expanding the geometric series allows us to conclude that \(\epsilon(M) \le f(s,M)\).

void zeta_ui_bernoulli(fmprb_t x, ulong n, long prec)

Computes \(\zeta(n)\) for even n via the corresponding Bernoulli number.

void zeta_ui_vec_borwein(fmprb_ptr z, ulong start, long num, ulong step, long prec)

Evaluates \(\zeta(s)\) at \(\mathrm{num}\) consecutive integers s beginning with start and proceeding in increments of step. Uses Borwein’s formula ([Bor2000], [GS2003]), implemented to support fast multi-evaluation (but also works well for a single s).

Requires \(\mathrm{start} \ge 2\). For efficiency, the largest s should be at most about as large as prec. Arguments approaching LONG_MAX will cause overflows. One should therefore only use this function for s up to about prec, and then switch to the Euler product.

The algorithm for single s is basically identical to the one used in MPFR (see [MPFR2012] for a detailed description). In particular, we evaluate the sum backwards to avoid storing more than one \(d_k\) coefficient, and use integer arithmetic throughout since it is convenient and the terms turn out to be slightly larger than \(2^\mathrm{prec}\). The only numerical error in the main loop comes from the division by \(k^s\), which adds less than 1 unit of error per term. For fast multi-evaluation, we repeatedly divide by \(k^{\mathrm{step}}\). Each division reduces the input error and adds at most 1 unit of additional rounding error, so by induction, the error per term is always smaller than 2 units.

void zeta_ui_borwein_bsplit(fmprb_t x, ulong s, long prec)

Computes \(\zeta(s)\) for arbitrary \(s \ge 2\) using a binary splitting implementation of Borwein’s algorithm. This has quasilinear complexity with respect to the precision (assuming that \(s\) is fixed). We have

\[\zeta(s) = \frac{1}{d_n (1-2^{1-s})} \sum_{k=0}^{n-1} \frac{(-1)^k(d_n-d_k)}{(k+1)^s} + \gamma_n(s)\]

where

\[d_k = n \sum_{i=0}^k \frac{(n+i-1)! 4^i}{(n-i)! (2i)!}.\]

On the domain of interest, \(|\gamma_n(s)| \le 3 / (3 + \sqrt 8)^n\).

We write the summation as a system of first-order recurrences for \((s_k, d_k, t_k)\) where \(t_k = d_k - d_{k-1}\). This system is described by the matrix equation

\[\begin{split}\begin{pmatrix} s_{k+1} \\ d_{k+2} \\ t_{k+3} \end{pmatrix} = \begin{pmatrix} 1 & (-1)^k (k+1)^{-s} & 0 \\ 0 & 1 & 1 \\ 0 & 0 & u(k) \end{pmatrix} \begin{pmatrix} s_k \\ d_{k+1} \\ t_{k+2} \end{pmatrix}.\end{split}\]

We derive the binary splitting scheme by considering a product of an arbitrary pair in the chain \(M_{n-1} M_{n-2} \cdots M_1 M_0\). This gives

\[\begin{split}\begin{pmatrix} 1 & A_L & B_L \\ 0 & 1 & C_L \\ 0 & 0 & D_L \end{pmatrix} \begin{pmatrix} 1 & A_R & B_R \\ 0 & 1 & C_R \\ 0 & 0 & D_R \end{pmatrix} = \begin{pmatrix} 1 & A_L+A_R & B_R+A_L C_R+B_L D_R \\ 0 & 1 & C_R+C_L D_R \\ 0 & 0 & D_L D_R \end{pmatrix}.\end{split}\]

The next step is to clear denominators. Instead of putting the whole matrix on a common denominator, we optimize by putting \(C, D\) on a denominator \(Q_1\) (the product of denominators of \(u\)) and \(A, B\) on a common denominator \(Q_3 = Q_1 Q_2\) (where \(Q_2\) is the product of \((k+1)^s\) factors). This gives a small efficiency improvement. Thus, we have

\[\begin{split}\begin{pmatrix} 1 & \dfrac{A_L}{Q_{3L}} & \dfrac{B_L}{Q_{3L}} \\[3ex] 0 & 1 & \dfrac{C_L}{Q_{1L}} \\[3ex] 0 & 0 & \dfrac{D_L}{Q_{1L}} \end{pmatrix} \begin{pmatrix} 1 & \dfrac{A_R}{Q_{3R}} & \dfrac{B_R}{Q_{3R}} \\[3ex] 0 & 1 & \dfrac{C_R}{Q_{1R}} \\[3ex] 0 & 0 & \dfrac{D_R}{Q_{1R}} \end{pmatrix} = \begin{pmatrix} 1 & \dfrac{Q_{3L} A_R + A_L Q_{3R}}{Q_{3L} Q_{3R}} & \dfrac{Q_{3L} B_R + A_L C_R Q_{2R} + B_L D_R Q_{2R}}{Q_{3L} Q_{3R}} \\[3ex] 0 & 1 & \dfrac{Q_{1L} C_R + C_L D_R}{Q_{1L} Q_{1R}} \\[3ex] 0 & 0 & \dfrac{D_L D_R}{Q_{1L} Q_{1R}} \end{pmatrix}.\end{split}\]

In the final matrix, we note that \(A / Q_3 = \sum_k (-1)^k (k+1)^{-s}\), and \(C / Q_1 = d_n\). Thus \((1 / d_n) \sum_k (-1)^k (k+1)^{-s} (d_n - d_k)\) is given by \(A/Q_3 - (B/Q_3) / (C/Q_1) = (A C - B Q_1) / (Q_3 C)\).

void zeta_ui(fmprb_t x, ulong s, long prec)

Computes \(\zeta(s)\) for nonnegative integer \(s \ne 1\), automatically choosing an appropriate algorithm.

void zeta_ui_vec(fmprb_ptr x, ulong start, long num, long prec)

void zeta_ui_vec_even(fmprb_ptr x, ulong start, long num, long prec)

void zeta_ui_vec_odd(fmprb_ptr x, ulong start, long num, long prec)

Computes \(\zeta(s)\) at num consecutive integers (respectively num even or num odd integers) beginning with \(s = \mathrm{start} \ge 2\), automatically choosing an appropriate algorithm.

Euler-Maclaurin summation

void zeta_series_em_sum(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, int deflate, ulong N, ulong M, long d, long prec)

void zeta_series(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, int deflate, long d, long prec)

Evaluates the truncated Euler-Maclaurin sum of order \(N, M\) for the length-d truncated Taylor series of the Hurwitz zeta function \(\zeta(s,a)\) at \(s\), using a working precision of prec bits. With \(a = 1\), this gives the usual Riemann zeta function.

If deflate is nonzero, \(\zeta(s,a) - 1/(s-1)\) is evaluated (which permits series expansion at \(s = 1\)).

The fmpcb_zeta_series function chooses default values for \(N, M\) using fmpcb_zeta_series_em_choose_param, targeting an absolute truncation error of \(2^{-\operatorname{prec}}\).

The Euler-Maclaurin (EM) formula states that

\[\sum_{k=N}^U f(k) = \int_N^U f(t) dt + \frac{1}{2} \left(f(N) + f(U)\right)\]\[ + \sum_{k=1}^{M} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(U) - f^{(2k-1)}(N) \right)\]\[ - \int_N^U \frac{\tilde B_{2M}(t)}{(2M)!} f^{(2M)}(t) dt\]

where \(f\) is a sufficiently differentiable function (for example, analytic), \(B_n\) is a Bernoulli number, and \(\tilde B_n(t) = B_n(t-\lfloor t\rfloor)\) is a periodic Bernoulli polynomial. If \(f\) decreases sufficiently rapidly, the formula remains valid after letting \(U \to \infty\).

To evaluate the Hurwitz zeta function, we set \(f(k) = (a + k)^{-s}\), giving \(\zeta(s,a) = \sum_{k=0}^{N-1} f(k) + \sum_{k=N}^{\infty} f(k)\), where EM summation is applied to the right sum. By choosing \(M\) and \(N\) large enough, and taking the standard logarithm branch cut, the EM formula gives an analytic continuation of \(\zeta(s,a)\) to all \(a, s \in \mathbb{C}\) (except for poles at \(s = 1\) and \(\mathrm{Re}(s) > 0, a = 0, -1, -2, \ldots\)). In order to evaluate derivatives with respect to \(s\) of \(\zeta(s,a)\), we substitute \(s \to s + x \in \mathbb{C}[[x]]\).

We choose \(N\) such that \(\Re(a+N) > 0\). Then the first integral is well-defined for \(s\) with \(\Re(s) > 1\) and has the closed form

\[\int_N^{\infty} f(t) dt = \int_N^{\infty} (a + t)^{-s}dt = \frac{(a+N)^{1-s}}{s-1},\]

providing analytic continuation of this term with respect to \(s\). Removing the singularity from this term also conveniently allows us to evaluate derivatives of \(\zeta(s,a) - 1/(s-1)\) at \(s = 1\).

The derivatives of \(f(k)\) are given by

\[f^{(r)}(k) = \frac{(-1)^r (s)_{r}}{(a+k)^{s+r}}\]

where \((s)_{r} = s (s+1) \cdots (s+r-1)\) denotes a rising factorial. Thus, the remainder integral becomes

\[R(s) = \int_N^{\infty} \frac{\tilde B_{2M}(t)}{(2M)!} \frac{(s)_{2M}}{(a+t)^{s+2M}} dt,\]

valid when \(\Re(a+N) > 0\) and \(\Re(s+2M-1) > 0\). We will use the stronger condition \(\Re(a+N) > 1\).

If \(F = \sum_k f_k x^k \in \mathbb{C}[[x]]\), define \(|F| = \sum_k |f_k| x^k\) and \(|F| \le |G|\) if \(\forall_k : |f_k| \le |g_k|\). It is easy to check that \(|F + G| \le |F| + |G|\) and \(|FG| \le |F||G|\). With this notation,

\[|R(s+x)| = \left|\int_N^{\infty} \frac{\tilde B_{2M}(t)}{(2M)!} \frac{(s+x)_{2M}}{(a+t)^{s+x+2M}} dt\right|\]\[\le \int_N^{\infty} \left| \frac{\tilde B_{2M}(t)}{(2M)!} \frac{(s+x)_{2M}}{(a+t)^{s+x+2M}} \right| dt\]\[\le \frac{4 \left| (s+x)_{2M} \right|}{(2 \pi)^{2M}} \int_N^{\infty} \left| \frac{1}{(a+t)^{s+x+2M}} \right| dt\]

where the fact that \(|\tilde B_{2M}(x)| < 4 (2M)! / (2\pi)^{2M}\) has been invoked. Thus it remains to bound the coefficients \(R_k\) satisfying

\[\int_N^{\infty} \left| \frac{1}{(a+t)^{s+x+2M}} \right| dt = \sum_k R_k x^k, \quad R_k = \int_N^{\infty} \frac{1}{k!} \left| \frac{\log(a+t)^k}{(a+t)^{s+2M}} \right| dt.\]

Writing \(a = \alpha + \beta i\), where by assumption \(\alpha + t \ge \alpha + N \ge 1\), we have

\[|\log(\alpha + \beta i + t)| = \left|\log(\alpha + t) + \log\left( 1 + \frac{\beta i}{\alpha + t}\right) \right| \le \log(\alpha + t) + \left|\log\left(1 + \frac{\beta i}{\alpha+t}\right)\right|\]\[= \log(\alpha+t) + \left|\frac{1}{2}\log\left(1+\frac{\beta^2}{(\alpha+t)^2}\right) + i\tan^{-1}\left(\frac{\beta}{\alpha + t}\right)\right| \le \log(\alpha + t) + C\]

where

\[C = \frac{1}{2}\log\left(1+\frac{\beta^2}{(\alpha+N)^2}\right) + \tan^{-1}\left(\frac{|\beta|}{\alpha+N}\right) \le \frac{\beta^2}{2 (\alpha+N)^2} + \frac{|\beta|}{(\alpha+N)}.\]

Also writing \(s = \sigma + \tau i\), where by assumption \(\sigma + 2M > 1\), we have

\[\frac{1}{|(\alpha+\beta i+t)^{\sigma+\tau i + 2M}|} = \frac{e^{\tau \operatorname{arg}(\alpha+\beta i+t)}}{|\alpha+\beta i+t|^{\sigma+2M}} \le \frac{K}{(\alpha + t)^{\sigma+2M}}\]

where \(K = \exp(\max(0, \tau \tan^{-1}(\beta / (\alpha + N))))\). Finally,

\[R_k \le \frac{K}{k!} \, I_k(N+\alpha, \sigma + 2M, C)\]

with the \(K\) and \(C\) defined above, where \(I_k(A,B,C)\) denotes the sequence of integrals

\[I_k(A,B,C) \equiv \int_A^{\infty} t^{-B} (C + \log t)^k dt\]

which can be evaluated as

\[I_k(A,B,C) = \frac{L_k}{(B-1)^{k+1} A^{B-1}}\]

where \(L_0 = 1\), \(L_k = k L_{k-1} + D^k\) and \(D = (B-1) (C + \log A)\).

void zeta_series_em_choose_param(fmpr_t bound, ulong * N, ulong * M, const fmpcb_t s, const fmpcb_t a, long d, long target, long prec)

Chooses N and M using a default algorithm.

Power sums

void zeta_log_ui_from_prev(fmprb_t log_k1, ulong k1, fmprb_t log_k0, ulong k0, long prec)

Computes \(\log(k_1)\), given \(\log(k_0)\) where \(k_0 < k_1\). At high precision, this function uses the formula \(\log(k_1) = \log(k_0) + 2 \operatorname{atanh}((k_1-k_0)/(k_1+k_0))\), evaluating the inverse hyperbolic tangent using binary splitting (for best efficiency, \(k_0\) should be large and \(k_1 - k_0\) should be small). Otherwise, it ignores \(\log(k_0)\) and evaluates the logarithm the usual way.

void zeta_powsum_series_naive(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, long n, long len, long prec)

void zeta_powsum_series_naive_threaded(fmpcb_ptr z, const fmpcb_t s, const fmpcb_t a, long n, long len, long prec)

Computes

\[z = S(s,a,n) = \sum_{k=0}^{n-1} \frac{1}{(k+a)^{s+t}}\]

as a power series in \(t\) truncated to length len. This function evaluates the sum naively term by term. The threaded version splits the computation over the number of threads returned by flint_get_num_threads().

void zeta_powsum_one_series_sieved(fmpcb_ptr z, const fmpcb_t s, long n, long len, long prec)

Computes

\[z = S(s,1,n) \sum_{k=1}^n \frac{1}{k^{s+t}}\]

as a power series in \(t\) truncated to length len. This function stores a table of powers that have already been calculated, computing \((ij)^r\) as \(i^r j^r\) whenever \(k = ij\) is composite. As a further optimization, it groups all even \(k\) and evaluates the sum as a polynomial in \(2^{-(s+t)}\). This scheme requires about \(n / \log n\) powers, \(n / 2\) multiplications, and temporary storage of \(n / 6\) power series. Due to the extra power series multiplications, it is only faster than the naive algorithm when len is small.