Zeta functions, polylogarithms and L-functions

The Hurwitz zeta function

and its analytic continuation to s \ne 1 are the core of a family of functions: the Riemann zeta function \zeta(s) = \zeta(s, 1), the Dirichlet L-functions, which are linear combinations of Hurwitz zeta values, and the polylogarithms, which are linear combinations of two Hurwitz zeta values with complex parameter. FLINT computes \zeta(s, a) and its derivatives with respect to s by Euler–Maclaurin summation with the rigorous error bounds of [Joh2013], and uses specialised methods where they are faster: exact formulas and the Euler product at integers, the Riemann–Siegel formula high on the critical line (Zeros of the Riemann zeta function), and a Taylor expansion in a when many values with the same s are needed.

Overview of methods#

Euler–Maclaurin summation#

Applying the Euler–Maclaurin formula to f(t) = (a + t)^{-s} on [N, \infty) gives the following representation, valid for all s \ne 1 by analytic continuation.

[Joh2013] extends the bound to complex a (with a factor depending on \arg(a + N)) and to the derivatives with respect to s, by treating s + x as a power series: the remainder of the expansion of \zeta(s + x, a) to length d is bounded coefficientwise by an explicit integral of (a+t)^{-\sigma-2M} \log^k(a+t), which is what _acb_poly_zeta_em_bound computes.

Choice of parameters. The terms of the Bernoulli sum first decrease like |s + 2j|^{2j}/(2\pi(a+N))^{2j} and then grow; the bound is small when N and M are both large compared with |s| and the precision. _acb_poly_zeta_em_choose_param takes M = \min(N, p + \min(N/100, 2000)) for a target of p bits and finds the smallest N for which the bound is below 2^{-p} times an estimate of |a^{-s}|, by doubling and bisection. For s of moderate size this gives N, M = O(p). The cost is dominated by the N powers (a+k)^{-s} and the M Bernoulli numbers (cached by bernoulli_cache_compute). The Bernoulli sum is evaluated term by term (_acb_poly_zeta_em_tail_naive), or by binary splitting (_acb_poly_zeta_em_tail_bsplit) for long power series (length at least 5 and at least M/10). For large |t|, the remainder bound requires N \gtrsim |t|/(2\pi), so the cost grows linearly with the height; this is why the Riemann–Siegel formula, with O(\sqrt{t}) terms, takes over there.

Power sums. The main sum \sum_{k<N} (a+k)^{-(s+x)} is computed by _acb_poly_zeta_em_sum with one of several methods:

• for a = 1 and length \le 2: acb_dirichlet_powsum_sieved, which computes k^{-s} only for prime k and obtains the others as products (ij)^{-s} = i^{-s} j^{-s} from a table, groups the even terms as a polynomial in 2^{-s} evaluated by Horner's rule, and so needs about N/\log N exponentials and N/2 multiplications;
• for a = 1 and length \le 4: acb_dirichlet_powsum_smooth, which sieves only by the primes 2, 3 and 5 and uses logarithmic memory;
• for longer series or general a: a product tree (_acb_poly_powsum_series_tree), or a direct sum, in parallel when several threads are available.

The powers k^{-s} = \exp(-s \log k) are computed with the logarithm of a nearby integer when possible (acb_dirichlet_powsum_term), and the precision of each term is lowered according to its magnitude.

Reflection. For \operatorname{Re}(s) < 0 the Riemann zeta function uses the functional equation \zeta(s) = 2 (2\pi)^{s-1} \sin(\pi s/2) \Gamma(1-s) \zeta(1-s), also as power series, since the cost of Euler–Maclaurin summation grows with |s|.

The zeta function at integers#

At even positive integers, \zeta(2n) = (-1)^{n+1} B_{2n} (2\pi)^{2n}/(2 (2n)!) is computed from the Bernoulli number when it is cheap (it is cached), and at negative integers \zeta(-n) = -B_{n+1}/(n+1). For large n the Euler product

truncated at the primes p < P, converges quickly: the relative error is at most about \zeta(n, P) \approx P^{1-n}/(n-1), so P \approx 2^{p/n} primes suffice (arb_zeta_ui_euler_product, which bounds the tail by mag_hurwitz_zeta_uiui). For odd n of moderate size, Borwein's algorithm [Bor2000]

has an error \varepsilon_n(s) of order (3 + \sqrt 8)^{-n} for real s \ge 1/2, so it needs n \approx 0.39 p terms; arb_zeta_ui_vec_borwein evaluates it for several s at once, and arb_zeta_ui_borwein_bsplit evaluates it by binary splitting at very high precision.

