リーマン予想をPythonで体験してみよう！前編

リーマン予想ってなんなん？って人ーっ！

　「リーマンゼータ関数のすべての非自明な零点の実部は 1 / 2である。」

　これがリーマン予想です。この記事では

• リーマン予想をプログラムを通して理解する
• 何をもって「解いた！」と言えるか理解する

が目的です。今回はプログラムにPythonを利用します。Python の哲学である、batteries included（バッテリーはすでにあるぜ） に従って、基本は標準ライブラリだけで記述していきます。

　ただ、さすがにこれしんどいな...って部分は外部ライブラリに頼ります（グラフとか）。

　Python で数学的計算をするというのに慣れていない人にもオススメかも！

　この記事では純粋な数値計算はコマンドラインで実行します。

>>> func = lambda x: x**2 - 2
>>> func(4) # func(4) は 4 * 4 - 2 = 14

みたいな感じ。
グラフを作成するプログラムなど、少し規模の大きいものは折り畳みしてあるので開いて見てください。

リーマンの論文【参考】

　リーマンの論文は非常にわかりやすく翻訳してくれた下記のサイトを見てください。

　実はリーマン予想はこの論文に現れる複数の予想の一つです。論文自体は「与えられた数より小さな素数の個数について」というタイトルにある通り、素数の個数を求めてみようって内容になっています。

　いずれはこの論文の内容を一から追う記事とか作りたいな ^^

リーマンゼータ関数

　リーマンゼータ関数とは、次のような関数のことです。

\zeta (s) = \sum_{}^{} \frac{1}{n^s} = \frac{1}{1^s} + \frac{1}{2^s} + \frac{1}{3^s} + \dots

　リーマン予想はこの関数の非自明な零点についての予想です。はてはて？？ってなりますね。

　ではこれから「リーマンゼータ関数の非自明な零点」を Python で求める旅に出かけましょう

関数の零点

　まずは関数の零点から。関数とは「何かしらの入力に対して、何かしらの出力を返すもの」です。
　この時、出力を 0 にする入力 を零点と言います。

　いくつか例を見ておきます。

• 一次関数 f(x) = 2x + 1 の零点

>>> func = lambda x: 2 * x + 1 # 一次関数 2x + 1 をラムダ式で宣言
>>> func(4) # func(4) は 9 なので零点じゃない
9
>>> func(-1/2) # func(-1/2) は 0 と出力されるので、1/2 は2x+1の零点！
0.0

• 二次関数 f(x) = x^2 - 2 の零点

>>> func = lambda x: x**2 - 2 # 二次関数 x^2 - 2 をラムダ式で宣言
>>> func(4) # func(4) は 4 * 4 - 2 = 14 で零点じゃない
14
>>> func(2**0.5) # func(2**0.5) は（ほぼ）0と出力されるので√2はx^2-2の零点！
4.440892098500626e-16
>>> func(-2**0.5) # func(-2**0.5) は（ほぼ）0と出力されるので-√2はx^2-2の零点！
4.440892098500626e-16

　二次関数の方で出力が0になっていないですが、計算機イプシロンの関係上これは仕方ない

　零点は簡単に見つかるものからなかなか見つからないものまでたくさんあります。リーマン予想はゼータ関数の零点を見つける問題らしいので、いくつか値を入れて 0 になるか試してみましょう。

>>> def zetaf(s, max): # ゼータ関数の定義。n が 1 から max まで計算します（ほんとは 1 から無限まで）
...     return sum([1.0 / (i**s) for i in range(1, max)])
...
>>> zetaf(2, 10000) # ζ(2) の値。π^2/6に近いですね。（バーゼル問題　参照）
1.6448340618480597
>>> zetaf(10, 10000) # s の値を大きくすると 1 に近づく。
1.000994575127818
>>> zetaf(0.7, 10000) # s <= 1 では発散する（今回は max 項で止めてるから有限値）
50.05059218064199
>>> zetaf(9 + 2j, 10000) # s は複素数でもいい。
(1.0003243487185522-0.0019624012386066296j)
>>> zetaf(0.1 + 5j, 10000) # s は複素数でもいいけど、s の実部（Re s）が 1 以下だとこれも発散する
(611.6871790461571-490.2696093745057j)

　なかなか見つかりませんね。てかこのままじゃ無理です。ゼータ関数には定義域があり、その定義域の範囲外では値を持ちません。

