USTC, School of Mathematics, Young

  黎曼猜想是德国数学家伯恩哈德·黎曼（Riemann）于1859年提出。它是数学中一个重要而又著名的未 解决的问题，有“猜想界皇冠”之称，多年来它吸引了许多出色的数学家为之绞尽脑汁。黎曼猜想的具体表述形式如下：

Riemann猜想(Riemann Hypothesis, RH)：函数 $$\zeta(s)=\sum_{n=1}^\infty \frac{1}{n^s} \tag{1}$$ 在解析延拓到 $\mathbb{C}\setminus{1}$ 后，其全部非平凡零点都位于直线 $\mathfrak{R}s=\frac{1}{2}$上。

其中 $\zeta(s)$ 称为Riemann $\zeta$函数，它的平凡零点是全体负偶数（ $-2, -4, \dots$ ）除此之外的零点称为非平凡零点。直线 $\mathfrak{R}s=\frac{1}{2}$ 也被称为特征线(critical line)；区域 ${z: 0;\mathfrak{R}z;1 }$ 被称为特征带(critical strip)。

  黎曼猜想之所以被认为是当代数学中一个重要的问题，主要是因为很多深入和重要的数学和物理结果都能在它成立的大前提下得到证明。大部分数学家也相信黎曼猜想的正确性。克雷数学研究所设立了1,000,000美元的奖金给予第一个得出正确证明的人。1

  尽管自从黎曼猜想至今仍未被证明，但是许多数值模拟验证了该猜想在大规模计算下的准确性。其中较为著名的有：2

  本项目基于python语言，验证了 $\zeta$ 函数在 $\mathfrak{I}s\leq 1000$ 内的非平凡零点(共计649个)均位于特征线上，并计算了前500个非平凡零点的近似值（精确到 $10^{-4}$ ）。其中用到的数学原理如下：

  由于 $\zeta$ 函数的原始定义域中并不包含特征带，我们不能直接使用定义中的级数进行计算。由交错级数的收敛理论，我们有如下的命题：3

Proposition: 令 $$e_k=\sum_{j=k}^n\binom{n}{j}$$ 则有： $$\zeta(s)=\frac{1}{1-2^{1-s}} \bigg ( \sum_{k=1}^n\frac{(-1)^{k-1}}{k^s} +\frac{1}{2^n}\sum_{k=n+1}^{2n}\frac{(-1)^{k-1}e_{k-n}}{k^s}\bigg )+\gamma_n(s) \tag{2}$$ 其中对于复变量 $s=\sigma+it$， $\sigma;0$ ，有如下的误差项估计： $$\left\lvert \gamma_n(s) \right\rvert \leq \frac{1}{8^n}\frac{(1+\left\lvert t/\sigma \right\rvert)e^{\left\lvert t \right\rvert \pi /2}}{\left\lvert 1-2^{1-s} \right\rvert} \tag{3}$$

  由上述命题，即可通过交错级数来逼近 $\zeta$ 函数在特征线上的值，并且误差项 $\gamma_n(s)$ 能够得到很好的控制。理论上，只要 $n$ 取得足够大，误差项可以被控制到任意指定的小量。

  由(3)可以得知，在给定误差要求的前提下，利用公式(2)需要累加的项数为 $O(\left\lvert t \right\rvert)$ 。

  除了1.2.1中给出的交错级数求和以外，还可以用Riemann-Siegel公式进行计算。

Riemann-Siegel公式：设 $x, y\in \mathbb{R}^+$ ，满足 $2\pi xy=\left\lvert t \right\rvert$ . 则对于 $s=\sigma+it, 0 \leq \sigma \leq 1$ ，有： $$\zeta(s)=\sum_{n\leq x}\frac{1}{n^s}+\chi(s)\sum_{n\leq y}\frac{1}{n^{1-s}}+O(x^{-\sigma})+O({\left\lvert t \right\rvert}^{1/2-\sigma}y^{\sigma-1}) \tag{4}$$ 其中 $$\chi(s)=2^s\pi^{s-1}\sin(\frac{\pi s}{2})\Gamma(1-s)$$ .

  特别地，在(4)中取 $x=y=\sqrt{\left\lvert t \right\rvert /(2\pi)}$ ，则有 $$\zeta(s)=\sum_{n\leq x}\frac{1}{n^s}+\chi(s)\sum_{n\leq y}\frac{1}{n^{1-s}}+E_m(s), m=\bigg\lfloor\Big(\frac{\left\lvert t \right\rvert}{2\pi}\Big)^{1/2}\bigg\rfloor \tag{5}$$ 其中， $$E_m(s)=O({\left\lvert t \right\rvert}^{-\sigma/2})$$ .

  对于 $s=\sigma+it$ 的情形，我们定义 $$\theta(t)=\arg\Big(\pi^{-it/2}\Gamma(\frac{1}{4}+i\frac{t}{2}\Big )) \tag{6}$$ 有 $$e^{i\theta(t)}=\chi(\frac{1}{2}+it)^{-1/2}$$

