🔢 mathematics

Simple and accurate approximations to the Riemann zeta function

This paper presents simple and accurate approximations for the Riemann zeta function and its derivative within the critical strip by combining the main sum of the Riemann-Siegel formula with an elementary-function-based remainder term derived via Gaussian quadrature, supported by extensive numerical evidence of their high precision.

Original authors: Alexey Kuznetsov

Published 2026-05-22

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Alexey Kuznetsov

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the Riemann zeta function as a massive, incredibly complex machine that mathematicians have been trying to understand for centuries. It's like a giant, invisible engine that powers some of the deepest mysteries in number theory. To "run" this machine and get a specific number out of it (a value for $\zeta(s)$ ), you usually have to feed it a lot of data.

For a long time, the best way to do this was like trying to count every single grain of sand on a beach to estimate the total volume. If the beach is small, it's easy. But if the beach is huge (which happens when the numbers get very large), counting every grain takes forever. This is the problem with the old methods: they require summing millions of terms just to get a single answer, making them slow and impractical for big numbers.

The New "Shortcut" Engine

In this paper, Alexey Kuznetsov introduces a new, clever shortcut. Instead of counting every single grain of sand, he builds a highly accurate model that predicts the total volume based on a few key measurements.

Think of the old method (the Riemann-Siegel formula) as a recipe that requires you to bake a huge cake, but then you have to throw away most of it because it's too messy. The "remainder" of the cake is hard to calculate, so you usually just guess it or use a very complicated, fussy method to get it right.

Kuznetsov's new method is like having a magic measuring cup. He realizes that the "messy remainder" part of the cake follows a very specific, predictable pattern. Instead of baking the whole thing, he pre-computes a set of 21 special numbers (for a specific version of the tool) that act like a secret code. These numbers are like a pre-tuned set of weights on a scale.

How the Magic Works

The Pre-Computed Weights: Before you even start, the author has calculated a list of special complex numbers (think of them as "magic coefficients"). He has done this for different levels of precision (like having a "small," "medium," and "large" version of the tool). You can download these numbers; you don't have to calculate them yourself every time.
The Simple Recipe: To get the answer, you take the main part of the calculation (which is easy) and add a correction term. This correction term doesn't require summing millions of items. Instead, it just uses those pre-computed "magic weights" combined with simple math functions (like exponentials and logarithms).
The Result: It's like swapping a 10-hour manual labor job for a 5-second button press. The paper shows that for very large numbers, this new method is incredibly accurate—often accurate to 100 or even 300 decimal places—without needing the heavy lifting of the old methods.

Why It's Better Than the Old Way

The paper compares this new method to the "gold standard" (the Riemann-Siegel formula).

The Old Way: To get high accuracy, you need to calculate many complicated correction terms. These terms are like trying to solve a puzzle where the pieces are shaped like jagged rocks; they are hard to fit together, and the formulas are messy and prone to errors.
The New Way: You only need a small, fixed list of pre-computed numbers (about 15 to 30 numbers depending on how precise you want to be). The math involved is smooth and simple, like fitting together smooth, round marbles.

The "Sweet Spot" Pattern

The author noticed something fascinating in the errors. When he tested the method, the mistakes weren't random. They danced in a specific pattern, like a heartbeat. The error would be tiny at certain regular intervals and slightly larger in between. This pattern actually helped him figure out how to choose the best "magic weights" in the first place. It's as if the machine "whispered" the secret to him through the pattern of its mistakes.

What This Means for Math

The paper doesn't claim this will cure diseases or predict the weather. Instead, it solves a specific, stubborn problem in pure mathematics: how to calculate this specific function quickly and accurately.

Speed: It allows computers to calculate these values hundreds or thousands of times faster than the built-in tools in software like MATLAB.
Precision: It allows mathematicians to see the function with extreme clarity (hundreds of decimal places), which is crucial for testing deep theories about prime numbers.
Simplicity: It replaces a complex, multi-step process with a simple formula that relies on a pre-made "cheat sheet" of numbers.

In short, Kuznetsov has built a high-precision, easy-to-use calculator for a function that was previously difficult and slow to compute, using a clever trick of pre-computed "magic numbers" to bypass the heavy lifting.

Technical Summary: Simple and Accurate Approximations to the Riemann Zeta Function

Problem Statement

The computation of the Riemann zeta function, $\zeta(s)$ , for large values of the imaginary part, $t = \text{Im}(s)$ , presents significant challenges. While the Euler-Maclaurin summation method offers rigorous error bounds, its computational complexity scales as $O(|s|)$ , rendering it impractical for large $t$ . The standard approach for large $t$ is the Riemann-Siegel formula, which reduces complexity to $O(t^{1/2})$ by utilizing a main sum and an asymptotic expansion for the remainder term. However, higher-order corrections in the Riemann-Siegel formula involve complex derivatives of special functions (specifically the function $F(z)$ ), which are difficult to evaluate numerically due to singularities and complicated formulas. Furthermore, existing high-precision algorithms, such as those by Galway, often require tuning multiple parameters to achieve optimal accuracy.

The paper addresses the need for approximations that are both computationally efficient (requiring few terms) and highly accurate across arbitrary vertical strips, including the critical strip, without the implementation difficulties associated with higher-order Riemann-Siegel corrections.

Methodology

The author develops a family of approximations, denoted $\zeta_p(s)$ , which combine the main sum of the Riemann-Siegel formula with a novel approximation of the remainder term. The core methodology relies on the following steps:

