Fast methods to compute the Riemann zeta function (Q640804)

This paper presents an algorithm which will compute \(\zeta(1/2+it)\) for \(t\geq 2\) in time \(O(t^{4/13+o(1)})\). This should be compared with the classical approach, via the Riemann--Siegel formula, taking time \(O(t^{1/2+o(1)})\), with the algorithm due to \textit{A. Schönhage} [Jahresber. Dtsch. Math.-Ver. 92, No. 1, 1--20 (1990; Zbl 0797.68090)], whose complexity is \(O(t^{3/8+o(1)})\), and with a method sketched informally by the reviewer, which uses \(O(t^{1/3+o(1)})\) operations. A precise statement of the result is that there is an absolute constant \(\kappa\) such that for any given \(\lambda>0\) there is a constant \(A(\lambda)\) with the following property. For any \(t\geq 2\) one can compute \(\zeta(1/2+it)\) to within \(\pm t^{\lambda}\), using at most \(A(\log t)^{\kappa}t^{3/14}\) operations on numbers of at most \(A(\log t)^2\) bits, and at most \(A(\log t)^{\kappa}t^{3/14}\) bits of storage. There is a sense in which the reviewer's algorithm, with complexity \(O(t^{1/3+o(1)})\), corresponds to the estimate \(\zeta(1/2+it)\ll t^{1/6}\log t\), proved by the method of Hardy and Littlewood. However it is noteworthy that the exponent \(3/14\) in the present paper is smaller than any current estimate for \(\mu(1/2)\). One explanation for this lies on the author's use of a precomputation to tabulate \(O((\log t)^{\kappa}t^{3/14})\) values of exponential sums -- a process which has no equivalent in estimating the size of \(\zeta(s)\). The algorithm builds on ideas in the author's paper [Ann. Math. (2) 174, No. 2, 859--889 (2011; Zbl 1243.11117)], in which an iterative algorithm for computing quadratic exponential sums is presented. In the present paper cubic sums \[ S:=\sum_{k\leq K}\exp\{2\pi i(ak+bk^2+ck^3)\} \] with small values of \(c\) are used. These arise essentially by splitting a sum \(\sum n^{-it}\) into sections of length \(K\) and approximating \(t\log(N+k)\) by the cubic terms in its power series. The values of \(a\) and \(b\) in \(S\) are normalized so that \(0\leq a\leq 1\) and \(0\leq b\leq 1/4\), and then the van der Corput B-process is applied. Providing \(c\) is sufficiently small in terms of \(K\) the \(B\)-process again produces a perturbation of a quadratic sum. The resulting error terms must be analyzed carefully, but the outcome is that the sum \(S\) is replaced by a similar one in which \(K\) is reduced by a factor of about \(1/2\), and \(c\) is increased. This process, of normalizing the values of \(a\) and \(b\) and then using the \(B\)-process, is repeated until the value of \(c\) becomes too large for the method to work. By this stage the value of \(K\) has been reduced sufficiently that \(S\) can be calculated by reference to the precomputed table.