与此同时，由Taylor公式，可以得到 $\theta(t)$ 的渐进公式： $$\theta(t)=\frac{t}{2}\log{\frac{t}{2\pi}}-\frac{t}{2}-\frac{\pi}{8}+\frac{1}{48t}+\frac{7}{5760t^3}+\dots \tag{7}$$

  我们再定义Riemann-Siegel Z-函数(也成为Hardy Z-函数)： $$Z(t)=e^{i\theta(t)}\zeta(\frac{1}{2}+it) \tag{8}$$ 可以证明， $Z(t)$ 是实数到实数的映射，从而是实函数。并且 $Z(t)$ 和 $\zeta(\frac{1}{2}+it)$ 有着相同的零点。

  结合(5)(6)(7)，我们可以得到 $$Z(t)=2\sum_{n=1}^m\frac{\cos(\theta(t)-t\log{n})}{\sqrt{n}}+(-1)^{m+1}\tau^{-1/2}\sum_{j=0}^M(-1)^j\tau^{-j}\Phi_j(z)+R_M(t) \tag{9}$$ 其中 $$R_M=O(t^{-(2M+3)/4}),\ \tau=\sqrt{\frac{t}{2\pi}}, \ m=\lfloor \tau\rfloor,\ z=2(\tau-m)-1 \tag{10}$$

余项中的函数 $$\Phi_0(z)=\frac{\cos(\frac{1}{2}\pi z^2+\frac{3}{8}\pi)}{\cos(\pi z)} \tag{11}$$ $$\Phi_1(z)=\frac{1}{12\pi^2}\Phi_0^{(3)}(z) \tag{12}$$ $$\dots\ \dots$$

对于(8)中的误差项 $R_M(t)$ ，当 $t\geq 200$ 时，有如下的估计：3 $$\left\lvert R_0(t) \right\rvert \leq 0.127t^{-3/4},\ \left\lvert R_1(t) \right\rvert \leq 0.053t^{-5/4}$$ 由此可加，在 $t$ 不太小的情况下，只需要取 $M=1$ 即可有效地控制误差。

  对于有界特征带 $C_T$ 表示实部在0到1，虚部在0到T中的复数， $\zeta$ 函数在其中的零点个数有如下的刻画：4

Backlund零点计数定理：用 $N(T)$ 表示 $\zeta$ 函数在 $C_T$ 上的零点个数，假设 $\zeta(s)$ 在线段 ${ \sigma+iT : 0 \leq \sigma \leq 1 }$ 上没有零点，则有 $N(T)=$ 和 $\frac{\theta(T)}{\pi}+1 \tag{13}$最接近的自然数。

从而，只需要计算出(6)中定义的 $\theta(T)$ ，即可得出在 $C_T$ 上的零点总数。

  对于实值函数，我们有如下的定理：

介值定理：设 $f$ 是 $\mathbb{R}$ 上的实值连续函数，如果存在 $a ; b, s.t. f(a)\cdot f(b) ; 0$, 则 $f$ 在区间 $(a, b)$ 上必有零点。

  有了上述结论，我们可以采用如下思路来数值验证Riemann猜想：

  需要补充说明的是，在步骤2中，我们注意到 $\zeta(\frac{1}{2}+it)$ 和 $Z(t)$ 有相同的零点，从而只需要计算 $Z(t)$ 的零点个数；而 $Z(t)$ 是连续的实函数，可以用介值定理来“探测”其零点。

  本项目的程序主要由3个模块组成：

定义的常量：

实现的函数（只列出主要的函数）：