\zeta (s) = \sum_{}^{} \frac{1}{n^s} = \frac{1}{1^s} + \frac{1}{2^s} + \frac{1}{3^s} + \dots　(定義域、Re(s) > 1)

　証明はしませんが、Re(s) <= 1 でゼータ関数は発散します（無限大になります）。有名なところで ζ(1) には調和級数という名前がついているのですが、非常に緩やかな速度で発散していきます。
　なお、Re(s) は複素数 s の実部って意味です。上記の範囲に零点がないのはグラフで見た方が分かりやすいかもですね。（ゼータ関数が定義される領域を薄く赤塗してます。）

上記のグラフを表示する Python のコード

　グラフの表示には matplotlib を使用します。ライブラリをインストールしていない場合は pip install matplotlib でインストールしてから下記のコードを入力してください。

riemann_zeta.py

__all__ = ["zetaf"]
__version__ = "0.1"

import matplotlib.pyplot as plt

def zetaf(s, max):
    return sum([1.0 / (i**s) for i in range(1, max)])

def make_zetaf_graph(s_max = 2):
    slist = [(i * 0.01 + 1) for i in range(1, s_max * 100)]
    zetavals = [zetaf(s, 500000) for s in slist]

    plt.rcParams["font.family"] = "MS Gothic"
    plt.title("ゼータ関数の挙動（50万項で打ち切り）")
    plt.xlabel("s の実部（虚部は 0）")
    plt.ylabel("ζ(s)")
    plt.axvspan(1, 1 + s_max, facecolor='red', alpha=0.15)
    plt.xlim(-1, 1 + s_max)
    plt.ylim(-1, 10)
    plt.axhline(0, linewidth=2, color="gray")
    plt.axvline(0, linewidth=2, color="gray")
    plt.plot(slist, zetavals, label="Zeta function")
    plt.grid()
    plt.legend()
    plt.show()

if __name__ == "__main__":
    make_zetaf_graph()

　ドキュメンテーション文字列は省略して、型ヒントもつけていませんが、実行自体は問題なくできます。（完全版はこの記事の最後にあります）

　ではどうやってゼータ関数の零点を見つければいいのでしょう？これには解析接続という理論を使用します。

　しかしこの記事の主題はリーマン予想を体験することなので数学的に難しいことは後回しにしておきます。難易度高めな数学のお話は、どっか別の記事で紹介します

　あまり難しくない範囲で、ゼータ関数の定義域をちょっとだけ拡張してみたいと思います。

0 < Re(s), s ≠ 0 に定義域を拡張（非自明な零点）

　導出過程は省略しますが、うまいこと変形すると下記の表示を得ることができます。

\zeta(s) = \frac{1}{s - 1} \sum_{n=1}^{\infty} \left( \frac{n}{(n+1)^s} - \frac{n - s}{n^s} \right)

　ゼータ関数が少しややこしくなっちゃいましたね。この形にしたときのメリットは「s の実部が 0 以上、かつ s ≠ 1 の範囲で発散せず値を持つ」というところです。下にグラフを載せますが、確かに新しいゼータ関数の定義の方ではちょっと定義域が広がっていますね。

上記のグラフを生成するプログラム

　riemann_zeta.py を修正する形でプログラムを完成させます。追記箇所にはコメントで「#追記」と書いてあります。また、プログラムは一部省略しています。

riemann_zeta.py【一部追記】

# 追記
__all__ = ["zetaf", "zetaf_other"]
__version__ = "0.1"

import matplotlib.pyplot as plt

def zetaf(s, max):
    return sum([1.0 / (i**s) for i in range(1, max)])

# 追記
def zetaf_other(s, max):
    return (1 / (s - 1)) * sum([(i / (i + 1)**s) - (i - s) / (i**s) for i in range(1, max)])

def make_zetaf_graph(s_max: int = 2):
    ~~~ 略 ~~~