Taylor expansion in the parameter#

To evaluate \zeta(s, a) at many a for a fixed s, as needed for the L-functions of all characters of a modulus, acb_dirichlet_hurwitz_precomp_init precomputes Taylor expansions at grid points. With \zeta(s, a) = \sum_{j=0}^{A-1} (a+j)^{-s} + \zeta(s, a + A) and the expansion

it stores K coefficients at each of the N points A + (2i+1)/(2N), 0 \le i < N, which leaves |x| \le 1/(2N) for any a \in (0, 1]. Each evaluation then costs a polynomial evaluation of length K and A powers. The truncation error is bounded as follows.

With a \ge A and |x| \le 1/(2N), C \approx 1/(2AN), so the error is O((2AN)^{-K}); acb_dirichlet_hurwitz_precomp_bound evaluates this bound and acb_dirichlet_hurwitz_precomp_choose_param minimises the total cost NK (precomputation) plus (number of evaluations) \times (K + A) subject to it.

Polylogarithms#

The polylogarithm \operatorname{Li}_s(z) = \sum_{k=1}^\infty z^k k^{-s} for |z| < 1, continued analytically with a branch cut on z \in (1, \infty), is computed by acb_polylog, and its expansion in s by acb_poly_polylog_series.

Small |z|. For |z| < 1/2, the series is summed directly, also for the derivatives with respect to s: \operatorname{Li}_{s+x}(z) = \sum_{d} x^d \frac{(-1)^d}{d!} \sum_{k \ge 1} T_d(k) with T_d(k) = z^k \log^d(k) k^{-s}. The tail after N terms is bounded by mag_polylog_tail: if U(k) is nonincreasing and T_d(k+1)/T_d(k) \le U(k), then \sum_{k \ge N} T_d(k) \le T_d(N) / (1 - U(N)) when U(N) < 1. One can take U(k) = |z| \, B(k, \max(0, -\sigma)) \, B(k\log k, d) with B(m, n) = (1 + 1/m)^n, because (k/(k+1))^{\sigma} \le (1 + 1/k)^{\max(0,-\sigma)} and (\log(k+1)/\log k)^d \le (1 + 1/(k \log k))^d.

General z. Otherwise, with s = 1 - v,

with the principal branch of \log(-z), which gives the conventional branch cut. The Hurwitz zeta functions are evaluated by Euler–Maclaurin summation with complex parameter. To expand in s, v is replaced by v + x, the right side becomes \Gamma(v+x) [E_1 Z_1 + E_2 Z_2] with E_1 = (i/(2\pi))^{v+x} and E_2 = (1/(2\pi i))^{v+x}, and x is mapped to -x at the end. At v = 1 the Hurwitz zeta functions have poles; with the deflated \tilde Z_j = Z_j - 1/x, E_1 Z_1 + E_2 Z_2 = (E_1 + E_2)/x + E_1 \tilde Z_1 + E_2 \tilde Z_2, where (E_1+E_2)/x is a power series since E_1 + E_2 vanishes at x = 0. At v = 0, -1, -2, \ldots, where \Gamma has a pole, the product x\Gamma(v+x) \cdot [E_1 Z_1 + E_2 Z_2]/x is formed instead. These cancellations require s to be an exact integer; a small ball around 0, 1, 2, \ldots gives an indeterminate result. The method is inefficient for |z| \to 0 (where the direct series is used) and for |z| \to \infty, where an asymptotic expansion would be better. For s = 2, acb_hypgeom_dilog is used instead (The dilogarithm).

Dirichlet L-functions#

For a character \chi modulo q,

and acb_dirichlet_l uses this decomposition (acb_dirichlet_l_hurwitz) in general, computing the \phi(q) Hurwitz zeta values with the Taylor expansion of Taylor expansion in the parameter when that is cheaper. At s = 1, the deflated values \zeta(s, a) - 1/(s-1) are used, since the poles cancel for non-principal \chi. acb_dirichlet_l_vec_hurwitz computes L(s, \chi) for all \chi modulo q at once: the map from the vector (\zeta(s, k/q))_k to (L(s, \chi))_\chi is a discrete Fourier transform on the group of characters (acb_dirichlet_dft_index), so the cost is \phi(q) Hurwitz zeta values plus O(\phi(q) \log \phi(q)) operations.

For large \operatorname{Re}(s), from about 8 + p/(2\log p), the Euler product 1/L(s,\chi) = \prod_p (1 - \chi(p) p^{-s}) is used; the truncation at p < P leaves an error at most \zeta(\operatorname{Re}(s), P) in 1/L, because the truncated product contains all terms \mu(k)\chi(k)k^{-s} with P-smooth k.