定义的常量：

实现的函数：

实现的函数：

分析用交错级数方法和Riemann-Siegel方法逼近 $\zeta(s)$ 或 $Z(t)$ 的时间复杂度，可以得知：在给定误差上界的情况下，式(2)中累加的项数为 $O(\left\lvert t \right\rvert)$ ，而式(9)中累加的项数为 $O(\sqrt{\left\lvert t \right\rvert})$ 。从而在 $T$ 较大时，使用Riemann-Siegel方法的计算效率会明显高于交错级数方法。

另一方面，使用式(2)计算时，误差项 $\gamma_n(s)$ 可以任意小，而式(9)中累加次数 $m$ 是由 $t$ 确定的，只有在 $\left\lvert t \right\rvert$ 较大时，才能通过误差项公式 $\left\lvert R_0(t) \right\rvert \leq 0.127t^{-3/4},\ \left\lvert R_1(t) \right\rvert \leq 0.053t^{-5/4}$ 有效地控制误差，因此从计算精度考虑，交错级数方法的计算误差能被更好地控制。

综合来看，当 $t$ 的模长较小时，我们使用交错级数方法，此时累加次数虽然是 $O(\left\lvert t \right\rvert)$ ，但是由于 $\left\lvert t \right\rvert$ 不太大，所需要地绝对时间仍然可以接受； $\left\lvert t \right\rvert$ 较大时，我们使用Riemann-Siegel公式，此时由于 $\left\lvert t \right\rvert$ 较大，误差项能被控制在可以接受的范围内，并且此时的计算效率也相对较高。

具体而言，我们选定两种计算方法的分界点为CHANGE_METHOD=200，此时用Riemann-Siegel方法的计算误差项满足 $$\left\lvert R_1(t) \right\rvert \leq 0.053t^{-5/4}\leq 7\times 10^{-5}$$ 足够我们的要求。

其次，本项目也从不同角度对 $\zeta$ 函数进行了可视化，包括 $\zeta$ 函数的实部、虚部在特征线上的动态变化图像，zeta函数的模长和Riemann-Siegel- $Z$ 函数的对比图像， $\zeta$ 函数在特征线上的实部、虚部分别的图像（具体见第4.3）。

在程序的最后输入语句：

check_RH(1000, 0.1)

运行即可。

注意：参数delta不能设置得过大，否则在采样时可能会漏掉零点；也不能设置得过小，否则运行时间会很长。

此程序会在当前目录生成一个名为"zeros_records.txt"的文本文件，其中每行存储了零点的编号、存在区间，例如

Zero No.1: (14.1, 14.2) Zero No.2: (21.0, 21.1) Zero No.3: (25.0, 25.1)

在此文本文件的最后存储了找到的零点个数和通过Backlund定理算出的零点个数，例如

Expecting 649.0 zeros. Find 649 zeros.

比较上面的Expecting 649.0 zeros和Find 649 zeros，可以得知，在 $C_{1000}$ 上Riemann猜想成功验证。

以及程序运行所用时间：

Total time cost: 280.6123475 seconds. Average time cost: 0.4323765 seconds per zero.

由于此程序运行时间较长，我们加入了进度条，从而可以直观地观测到程序的运行进度。

在程序的最后输入语句

compute_from_file(FILE_NAME, 500)

运行即可。

此程序会从文件FILE_NAME中读取零点的存在区间，并计算零点的精确值（精确到 $10^{-5}$ ），并在当前目录生成一个名为"compute_zeros.txt"的文本文件，其中存储了计算得出的前500个零点的数值，例如

Zero No 1: 14.13473 Zero No 2: 21.02204 Zero No 3: 25.01086

在文件的最后存储了计算的精度和运行时间，例如

With accuracy: 1e-05 Total time cost: 962.0139017 seconds. Average time cost: 1.9240278034 seconds per zero.

可以通过修改全局变量ACCURACY的值来调节计算精度（默认为1e-5）。

为了检验零点计算的准确性，我们可以将文件compute_zeros.txt中的计算结果和LMFDB(The L-functions and modular forms database)中的 zeta函数零点表 对比，可见计算结果的小数点后前4位均为准确值。

在程序的最后输入语句

draw_zeta_critical_line(1, 50)

并运行，可以得到在参数1