# 追記
def make_zetaf_other_graph(s_max: int = 2):
    slist01 = [i * 0.01 for i in range(1, 100)]
    slist = [(i * 0.01 + 1) for i in range(1, s_max * 100)]
    zetavals = [zetaf(s, 500000) for s in slist]
    zetavals_o = [zetaf_other(s, 500000) for s in slist01] + [zetaf_other(s, 500000) for s in slist]

    plt.rcParams["font.family"] = "MS Gothic"
    plt.title("ゼータ関数の挙動（50万項で打ち切り）")
    plt.xlabel("s の実部（虚部は 0）")
    plt.ylabel("ζ(s)")
    plt.axvspan(0, 1, facecolor='red', alpha=0.1)
    plt.axvspan(1, 1 + s_max, facecolor='red', alpha=0.2)
    plt.xlim(-1, 1 + s_max)
    plt.ylim(-10, 10)
    plt.axhline(0, linewidth=2, color="gray")
    plt.axvline(0, linewidth=2, color="gray")
    plt.plot(slist, zetavals, label="Zeta function")
    plt.plot(slist01 + slist, zetavals_o, label="Other Zeta function")
    plt.grid()
    plt.legend()
    plt.show()

# 変更
if __name__ == "__main__":
    make_zetaf_other_graph()

　グラフでは元のゼータ関数（青線）とちょっとずれていますが、これはどちらのゼータ関数の表示も、足す項数を50万で打ち切っているからです。

\zeta(s) = \frac{1}{s - 1} \sum_{n=1}^{\infty} \left( \frac{n}{(n+1)^s} - \frac{n - s}{n^s} \right),　\zeta (s) = \sum_{n=1}^{\infty} \frac{1}{n^s}　　(理想)

\zeta(s) = \frac{1}{s - 1} \sum_{n=1}^{500000} \left( \frac{n}{(n+1)^s} - \frac{n - s}{n^s} \right),　\zeta (s) = \sum_{n=1}^{500000} \frac{1}{n^s}　　(現実)

　まぁ理想と現実の差は置いといて、とりあえずこの表示によってゼータ関数の定義域が広がりました！

　リーマン予想の主張は 「リーマンゼータ関数のすべての非自明な零点の実部は 1 / 2である。」 です。今見たゼータ関数の定義域は　0 < Re(s), s ≠ 0　なので、この新しいゼータ関数の形を利用すれば零点が見つかるわけです！

　ゼータ関数の零点は、見てわかる通り実数の領域だけでは見つかりません。s が複素数の時の値を見なきゃいけないわけですね。

　いきなりですが、s = 1/2 + i14.134725... は非自明な零点の一つです。これを式に代入してみます。

>>> def new_zetaf(s, max):
...     return (1 / (s - 1)) * sum([(i / (i + 1)**s) - (i - s) / (i**s) for i in range(1, max)])
...
>>> new_zetaf(0.5 + 14.134725j, 50000000)
(-5.812375642782399e-05-4.097708090742624e-05j)

　結果を見ると、確かにかなり 0 に近いですね！（残念ながら桁数が足りてないので正確に 0 にはなってないですが）

　ほかにも零点がありますが、一つずづ見るよりもグラフにした方が分かりやすいのでグラフにしてみます。ゼータ関数の引数が複素数なので、次のようにして３次元グラフにします。

\zeta(s) = \zeta(x, y) = \frac{1}{x + iy- 1} \sum_{n=1}^{\infty} \left( \frac{n}{(n+1)^{x + iy}} - \frac{n - (x + iy)}{n^{x + iy}} \right)

　x 軸、y 軸はわかりやすいですが、z 軸にはゼータ関数の絶対値を取りたいと思います。実際にはゼータ関数の値は複素数になりますが、今回着目したいのはζ関数の値が零になる所なのでこうします。

　下記の図が先ほどの式をグラフにしたものです。

グラフを表示するPythonプログラム

　riemann_zeta.py に追記します。わかりにくければ後回しにして、最後の完全版を見てください。

riemann_zeta.py【一部追記】

# ~~~略~~~
from mpl_toolkits.mplot3d import Axes3D

# ~~~略~~~

# 追記
def make_zetaf_3d():
    x_vals = np.linspace(0.01, 0.99, 100)
    y_vals = np.linspace(0.01, 40, 100)
    X, Y = np.meshgrid(x_vals, y_vals)
    Z = np.zeros_like(X, dtype=np.float64)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            s = complex(X[i, j], Y[i, j])
            try:
                val = zetaf_other(s, max=100)
                Z[i, j] = min(abs(val), 6)
            except Exception:
                Z[i, j] = np.nan
    fig = plt.figure(figsize=(10, 7))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap='YlGnBu', edgecolor='none', alpha=0.7)
    x_plane = 0.5
    y_plane = np.linspace(0.01, 40, 50)
    z_plane = np.linspace(0, 6, 50)
    Yp, Zp = np.meshgrid(y_plane, z_plane)
    Xp = np.full_like(Yp, x_plane)
    ax.plot_surface(Xp, Yp, Zp, color='firebrick', alpha=0.6, edgecolor='none')
    ax.set_xlabel('Re(s)')
    ax.set_ylabel('Im(s)')
    ax.set_zlabel('|ζ_other(s)|')
    ax.set_title('3D Plot of zetaf_other(s) with x = 1/2 Plane')
    fig.colorbar(surf, shrink=0.5, aspect=10)
    plt.show()


