非自明な零点の使い道
リーマン予想とは
「リーマンゼータ関数のすべての非自明な零点の実部は 1 / 2である。」
という予想です。この記事は、
「そもそもなんでリーマンゼータ関数の非自明な零点なんかに興味関心があるんだろ？」
って人のための記事です。要はこの記事で非自明な零点の使い道を説明したいって感じです。前提として「リーマンゼータ関数の非自明な零点が求まったら～」なので、実際に零点を求めたりしません。
シミュレーションに使用する言語は Python にしようと思います。Haskell も見やすくていいなと思ったんですが、ユーザの多さから Python を選択しました。コピペで動くように折り畳みの中にソースコード書いてます。
非自明な零点の求め方とかについては下記の記事を参考にしてください。
ある数以下の素数の個数
数学という学問の歴史は非常に長いですが、その長い歴史の中でも素数というテーマは常にと言っていいほど研究されてきました。この記事で取り上げるリーマン予想は「ある数以下の素数の個数」に関する研究で出てきた予想です。
ある整数 k 以下の素数の個数は慣習的に π(k) と書きます。円周率でよく使うあのパイです。いくつか具体例を挙げておきます。
\begin{align}
\pi(8) &= 4 (2, 3, 5, 7 の 四個が素数)\\
\pi(15) &= 6 (2, 3, 5, 7, 11, 13 の 六個が素数)\\
\pi(30) &= 10 (2, 3, 5, 7, 11, 13, 17, 19, 23, 29 の 十個が素数)
\end{align}
関数 π(k) の挙動はなんとなく理解しやすいですね。横軸を整数 k 、縦軸を π(k) として、整数 k が 100 までの区間でグラフ化すると次のようになります。見た目が階段みたいなので、素数階段とも呼ばれるグラフです。
上記を表示するプログラムはここ！
from math import isqrt
from typing import List
def count_primes_up_to(n: int) -> int:
"""
素数計数関数 π(n) を計算する。
引数:
n (int): 素数を数える上限の整数（n 以下の素数の個数を求める）
戻り値:
int: n 以下の素数の個数
例:
>>> count_primes_up_to(30)
10 # 素数: 2, 3, 5, 7, 11, 13, 17, 19, 23, 29
"""
if n < 2:
return 0
# エラトステネスの篩を使って素数を列挙
sieve: List[bool] = [True] * (n + 1)
sieve[0] = sieve[1] = False
for p in range(2, isqrt(n) + 1):
if sieve[p]:
for multiple in range(p * p, n + 1, p):
sieve[multiple] = False
return sum(sieve)
import matplotlib.pyplot as plt
from math import isqrt
from typing import List
def is_prime(k: int) -> bool:
"""
与えられた整数が素数かどうかを判定する。
引数:
k (int): 判定対象の整数
戻り値:
bool: 素数なら True、そうでなければ False
"""
if k < 2:
return False
for i in range(2, isqrt(k) + 1):
if k % i == 0:
return False
return True
def prime_staircase(n: int) -> List[int]:
"""
整数 k = 0 から n までの素数階段 π(k) をリストとして返す。
引数:
n (int): 上限値。0 から n までの素数階段を計算する。
戻り値:
List[int]: 各 k に対する π(k) の値のリスト（長さ n+1）
"""
count = 0
result: List[int] = []
for k in range(n + 1):
if is_prime(k):
count += 1
result.append(count)
return result
def plot_prime_staircase(n: int) -> None:
"""
素数階段 π(k) をグラフとして表示する。
引数:
n (int): 上限値。0 から n までの素数階段を描画する。
戻り値:
None
"""
staircase = prime_staircase(n)
x = list(range(n + 1))
y = staircase
plt.figure(figsize=(10, 5))
plt.step(x, y, where='post', color='blue', linewidth=2)
plt.title(f"k以下の素数の個数 π(k) のグラフ！素数階段と呼ばれます(k ≤ {n})", font="MS Gothic")
plt.xlabel("k")
plt.ylabel("π(k)")
plt.grid(True)
plt.tight_layout()
plt.show()
def main():
"""
作成した関数の動作確認はここでやっときます(^^♪
"""
print("整数", a := 100, f"以下の素数の個数は {count_primes_up_to(a)} 個です！")
plot_prime_staircase(a)
if __name__ == "__main__":
main()
（関数 π(x) は、正確には素数計数関数と呼ばれます。でも素数階段 π(x) って言い方が好きなので、この記事では素数階段で通します。）
たくさんの数学者がこの素数階段 π(k) を明示的な関数にしようと挑戦してきました。ちょっとそれらも見ておきましょう。
素数階段 π(x) の近似式
ここでは二つの π(k) の近似式を紹介します。ついでに今までは「ある整数以下の素数の個数」でしたが、実数でも当然正しく計算できる（例えば、π(30.1) は 30.1 以下の素数の個数って意味）ので、今後は実数らしさを強調するために π(x) と書きます。
素数階段の近似式として有名なのが下記の二つです。
\begin{align}
\pi(x) &\approx \frac{x}{\log x}&...①\\
\pi(x) &\approx \int_{2}^{x}\frac{dt}{\log t} = \int_{0}^{x}\frac{dt}{\log t} - \int_{0}^{2}\frac{dt}{\log t} = \mathrm{li}{(x)} - \mathrm{li}{(2)}&...②
\end{align}
これらの関数は厳密に π(x) と同じにはならないですが、かなりいい感じの挙動をします。グラフを見た方がわかりやすいですかね！
（補足。実は、①の公式は②の公式から導くことができます）
上記を表示するプログラムはここ！
コピペで実行するだけなら if __name__ == "__main__": とかいらないなって気づきました(^^♪
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from math import log
from typing import List
def prime_pi(x: float) -> int:
"""
素数階段関数 π(x): x 以下の素数の個数を返す。
Parameters:
x (int): 入力値（x ≥ 2）
Returns:
int: x 以下の素数の個数
"""
x = int(x)
if x < 2:
return 0
# エラトステネスの篩で素数を列挙
sieve = np.ones(x + 1, dtype=bool)
sieve[:2] = False
for i in range(2, int(x**0.5) + 1):
if sieve[i]:
sieve[i*i : x+1 : i] = False
return np.count_nonzero(sieve)
def gauss_approx(x: float) -> float:
"""
ガウス近似 π(x) ≈ x / log(x)
Parameters:
x (float): 入力値（x > 1）
Returns:
float: 近似値
"""
return x / log(x) if x > 1 else 0.0
def riemann_approx(x: float) -> float:
"""
リーマン近似 π(x) ≈ ∫₂ˣ dt / log(t)
Parameters:
x (float): 入力値（x ≥ 2）
Returns:
float: 近似値（数値積分）
"""
if x < 2:
return 0.0
result, _ = quad(lambda t: 1 / log(t), 2, x)
return result
# x軸の値（整数）
x_vals = np.arange(2, 10001)
# 各関数の値を計算
pi_vals = [prime_pi(x) for x in x_vals]
gauss_vals = [gauss_approx(x) for x in x_vals]
riemann_vals = [riemann_approx(x) for x in x_vals]
# グラフ描画
plt.figure(figsize=(12, 6))
plt.plot(x_vals, pi_vals, label=r"$\pi(x)$ (Prime count)", color='black', linewidth=2)
plt.plot(x_vals, gauss_vals, label=r"$\frac{x}{\log x}$ (Gauss)", color='blue', linestyle='--')
plt.plot(x_vals, riemann_vals, label=r"$\int_2^x \frac{dt}{\log t}$ (Riemann)", color='green', linestyle='-.')
plt.title("素数階段とその近似", font="MS Gothic")
plt.xlabel("x")
plt.ylabel("値", font="MS Gothic")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
結果を見てみると、②の公式の方がだいぶいい精度で素数階段 π(x) を近似できていることがわかります！ちなみに②に出てくる li(x) は対数積分と呼ばれています。残念ながら対数積分は初等関数（三角関数とか指数関数とか）で表現できません。
さぁ、うまいこと近似しているのはわかりますが、これらは結局のところ近似式です。厳密にあのカクカクな素数階段について、近似 ≂ ではなく等号 = で書けるような式はないのでしょうか！
リーマンの素数公式
というわけで本題です。リーマンはかの有名な論文「与えられた数より小さな素数の個数について」の中で、現在ではリーマンの素数公式という名がついている π(x) に関する式を求めました。
その導出過程が非常に美しいのですが、複雑な数学の話はイタズラに読む気を失せさせるので一気に答えだけ出します。
いろんな表示の仕方があるんですが、ここでは下記のように表示します。Wikipediaがこれだったので^^
\pi(x) = \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \left( \mathrm{li}(x^{1/m}) - \sum_{\rho} \mathrm{li}(x^{\rho/m}) - \log 2 + \int_{x^{1/m}}^{\infty} \frac{dt}{t(t^2 - 1)\log t} \right)
素数階段 π(x) がイコールで結ばれてます！ま、さっきまで見た素数公式より数段難しくなりましたがね！！がはは
式の意味を説明する前に、さすがにもう少しわかりやすくしたいのでリーマンの素数公式を分解します。と言っても、単純に分配法則でかっこを外すだけですが。
\begin{align}
\pi(x) &= \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \left( \mathrm{li}(x^{1/m}) - \sum_{\rho} \mathrm{li}(x^{\rho/m}) - \log 2 + \int_{x^{1/m}}^{\infty} \frac{dt}{t(t^2 - 1)\log t} \right)\\
&= \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \mathrm{li}(x^{1/m}) - \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \sum_{\rho} \mathrm{li}(x^{\rho/m}) + \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \left(- \log 2 + \int_{x^{1/m}}^{\infty} \frac{dt}{t(t^2 - 1)\log t} \right)\\
&= R(x) - S(x) + I(x)
\end{align}
最後の文字の置き換えは下記の式で対応しています。
\begin{align}
R(x) &= \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \cdot \mathrm{li}(x^{1/m}) &...(A)\\
S(x) &= \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \sum_{\rho} \mathrm{li}(x^{\rho/m}) &...(B)\\
I(x) &= \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \left(- \log 2 + \int_{x^{1/m}}^{\infty} \frac{dt}{t(t^2 - 1)\log t} \right) &...(C)
\end{align}
では、式の意味を探っていきましょう。各項に (A), (B), (C) とつけたので、その順番に見ていきます。
(A)の項 R(x)
R(x) = \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \cdot \mathrm{li}(x^{1/m})
意味不明なのは μ(x) ですかね。これはメビウス関数と呼ばれる関数です。詳しくはWikipediaでも見てもらって、ここでは簡単にメビウス関数の定義だけ話しておきます。
メビウス関数 μ(x) は、「x を素因数分解したときに、どこかに素数の 2 乗以上の項があれば 0、ない場合は構成する素数の個数が奇数なら -1 、偶数なら 1 を返す」ような関数です。わかりにくいですね、具体例を出しておきます。
\begin{align}
\mu(28) &= 0&(28=2^2 \times 7 で素数の2乗があるから 0)\\
\mu(15) &= 1&(15=3 \times 5 で偶数個の素数の積だから 1)\\
\mu(30) &= -1&(30=2 \times 3 \times 5 で奇数個の素数の積だから -1)\\
\end{align}
って感じです。これを踏まえて、ちょっと R(x) の具体的な値を計算してみましょう。ためしに R(30) を計算します。
\begin{align}
R(30) &= \sum_{m=1}^{\log_2 30} \frac{\mu(m)}{m} \cdot \mathrm{li}(30^{1/m})\\
&= \sum_{m=1}^{4} \frac{\mu(m)}{m} \cdot \mathrm{li}(30^{1/m}) (\log_2 30=4.907...なので足すのは4項だけ)\\
&= \frac{\mu(1)}{1} \cdot \mathrm{li}(30^{1/1}) + \frac{\mu(2)}{2} \cdot \mathrm{li}(30^{1/2}) + \frac{\mu(3)}{3} \cdot \mathrm{li}(30^{1/3}) + \frac{\mu(4)}{4} \cdot \mathrm{li}(30^{1/4})\\
&(メビウス関数の定義、\mu(1) = 1、\mu(2) = -1、\mu(3) = -1、\mu(4) = 0 を代入)\\
&= \mathrm{li}(30^{1/1}) - \frac{1}{2} \cdot \mathrm{li}(30^{1/2}) - \frac{1}{3} \cdot \mathrm{li}(30^{1/3}) + \frac{0}{4} \cdot \mathrm{li}(30^{1/4})\\
&= 13.02 - \frac{1}{2} \cdot 3.923 - \frac{1}{3} \cdot 2.260 (有効桁４桁で切ってます)\\
&= 10.30 (これが答え！)
\end{align}
実際の素数階段の答えは π(30) = 10 なことを考えるとだいぶ近いですね！実際、この R(x) だけですでに相当イイ近似を得ることができます。
上記を表示するプログラムはここ！
R(x)を表示するプログラム
import numpy as np
import matplotlib.pyplot as plt
from sympy import mobius, li
import math
def prime_pi(x):
if x < 2:
return 0
sieve = [True] * (x + 1)
sieve[0:2] = [False, False]
for i in range(2, int(x**0.5) + 1):
if sieve[i]:
for j in range(i*i, x + 1, i):
sieve[j] = False
return sum(sieve[:x+1])
def R(x):
if x < 2:
return 0
total = 0
for m in range(1, int(math.log2(x)) + 1):
mu = mobius(m)
if mu != 0:
total += mu / m * float(li(x ** (1 / m)))
return total
# x軸の値（0〜100）
x_vals = np.arange(2, 101)
r_vals = [R(x) for x in x_vals]
pi_vals = [prime_pi(x) for x in x_vals]
# グラフ描画
plt.rcParams['font.family'] = 'MS Gothic'
plt.figure(figsize=(10, 6))
plt.step(x_vals, pi_vals, label=r'$\pi(x)$ (厳密な素数階段)', color='blue', linewidth=2, where='post')
plt.plot(x_vals, r_vals, label=r'$R(x)$ (主要項のみ)', color='crimson', linestyle='--', linewidth=2.5, marker='o', markersize=4)
plt.xlabel('x', fontsize=12)
plt.ylabel('素数の個数', fontsize=12)
plt.title('素数個数関数 π(x) とリーマン主要項 R(x) の比較', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()
この時点でだいぶイイ感じですが、まだカクカクした階段を再現できていないですね！
(B)の項 S(x)
S(x) = \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \sum_{\rho} \mathrm{li}(x^{\rho/m})
見た感じは(A)の式と似ていますね。実際、計算のやり方は(A)の方と変わらないです。ただし(A)の時と決定的に違うのが、対数積分 li が変数 ρ の和になっています。
実はこの変数 ρ こそがリーマンゼータ関数の非自明な零点なのですね。非自明な零点全てで和をとるって意味です。
さて早速プログラムに落としましょう！って思ったけどちょっと大変なんですね、この計算。リーマンゼータ関数の非自明な零点って複素数だし、ちょっと自分のくそ雑魚PCで愚直に計算するのは大変なので、下記のような近似を行いたいと思います。
まずメビウス関数についての和の部分に着目します。
\sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m}
これ、m がデカくなればなるほど終盤の寄与はだいぶ小さくなりますね。さらに対数積分の引数の x^(ρ/m) も、 m がデカくなればなるほど小さくなります。なので、 x が十分デカい時には m = 1 の寄与が支配的になるのですね。
今回はこれを利用して、下記のような大胆な近似をします。
S(x) = \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \sum_{\rho} \mathrm{li}(x^{\rho/m}) \approx \sum_{\rho} \mathrm{li}(x^{\rho}) (m=1だけ使用)
ついでに対数積分 li(x) も大胆に漸近展開して近似して、最終的に下記のようにします。
S(x) = \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \sum_{\rho} \mathrm{li}(x^{\rho/m}) \quad \approx \quad -\sum_{\rho} \frac{x^{\rho}}{\rho \log x}
これでも、ある程度大きい x では十分に動きを確認できます。下記でグラフを書いておきます。（零点最初の〇〇個って書いてありますが、正確には共役な零点もあるので〇〇×２個の零点を使用しています）
ちなみに、すでに見つかっている〇点のリストは下記のサイトから見ることができます！
- 零点最初の30個
- 零点最初の60個
- 零点最初の100個
上記を表示するプログラムはここ！
import numpy as np
import matplotlib.pyplot as plt
from sympy import mobius, li
from mpmath import zetazero, li as mpli
import math
from sympy import mobius, li, I
from sympy import re, im
# 厳密な素数個数関数 π(x)
def prime_pi(x):
if x < 2:
return 0
sieve = [True] * (x + 1)
sieve[0:2] = [False, False]
for i in range(2, int(x**0.5) + 1):
if sieve[i]:
for j in range(i*i, x + 1, i):
sieve[j] = False
return sum(sieve[:x+1])
# 主要項 R(x)
def R(x):
if x < 2:
return 0
total = 0
for m in range(1, int(math.log2(x)) + 1):
mu = mobius(m)
if mu != 0:
total += mu / m * float(li(x ** (1 / m)))
return total
zeta_zeros = [
14.134725141, 21.022039639, 25.010857580, 30.424876125, 32.935061588,
37.586178159, 40.918719012, 43.327073281, 48.005150881, 49.773832478,
52.970321478, 56.446247697, 59.347044003, 60.831778525, 65.112544048,
67.079810529, 69.546401711, 72.067157674, 75.704690699, 77.144840069,
79.337375020, 82.910380854, 84.735492981, 87.425274613, 88.809111208,
92.491899271, 94.651344041, 95.870634228, 98.831194218, 101.317141358,
103.725538040, 105.446623052, 107.168611184, 111.029535543, 111.874659176,
114.320220915, 116.226680321, 118.790782865, 121.370125002, 122.946829294,
124.256818554, 127.516683880, 129.578704199, 131.087688530, 133.497737203,
134.756509753, 138.116042055, 139.736208952, 141.123707404, 143.111845808,
146.000982487, 147.422765343, 150.053520421, 150.925257612, 153.024693811,
156.112909294, 157.597591817, 158.849988171, 161.188964138, 163.030709687,
165.537069188, 167.184439978, 169.094515416, 169.911976479, 173.411536519,
174.754191523, 176.441434297, 178.377407776, 179.916484020, 182.207078484,
184.874458461, 185.598783678, 187.228922584, 189.416158656, 192.026656361,
193.630345505, 195.265396679, 196.876481841, 198.831945345, 201.264751944,
202.493594514, 204.189671803, 205.394697202, 207.906258888, 209.576509716,
211.690862595, 213.347919360, 214.547044783, 216.169538508, 219.067596349,
220.714918839, 222.426526841, 224.007000255, 226.378377154, 228.873973029,
229.337413306, 231.250188700, 234.265375456, 236.524229666, 237.769820548,
]
# S(x) の定義
def S(x):
sum_val = 0
for gamma in zeta_zeros:
rho1 = 0.5 + I * gamma
rho2 = 0.5 - I * gamma
try:
term1 = x ** rho1 / (rho1 * np.log(x))
term2 = x ** rho2 / (rho2 * np.log(x))
sum_val += re(term1 + term2)
except Exception:
continue
return -sum_val
# x軸の値
x_vals = np.arange(2, 151)
# 各関数の値を計算
pi_vals = [prime_pi(x) for x in x_vals]
r_vals = [R(x) for x in x_vals]
# R(x) + S(x) の計算
rs_vals = [R(x) + S(x) for x in x_vals]
# グラフ描画
plt.rcParams['font.family'] = 'MS Gothic'
# グラフに追加
plt.figure(figsize=(10, 6))
plt.step(x_vals, pi_vals, label=r'$\pi(x)$ (厳密な素数階段)', color='blue' , linewidth=2, where='post')
plt.plot(x_vals, r_vals, label=r'$R(x)$ (主要項のみ)', color='crimson', linestyle='--', linewidth=2.5)
plt.plot(x_vals, rs_vals, label=r'$R(x) + S(x)$ (零点補正付き)', color='green', linewidth=3.0)
plt.xlabel('x', fontsize=12)
plt.ylabel('素数の個数', fontsize=12)
plt.title('素数個数関数 π(x)、R(x)、R(x)+S(x) の比較(零点100個)', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()
驚くほど元の階段関数 π(x) に近づいていきますね！！
(C)の項 I(x)
最後！(C)の項である I(x) ですが、次のような定義になっています。
I(x) = \sum_{m=1}^{\log_2 x} \frac{\mu(m)}{m} \left(- \log 2 + \int_{x^{1/m}}^{\infty} \frac{dt}{t(t^2 - 1)\log t} \right)
ややこしいしきっと大事な項目なのでしょう！って思いそうですが、こいつはリーマンの素数公式と階段関数を厳密な等式とするための誤差項としての見方が正しいです。誤差です。まじで。
一応この式を適用した場合のグラフを下記に示します。
ほぼ変わらんですね。差が見えないのデ拡大してみると、若干よくなってるように見える気がします。
まぁ微差も微差なので、ここでは誤差みたいなものとして紹介しました。
プログラムはこちら！
import numpy as np
import matplotlib.pyplot as plt
import math
from sympy import mobius, li, I, re
from mpmath import quad, log as mplog, mp
mp.dps = 30 # 高精度設定
def prime_pi(x: int) -> int:
"""
厳密な素数個数関数 π(x) をエラトステネスの篩で計算する。
x 以下の素数の個数を返す。
"""
if x < 2:
return 0
sieve = [True] * (x + 1)
sieve[0:2] = [False, False]
for i in range(2, int(x**0.5) + 1):
if sieve[i]:
for j in range(i*i, x + 1, i):
sieve[j] = False
return sum(sieve[:x+1])
def R(x: float) -> float:
"""
リーマンの主要項 R(x) を計算する。
Möbius 関数と対数積分 li(x) を用いた近似。
"""
if x < 2:
return 0
total = 0
for m in range(1, int(math.log2(x)) + 1):
mu = mobius(m)
if mu != 0:
total += mu / m * float(li(x ** (1 / m)))
return total
def S(x: float) -> float:
"""
零点補正項 S(x) を計算する。
Riemannゼータ関数の非自明な零点を用いて振動項を加える。
"""
sum_val = 0
for gamma in zeta_zeros:
rho1 = 0.5 + I * gamma
rho2 = 0.5 - I * gamma
try:
term1 = x ** rho1 / (rho1 * np.log(x))
term2 = x ** rho2 / (rho2 * np.log(x))
sum_val += re(term1 + term2)
except Exception:
continue
return -sum_val
def I_term(x: float) -> float:
"""
積分補正項 I(x) を計算する。
対数積分の補正として、無限積分を数値的に評価する。
"""
if x < 2:
return 0
total = 0
for m in range(1, int(math.log2(x)) + 1):
mu = mobius(m)
if mu == 0:
continue
lower = x ** (1 / m)
try:
integral = quad(lambda t: 1 / (t * (t**2 - 1) * mplog(t)), [lower, mp.inf])
term = mu / m * (-mp.log(2) + integral)
total += float(term)
except Exception:
continue
return total
zeta_zeros = [
14.134725141, 21.022039639, 25.010857580, 30.424876125, 32.935061588,
37.586178159, 40.918719012, 43.327073281, 48.005150881, 49.773832478,
52.970321478, 56.446247697, 59.347044003, 60.831778525, 65.112544048,
67.079810529, 69.546401711, 72.067157674, 75.704690699, 77.144840069,
79.337375020, 82.910380854, 84.735492981, 87.425274613, 88.809111208,
92.491899271, 94.651344041, 95.870634228, 98.831194218, 101.317141358,
103.725538040, 105.446623052, 107.168611184, 111.029535543, 111.874659176,
114.320220915, 116.226680321, 118.790782865, 121.370125002, 122.946829294,
124.256818554, 127.516683880, 129.578704199, 131.087688530, 133.497737203,
134.756509753, 138.116042055, 139.736208952, 141.123707404, 143.111845808,
146.000982487, 147.422765343, 150.053520421, 150.925257612, 153.024693811,
156.112909294, 157.597591817, 158.849988171, 161.188964138, 163.030709687,
165.537069188, 167.184439978, 169.094515416, 169.911976479, 173.411536519,
174.754191523, 176.441434297, 178.377407776, 179.916484020, 182.207078484,
184.874458461, 185.598783678, 187.228922584, 189.416158656, 192.026656361,
193.630345505, 195.265396679, 196.876481841, 198.831945345, 201.264751944,
202.493594514, 204.189671803, 205.394697202, 207.906258888, 209.576509716,
211.690862595, 213.347919360, 214.547044783, 216.169538508, 219.067596349,
220.714918839, 222.426526841, 224.007000255, 226.378377154, 228.873973029,
229.337413306, 231.250188700, 234.265375456, 236.524229666, 237.769820548,
]
# x軸の値
x_vals = np.arange(2, 151)
# 各関数の値を計算
pi_vals = [prime_pi(x) for x in x_vals]
r_vals = [R(x) for x in x_vals]
# R(x) + S(x) の計算
rs_vals = [R(x) + S(x) for x in x_vals]
rsi_vals = [R(x) + S(x) + I_term(x) for x in x_vals]
# グラフ描画
plt.rcParams['font.family'] = 'MS Gothic'
# グラフに追加
plt.figure(figsize=(10, 6))
plt.step(x_vals, pi_vals, label=r'$\pi(x)$ (厳密な素数階段)', color='blue' , linewidth=2, where='post')
plt.plot(x_vals, r_vals, label=r'$R(x)$ (主要項のみ)', color='crimson', linestyle='--', linewidth=2.5)
plt.plot(x_vals, rs_vals, label=r'$R(x) + S(x)$ (零点補正付き)', color='green', linewidth=3.0)
plt.plot(x_vals, rsi_vals, label=r'$R(x) + S(x) + I(x)$ (リーマンの素数公式)', linewidth=3.0)
plt.xlabel('x', fontsize=12)
plt.ylabel('素数の個数', fontsize=12)
plt.title('素数個数関数 π(x)、R(x)、R(x)+S(x) の比較(零点100個)', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True)
plt.tight_layout()
plt.show()
まとめ
リーマンってすごい