For primitive \chi of parity \delta, the completed function \Lambda(s, \chi) = (q/\pi)^{(s+\delta)/2} \Gamma((s+\delta)/2) L(s,\chi) satisfies \Lambda(s, \chi) = \epsilon(\chi) \Lambda(1 - s, \overline\chi), which moves evaluations from \operatorname{Re}(s) < -1/2 to the right half-plane. The root number is \epsilon(\chi) = G(\chi)/(i^\delta \sqrt{q}), where the Gauss sum G(\chi) = \sum_x \chi(x) e^{2\pi i x/q} factors over the prime powers dividing q by the Chinese remainder theorem. It can also be obtained from the theta series \Theta(\chi, t) = \sum_{n \in \mathbb{Z}} \chi(n) n^\delta e^{-\pi n^2 t^2/q}, which satisfies t^{1/2+\delta}\Theta(\chi, t) = \epsilon(\chi)\Theta(\overline{\chi}, 1/t), so that \epsilon(\chi) = \Theta(\chi, 1)/\overline{\Theta(\chi, 1)} whenever \Theta(\chi, 1) \ne 0; this needs O(\sqrt{pq}) terms.

The Mellin transform of the theta series gives the approximate functional equation, a sum of incomplete gamma functions,

whose terms decrease like e^{-\pi n^2/q}, so O(\sqrt{pq}) terms suffice at precision p [Coh2007b]. acb_dirichlet_l_fmpq_afe uses it for rational s, for which the upper incomplete gamma function with rational parameter can be computed efficiently (by its asymptotic series when possible), with each term computed at the precision its magnitude requires; it is faster than Euler–Maclaurin summation at very high precision or for large q.

Other functions#

Lerch transcendent. \Phi(z, s, a) = \sum_{k \ge 0} z^k (k+a)^{-s} reduces to the Hurwitz zeta function for z = 1, to the Dirichlet eta function, digamma or Hurwitz zeta functions for z = -1, and to polylogarithms for a = 1. Otherwise acb_dirichlet_lerch_phi sums the series directly when |z| \le 0.875, and else evaluates the Hankel contour integral \Phi(z,s,a) = -\frac{\Gamma(1-s)}{2 \pi i} \int_C \frac{(-t)^{s-1} e^{-a t}}{1 - z e^{-t}} dt with acb_calc_integrate, deforming the contour around the poles t = \log z + 2\pi i k.

Stieltjes constants. The generalised Stieltjes constants \gamma_n(a) are the coefficients of the Laurent expansion of \zeta(s, a) at s = 1. For moderate n, they are read off from the Euler–Maclaurin expansion of \zeta(s + x, a) at s = 1 with deflation; for large n, acb_dirichlet_stieltjes integrates a representation along a contour that avoids the cancellation of the power series method [JB2018], which makes very large n, such as 10^{100}, feasible.

Use in FLINT#

• acb_dirichlet_zeta, acb_dirichlet_zeta_jet, acb_dirichlet_hurwitz, acb_dirichlet_eta, acb_dirichlet_xi (acb_dirichlet); acb_zeta and acb_hurwitz_zeta are wrappers.
• _acb_poly_zeta_cpx_series, _acb_poly_zeta_em_choose_param, _acb_poly_zeta_em_bound, _acb_poly_zeta_em_sum, _acb_poly_zeta_em_tail_bsplit, acb_poly_zeta_series (acb_poly): Euler–Maclaurin summation.
• acb_dirichlet_powsum_sieved, acb_dirichlet_powsum_smooth, _acb_poly_powsum_series_tree: power sums.
• arb_zeta_ui, arb_zeta_ui_euler_product, arb_zeta_ui_vec_borwein, arb_zeta_ui_borwein_bsplit, arb_zeta_ui_bernoulli (arb).
• acb_dirichlet_hurwitz_precomp_init, acb_dirichlet_hurwitz_precomp_bound, acb_dirichlet_hurwitz_precomp_eval.
• acb_polylog, acb_poly_polylog_series, mag_polylog_tail.
• acb_dirichlet_l, acb_dirichlet_l_hurwitz, acb_dirichlet_l_euler_product, acb_dirichlet_l_vec_hurwitz, acb_dirichlet_l_fmpq_afe, acb_dirichlet_root_number, acb_dirichlet_gauss_sum, acb_dirichlet_theta_arb.
• acb_dirichlet_lerch_phi, acb_dirichlet_stieltjes.

The tutorial The Riemann zeta function, its zeros and L-functions evaluates the zeta function, its zeros and L-functions in complete programs.