# 変更
if __name__ == "__main__":
    make_zetaf_3d()

　めっちゃ見づらいですね。でも、ちょっとだけ傾向をつかむことはできそうです。なんとなくしわがある絨毯みたいに見えますね？そしてしわの沈む部分はゼータ関数の絶対値が小さい＝０に近いです。
　そう、これがゼータ関数の零点（のヒント）です。

　なんかぱっと見、0 < Re(s) < 1 のすべてで零になっていそうですが、厳密には Re(s) = 0.5 でしか零になりません。

　今作った３次元グラフを真横から眺めてみます。

　とりあえず、上のグラフから零点の候補になりそうな点を視覚的に確認することができます（赤丸で囲っているところ）。過去の研究の結果で出ている零点をいくつか列挙しておきます。

\frac{1}{2} + i \cdot 14.134725 \dots

\frac{1}{2} + i \cdot 21.022040 \dots

\frac{1}{2} + i \cdot 25.010858 \dots

\frac{1}{2} + i \cdot 30.424876 \dots

\frac{1}{2} + i \cdot 32.935062 \dots

\frac{1}{2} + i \cdot 37.586178 \dots

　ちょうどさっきのしわの部分に近いですね！（Im(s) の値に注目！）
　ついでに、やっぱ３次元グラフだと見にくいので、s の実部が 1/2 の時のゼータ関数の実部・虚部・絶対値を見ておきましょう。真に興味があるのは絶対値の値で、下図の一番右のグラフです。

　図を見ると、一番右の絶対値が零になる点が、過去の研究で見つかっている零点とほとんど完全に一致しました！

　ちなみに現段階で非自明な零点は 0 < Re(s) < 1 の範囲にしかないことがわかっています。念のため、実部がほかの値でも確認しておきましょう。

• s = 1/3 + it のグラフ

• s = 1/√2 + it のグラフ

　確かにぱっと見、零点は無さそうに見えますが「零点は s = 1/2 + it 上にしかない！」と言い切るのは難しそうですね。（難しいどころじゃないか）

上記のグラフを表示するプログラム

　riemann_zeta.py に追記します。わかりにくければ後回しにして、最後の完全版を見てください。

riemann_zeta.py【一部追記】

# ~~~略~~~
from mpl_toolkits.mplot3d import Axes3D

# ~~~略~~~

# 追記
def make_zetaf_2d(re_s = 0.5):
    t_values = np.linspace(0.1, 40, 500)
    zeta_values = []
    for t in t_values:
        s = re_s + 1j * t
        z = zetaf_other(s, 100)
        zeta_values.append(z)
    real_parts = [z.real for z in zeta_values]
    imag_parts = [z.imag for z in zeta_values]
    abs_values = [abs(z) for z in zeta_values]
    plt.figure(figsize=(15, 5))
    plt.subplot(1, 3, 1)
    plt.plot(t_values, real_parts, color='blue')
    plt.axhline(0, color='black', linestyle='--', linewidth=1)  # y=0 ライン
    plt.grid(True)
    plt.title('Re(ζ(s)) vs t')
    plt.xlabel('t')
    plt.ylabel('Re(ζ(s))')
    plt.subplot(1, 3, 2)
    plt.plot(t_values, imag_parts, color='green')
    plt.axhline(0, color='black', linestyle='--', linewidth=1)
    plt.grid(True)
    plt.title('Im(ζ(s)) vs t')
    plt.xlabel('t')
    plt.ylabel('Im(ζ(s))')
    plt.subplot(1, 3, 3)
    plt.plot(t_values, abs_values, color='red')
    plt.axhline(0, color='black', linestyle='--', linewidth=1)
    plt.grid(True)
    plt.title('|ζ(s)| vs t')
    plt.xlabel('t')
    plt.ylabel('|ζ(s)|')
    plt.tight_layout()
    plt.show()


