New, high-performance algorithms for the computation of the Hardy Function

The Riemann zeta-function ζ (𝑠) for 𝑠 ∈ ℂ, is a complex function fundamental to prime number theory. Computation of ζ (𝑠) along the critical line 𝑠 = 1/2 + I𝑡, 𝑡 ∈ ℝ, is important for the distribution of primes, since all the zeta-zeros are predicted to lie here (the Riemann Hypothesis). Hardy’s Z-function 𝑍 (𝑡) = 𝑒I𝜃 𝑡 ζ (1/2 + I𝑡), 𝜃 (𝑡) ∈ ℝ, defines the amplitude of the zeta-function along the critical line. Such computations also provide a means of Riemann Hypothesis verification, within closed intervals, and estimates on the bounds of ζ (1/2 + I𝑡), central to the Lindelöf Hypothesis. Until recently, the most efficient means of computing 𝑍 (𝑡) employed the Riemann-Siegel (RS) formula, an 𝑂 (√ 𝑡) operational method. In 2011, this was superseded by an algorithm devised by G. A. Hiary, requiring just 𝑂 (𝑡1/3𝑙𝑜𝑔 𝑡 𝜅) operations. The methodology involved the sub-division of the RS formula into sequences of quadratic Gauss sums of length 𝑁. Such sums are amenable to rapid computation, in 𝑂 (𝑙𝑜𝑔 (𝑁) ) operations, using a recursive scheme. Central to this work is a new asymptotic formula for 𝑍 (𝑡), with some interesting analytical properties. Computationally, this new formulation allows one to improve on the quadratic sum expansion and express 𝑍 (𝑡) in terms of sub-sequences of 𝑚𝑡ℎ-order Gauss sums. In the cubic 𝑚 = 3 case, these sums can be computed rapidly, utilising a similar scheme to that used for the quadratics. The result is a more efficient algorithm, requiring only 𝑂 ( (𝑡 /𝜀𝑡) 1/4 𝑙𝑜𝑔 (𝑡) 3) operations. Sample computations, using the open-source code developed for this eCSE Project and available on GitHub repository, support these findings.