Integral Representation: The remainder term of the exact Riemann-Siegel formula is expressed as an integral $I_M(s)$ involving a complex-valued measure $\eta(dx) = \frac{e^{-\pi x^2}}{\cosh(\pi \theta x)} dx$ , where $\theta = e^{-i\pi/4}$ .
Gaussian Quadrature: Instead of using asymptotic series expansions, the author approximates this integral using a Gaussian quadrature rule. The goal is to find a set of nodes $\{x_{p,j}\}$ and weights $\{\omega_{p,j}\}$ such that the quadrature is exact for a specific set of $4p+2$ test functions.
Coefficient Derivation: The coefficients $\{\omega_{p,j}\}$ ${ω_{p, j}}$ and $\{\lambda_{p,j}\}$ ${λ_{p, j}}$ (where $\lambda_{p,j} = \theta x_{p,j}$ $λ_{p, j} = θ x_{p, j}$ ) are derived by solving a system of linear equations. This system is constructed by requiring the quadrature to be exact for functions of the form $f_k(x) = e^{-\pi x^2 + 2\pi \theta y_{p,k} x}$ $f_{k} (x) = e^{- π x^{2} + 2 π θ y_{p, k} x}$ , where $y_{p,k}$ $y_{p, k}$ are equally spaced points in $[-1, 1]$ $[- 1, 1]$ .
- The problem is transformed into finding the nodes and weights of a Gaussian quadrature for a linear functional $L$ defined by the moments of the Mordell integral $H(y)$ .
- These nodes and weights are computed using orthogonal polynomials generated via a three-term recurrence relation.
Approximation Formula: The final approximation is defined as:
$\zeta_p(s) = \sum_{n=1}^{N_t} n^{-s} + \chi(s) \sum_{n=1}^{N_t} n^{s-1} - \frac{(-1)^{N_t}}{2} \left[ I_{N_t+1/2, p}(s) + \chi(s) I_{N_t+1/2, p}(1-s) \right]$
where $I_{M,p}(s)$ is the quadrature approximation of the integral term, involving only elementary functions (exponentials, logarithms) and precomputed complex coefficients.

Key Contributions

Simplified Implementation: The proposed approximations depend on a single integer parameter $p$ and a set of precomputed coefficients $\{\omega_{p,j}\}$ and $\{\lambda_{p,j}\}$ . Unlike the Riemann-Siegel formula, which requires evaluating complex derivatives of special functions for higher-order terms, $\zeta_p(s)$ requires only elementary functions and the Gamma function.
Precomputed Coefficients: The author provides precomputed coefficients for $1 \le p \le 30$ and selected values up to $p=150$ . These coefficients are shown to be of moderate magnitude and exhibit rapid decay, facilitating stable numerical evaluation.
Derivative Approximation: The paper extends the methodology to approximate the derivative $\zeta'(s)$ by differentiating the approximation formula with respect to $s$ (evaluated at fixed $N_t$ ), avoiding the discontinuities associated with differentiating the floor function in $N_t$ .
Hybrid Algorithm: A strategy is proposed for computing $\zeta(s)$ across the entire complex plane by switching between direct summation, Euler-Maclaurin summation, and the $\zeta_p(s)$ approximation based on the region of $s$ .

Results

The paper presents extensive numerical evidence regarding the accuracy of $\zeta_p(s)$ :

Accuracy in the Critical Strip: Numerical experiments in the critical strip ( $0 \le \sigma \le 1$ $0 \leq σ \leq 1$ ) demonstrate that the error $|\zeta_p(s) - \zeta(s)|$ $∣ ζ_{p} (s) - ζ (s) ∣$ decreases significantly as $p$ $p$ increases.
- For $p=10$ , errors are $< 10^{-15}$ for $t > 250$ .
- For $p=50$ , errors are $< 10^{-100}$ for $t > 4000$ .
- For $p=150$ , errors are $< 10^{-300}$ for $t > 6900$ .
Error Patterns: The error exhibits a distinct pattern where it vanishes at the endpoints of intervals $[t_n, t_{n+1}]$ (where $t_n = 2\pi n^2$ ) and at $4p$ equally spaced points within these intervals. This behavior is analytically explained by the fact that the dominant error term is proportional to $H_p(B(t)) - H(B(t))$ , which is zero at the quadrature nodes.
Comparison with Riemann-Siegel: For moderate $t$ (e.g., $10^2 \le t \le 10^4$ ), the approximation $\zeta_7(s)$ is shown to be more accurate than the Riemann-Siegel approximation with $K=5$ ( $\zeta^{RS}_5$ ). While $\zeta^{RS}_5$ has a faster asymptotic decay rate ( $O(t^{-13/4})$ ) compared to $\zeta_7$ ( $O(t^{-1/4})$ ), $\zeta_7$ requires significantly fewer precomputed coefficients (15 vs. ~100–170 for high precision) and is easier to implement.
Precision Limits: In double-precision implementations, rounding errors dominate the approximation error for $t > 200$ . However, in quadruple precision, the approximations maintain high accuracy up to much larger $t$ (e.g., $t > 10^4$ with errors $< 10^{-29}$ ).
Performance: The MATLAB implementation of the hybrid algorithm is reported to be hundreds to thousands of times faster than the built-in MATLAB zeta function.

Significance and Claims

The paper claims that these approximations enable high-precision computation of $\zeta(s)$ and $\zeta'(s)$ in the critical strip and other vertical strips for large $t$ . The significance lies in the trade-off between implementation complexity and accuracy: the method offers a "simple" alternative to the Riemann-Siegel formula that avoids the computational burden of high-order correction terms while maintaining competitive or superior accuracy for a wide range of $t$ .

The author notes that while rigorous error bounds for $\zeta_p(s)$ are not provided in this work, the numerical evidence is extensive. The method is presented as a practical tool for high-precision computations, which are crucial for mathematical discoveries such as the disproof of the Mertens conjecture. The paper concludes by providing open-source implementations in MATLAB, Python, and Fortran, facilitating immediate use by the research community.