# 変更
if __name__ == "__main__":
    make_zetaf_2d()

s ≠ 0 の複素数全体に定義域を拡張（自明な零点）

　ところで非自明な零点は見つかりましたが、自明な零点はどこにあるのでしょうか？？
　自明な零点がどこにあるか知るにはゼータ関数をさらに解析接続して、もっと広範囲で定義される関数にする必要があります。
　これまた導出過程と証明は省略しますが、下記のように書くと非自明な零点の場所がわかります。
　（リーマンが見つけた式で関数等式とかって言われます。元論文とは形を変えています。）

\zeta(s) = 2^s \pi^{s-1} \sin (\frac{\pi s}{2}) \Gamma (1-s) \zeta (1-s)　　(定義域　s \neq 1)

　元のゼータ関数はどっか行っちゃいましたね。この式の形にすると、s = 1 + 0i でだけ発散しますがそれ以外のすべての複素数で有限の値をとることができます。なお、Γ はガンマ関数と呼ばれる関数で、階乗という概念を任意の複素数で成り立つように解析接続されたものです。ガンマ関数については次の式だけ知っておけばいいかなと思います。

\Gamma(n) = (n - 1)! = (n-1)\times(n-2)\times(n-3)\times\dots\times3\times2\times1

　そしてもう一つ重要な性質として、ガンマ関数は負の整数で発散します。
　上の式で着目したいのはゼータ関数の右辺に現れる sin(πs/2) という量です。sin は正弦関数と呼ばれ、下記の法則が成り立ちます。

\begin{align}
①　\sin(x + 2\pi) = \sin(x) \\
②　\sin(0) = \sin(\pi) = 0
\end{align}

　この法則から、sin(πs/2) は s = 2n (nは任意の整数) で 0 という値をとることがわかります。しかし、正の偶数（例えば s = 2）では 0 出ない有限の値が出ていますね。
　これは s が正の偶数の時、ガンマ関数は無限大に発散して、sin(πs/2) は 0 になるのでうまいこと打ち消しあうことにより、0 出ない有限の値をとるからです。（かなり正確さにかけた説明になってしまいますが、とにかく正の偶数では 0 にならないということ）

　というわけで、先ほどのゼータ関数は s = -2n (nは任意の自然数) で 0 になることが式の形から自明なのですね。

　以上のことから、ゼータ関数 ζ(s) の s = -2n (nは任意の自然数) の点たちを「自明な零点」と呼びます。

　今回、自明な零点を見るためのゼータ関数はプログラムで実装しません。（めんどっちぃのと負の領域の零点はグラフで見なくても自明なので...言い訳...）

プログラムの全貌

　今回のもろもろを実行するのに使用したファイルは次のようになっています。（ドキュメンテーションとかは省略して、コメントもほぼないです。型ヒントもない😥。いずれ修正して、細かい解説つけようと思います）

完成したプログラム！

riemann_zeta.py

"""ゼータ関数が定義されたモジュールです
"""

__all__ = ["zetaf", "zetaf_other"]
__version__ = "0.1"

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D  # 3D描画用
from scipy.special import gamma

def zetaf(s, max):
    return sum([1.0 / (i**s) for i in range(1, max)])

def zetaf_other(s, max):
    return (1 / (s - 1)) * sum([(i / (i + 1)**s) - (i - s) / (i**s) for i in range(1, max)])

def make_zetaf_graph(s_max = 2):
    slist = [(i * 0.01 + 1) for i in range(1, s_max * 100)]
    zetavals = [zetaf(s, 500000) for s in slist]

    plt.rcParams["font.family"] = "MS Gothic"
    plt.title("ゼータ関数の挙動（50万項で打ち切り）")
    plt.xlabel("s の実部（虚部は 0）")
    plt.ylabel("ζ(s)")
    plt.axvspan(1, 1 + s_max, facecolor='red', alpha=0.15)
    plt.xlim(-1, 1 + s_max)
    plt.ylim(-1, 10)
    plt.axhline(0, linewidth=2, color="gray")
    plt.axvline(0, linewidth=2, color="gray")
    plt.plot(slist, zetavals, label="Zeta function")
    plt.grid()
    plt.legend()
    plt.show()

def make_zetaf_other_graph(s_max: int = 2):
    slist01 = [i * 0.01 for i in range(1, 100)]
    slist = [(i * 0.01 + 1) for i in range(1, s_max * 100)]
    zetavals = [zetaf(s, 500000) for s in slist]
    zetavals_o = [zetaf_other(s, 500000) for s in slist01] + [zetaf_other(s, 500000) for s in slist]

    plt.rcParams["font.family"] = "MS Gothic"
    plt.title("ゼータ関数の挙動（50万項で打ち切り）")
    plt.xlabel("s の実部（虚部は 0）")
    plt.ylabel("ζ(s)")
    plt.axvspan(0, 1, facecolor='red', alpha=0.1)
    plt.axvspan(1, 1 + s_max, facecolor='red', alpha=0.2)
    plt.xlim(-1, 1 + s_max)
    plt.ylim(-10, 10)
    plt.axhline(0, linewidth=2, color="gray")
    plt.axvline(0, linewidth=2, color="gray")
    plt.plot(slist, zetavals, label="Zeta function")
    plt.plot(slist01 + slist, zetavals_o, label="Other Zeta function")
    plt.grid()
    plt.legend()
    plt.show()

def make_zetaf_3d():
    x_vals = np.linspace(0.01, 0.99, 100)
    y_vals = np.linspace(0.01, 40, 100)
    X, Y = np.meshgrid(x_vals, y_vals)
    Z = np.zeros_like(X, dtype=np.float64)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            s = complex(X[i, j], Y[i, j])
            try:
                val = zetaf_other(s, max=100)
                Z[i, j] = min(abs(val), 6)
            except Exception:
                Z[i, j] = np.nan
    fig = plt.figure(figsize=(10, 7))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap='YlGnBu', edgecolor='none', alpha=0.7)
    x_plane = 0.5
    y_plane = np.linspace(0.01, 40, 50)
    z_plane = np.linspace(0, 6, 50)
    Yp, Zp = np.meshgrid(y_plane, z_plane)
    Xp = np.full_like(Yp, x_plane)
    ax.plot_surface(Xp, Yp, Zp, color='firebrick', alpha=0.6, edgecolor='none')
    ax.set_xlabel('Re(s)')
    ax.set_ylabel('Im(s)')
    ax.set_zlabel('|ζ_other(s)|')
    ax.set_title('3D Plot of zetaf_other(s) with x = 1/2 Plane')
    fig.colorbar(surf, shrink=0.5, aspect=10)
    plt.show()

def make_zetaf_2d(re_s = 0.5):
    t_values = np.linspace(0.1, 40, 500)
    zeta_values = []
    for t in t_values:
        s = re_s + 1j * t
        z = zetaf_other(s, 100)
        zeta_values.append(z)
    real_parts = [z.real for z in zeta_values]
    imag_parts = [z.imag for z in zeta_values]
    abs_values = [abs(z) for z in zeta_values]
    plt.figure(figsize=(15, 5))
    plt.subplot(1, 3, 1)
    plt.plot(t_values, real_parts, color='blue')
    plt.axhline(0, color='black', linestyle='--', linewidth=1)  # y=0 ライン
    plt.grid(True)
    plt.title('Re(ζ(s)) vs t')
    plt.xlabel('t')
    plt.ylabel('Re(ζ(s))')
    plt.subplot(1, 3, 2)
    plt.plot(t_values, imag_parts, color='green')
    plt.axhline(0, color='black', linestyle='--', linewidth=1)
    plt.grid(True)
    plt.title('Im(ζ(s)) vs t')
    plt.xlabel('t')
    plt.ylabel('Im(ζ(s))')
    plt.subplot(1, 3, 3)
    plt.plot(t_values, abs_values, color='red')
    plt.axhline(0, color='black', linestyle='--', linewidth=1)
    plt.grid(True)
    plt.title('|ζ(s)| vs t')
    plt.xlabel('t')
    plt.ylabel('|ζ(s)|')
    plt.tight_layout()
    plt.show()

# 変更
if __name__ == "__main__":
    make_zetaf_2d()

まとめと後編

　以上で、何とかリーマン予想に入門することはできた気がします！（一応リーマン予想の意味は分かったし！）

　しかし最も重要な部分である、「ゼータ関数の非自明な零点の求め方」が抜けちゃってますね（さすがに「グラフにすればいいジャン！」は無理があります。）。

　後編の記事では、ゼータ関数の非自明な零点の求め方に焦点を絞っていきたいと思います。

　リーマン予想をPythonで体験してみよう！後編　作成中 . . .