電卓でリーマンの素数個数関数を計算した

素数の個数を求める式

かなり昔から、ある数nより小さい素数の個数$p$は、次の簡単な式で近似できることが知られていた。

$$p\sim\frac{n}{log(n)}$$

素数と自然対数には深い関係があるらしい。それは何か。素数は無限にあるけれどこの関係はどこまでも続くのか。当時の学者もそう思って研究したようだが、なかなか解決できなかった。その中、1859年に数学者リーマンが次の式を発見した。

$$\pi(x)=\sum_{n=1}\frac{\mu(n)}{n}J(\sqrt[n]{x})　ただし　\sqrt[n+1]{x}<2　まで$$

$\pi(x)$は素数の個数を表す関数で、どんな大きな$x$に対しても$x$より小さい素数の個数が得られる。見れば見るほど不思議な式だ。$x$は自然数ではなくて実数、右辺には$n$乗根もあるのに、和は素数の個数という自然数に収束するというのだから。$\mu(n)$はメビウス関数という変な関数。$n$がどのような素数の積かによって$-1,0,1$のどれかの値になる。$J(x)$は、次のように定義されてる。

$$J(x)=Li(x)-\sum_{\rho}(Li(x^\rho)+Li(x^{1-\rho}))-\log2+\int_x^\infty\frac{dt}{t(t^2-1)\log t}$$

$Li(x)$や$\rho$など見慣れない記号があるが、調べてみると電卓でも計算できそう。2番めの項の$\rho$が複素数なので、複素数の関数計算ができる電卓が必要だ。式の詳細は理解できないから計算のための情報は書籍やネット頼み。途中が怪しくても、それなりの計算結果が得られればよし、としよう。

なぜ電卓？

実際に素数の個数を知ることは簡単だ。https://ja.wolframalpha.com/へ行って、$PrimePi(1000,000)$と入力すれば100万以下の素数が78498個あることがすぐに分かる。今回の興味は、今も数学者を悩ます歴史的な数式を電卓で計算する、ということにある。上の式の複素数$\rho$は、計算機の父チューリングも計算したとのこと。複素数で$\log(x)$や$e^x$の計算ができるHP35sを購入したので、摩訶不思議な数式を電卓で計算できたら面白い！というノリだ。電卓の複素数計算自体が面白くて、例えば虚数単位$i$の$i$乗を計算したのが右の写真。

$$i^i=e^{-\frac{\pi}{2}}\sim0.2078795763507619085469...$$

これは実数しかも超越数。

近似計算

$J(x)$の計算がプログラムの中心。第2項の$\rho$は無限個あることが証明されているから、加算も無限回でスパコンでも完全な計算はできない。それに電卓の計算速度は極端に遅いから、実際に加算できるのは10項以下だろう。正確な値を求めるプログラムにはならないが、そこは気楽な素人の遊びということで。

第１項　$Li(x)$

この項はネットで見つけた計算式で求めた。$\log$は自然対数なので電卓ではLNだ。この項は複素数ではないので普通の電卓でも計算できる。

$$Li(x)=\frac{x}{\log x}\sum_{k=0}^n\frac{k!}{\log(x^k)}$$

正しい値と比べてみると$x$の値に対して必要な精度が得られる項数$n$が異なる。プログラムでは必要な精度が得られるように$x$の値から$n$を適当に計算した。

第２項　$\sum Li(x^\rho)$

この項の計算が今回の中心だ。$\rho$はリーマンの論文に出てくるゼータ関数 $\zeta(s)$とやらがゼロとなる値で複素数。

$$\rho=\frac{1}{2}+i t　ただし　\zeta(\rho)=0$$

$\rho$の計算は大変複雑らしいが、先程のhttps://ja.wolframalpha.com/でN[ZetaZero[1],20]と入力すれば、１番目の$\rho_1$が20桁まで計算できる（感謝）。

$\rho_1=0.50000000000000000000 + 14.134725141734693791 i$

このサイトでプログラムに使う$\rho_1\sim\rho_8$を計算した。

$$Li(x^\rho)=Ei(\rho\log x)$$

$$Ei(z)=\frac{e^z}{z}\sum_{k=0}^n\frac{k!}{z^k}$$

第１項と似たような計算式だが、$\rho$が複素数なので35sの出番。調べた$\rho_1\sim\rho_8$の複素数値をそのまま使って関数計算する。$Ei(z)$の近似値を求める和は、実部はかなり正確だったが虚部は誤差が大きくなった。でも最終的な計算では実部の値しか使わないので助かった。リーマン先生は、全ての$\rho_n$の実部は$\frac{1}{2}$である、と予想したが、論文発表から160年が経過する現在でも証明されていない。証明に100万ドルの懸賞金がかかっている。

第３，４項

この２項は大きな値にならないとのことなので省略した。35sには数値積分の機能があるけれど時間がかかりすぎて使えそうにない。

プログラム

複雑な計算なので、プログラムは複数の部分に分けた。HP35sを持っていないと意味はないが、記録のために記載。こうして記載してみるとそれほど長いプログラムではない印象。

• $Ei(z)$の計算

E001	LBL E
E002	1
E003	XEQ E009	E009以下をサブプログラムとして実行
E004	RCL÷ Z
E005	RCL Z
E006	$e^x$
E007	×
E008	RTN
E009	STO O
E010	R↓
E011	STO Z
E012	1
E013	STO S
E014	STO K
E015	RCL÷ Z
E016	RCL× K
E017	STO+ S
E018	ISG K
E019	ENTER
E020	DSE O
E021	GTO E015
E022	RCL S
E023	RTN

• $Li(x)$の計算

L001	LBL L
L002	STO Y
L003	LN
L004	ENTER
L005	IP
L006	1
L007	-
L008	XEQ E009
L009	RCL× Y
L010	RCL Y
L011	LN
L012	÷
L013	RTN

• $J(x)$の計算

第2項の$\rho$による和は計算時間を勘案して8項までとした。

$\rho=\frac{1}{2}+it$は、最初の8個をメモリに入れ実行時に参照した。

$t=\{14.1347251417,\;21.0220396388,\;25.0108575801,\;\cdots\}$だ。

HP35sでは、$\rho_1$を $0.5\;i\;14.1347251417$ と入力する。

J001	LBL J
J002	LN
J003	STO X	変数Xに$\log{x}$を入れておく
J004	LASTx
J005	XEQ L001	$Li(x)$を計算
J006	STO Q	メモリQに格納
J007	0	$\sum_{\rho} Li(x^\rho)$計算開始
J008	STO T	和を入れる変数Tをクリア
J009	30.037	メモリの30〜37番地に$\rho$の値が8個格納されている
J010	STO I	ループカウンタを初期化
J011	RCL(I)	$\rho$（複素数）をメモリから得る
J012	RCL× X	$\rho\log{x}$を計算
J013	XEQ E001	$Li(x^\rho)$を$Ei(\rho\log{x})$で計算
J014	STO+ T	結果をTに加算
J015	ISG I	ループ終了を検査
J016	GTO J008
J017	RCL Q	$Li(x)$を呼び出す
J018	RCL T	$\sum(Li(x^\rho)+Li(x^{1-\rho}))=2Re(\sum Li(x^\rho))$
J019	2	虚部が相殺されるので和は実部の２倍
J020	×
J021	ENTER
J022	ABS	実部を計算（35sには実部を得る機能がない）
J023	x<>y	$Re(z)=|z|\cos(\arg(z))$
J024	ARG
J025	COS
J026	×
J027	-	$J(x)=Li(x)-\sum_{\rho}Li(x^\rho)$が得られる
J028	RTN

• $\pi(x)$の計算

メビウス関数$\mu(n)$は事前に計算した値をメモリに入れ実行時に参照。調べると$\pi(10^{12})$まで計算するには$\mu(40)$程度までの値が必要と思われたが、P026の短縮により実際の計算は$\mu(10)$程度で終了する。

P001	LBL P
P002	STO P
P003	0
P004	STO W
P005	STO V	ここまでは内部変数等の初期化
P006	40	メモリの40番地から$\mu(n)$の値が格納されている
P007	STO J
P008	1	$\pi(x)$計算のループ先頭
P009	STO+ V
P010	STO+ J
P011	RCL(J)	$\mu(n)$を得る
P012	x=0?
P013	GTO P009	$\mu(n)=0$なら$J(x)$は計算不要
P014	2
P015	RCL P
P016	RCL V
P017	PSE	途中経過を表示するため１秒停止
P018	$\sqrt[x]{y}$	$\sqrt[n]{x}$
P019	$x\le y$?	$\sqrt[n]{x}<2$なら計算終了
P020	GTO P030
P021	XEQ J001	$J(\sqrt[n]{x})$を計算
P022	RCL÷ V	$/n$
P023	RCL×(J)	$\frac{\mu(n)}{n}$
P024	STO+ W	W $=\sum \frac{\mu(n)}{n}J(\sqrt[n]{x})$
P025	ABS
P026	1.0	計算時間短縮のため 項<１ なら計算中止
P027	x<y?
P028	GTOP009
P029	RCL W
P030	ENTER
P031	INTG
P032	RTN

計算結果

プログラムの実行結果は次のとおり。肝心の$\rho$についての和は8項にすぎないが、$\rho$を使わない計算と比べると計算精度が改善されている。下表のように加算する項数が0（$\rho$を使わない）,5,8と増えると誤差は減少する。10億未満の素数の個数は5000万個以上あるが、それが9個の誤差で得られているのはいい感じだ。

$x$	$\pi(x)$	$\pi(x)$と計算結果との差			計算時間 （秒）
		0項	5項	8項
100	25	-1	-1	-1	35
1,000	168	-1	-2	-2	28
100,000	9,592	-4	1	1	35
1,000,000	78,498	29	19	19	42
1,000,000,000	50,847,534	-79	-66	9	50
10,000,000,000	455,052,511	-1,827	-817	-335	51
1,000,000,000,000	37,607,912,018	-1,475	-2,629	-3478	51

表末尾の$x=10^{12}$（10兆）では誤差の大小が逆転する。近似式の問題だろうか、それとも35sの有効桁数に近づいて計算誤差が累積しただろうか。計算時間は8項加算の時間だが、直線的に増えないのは、前述の$\mu(n)$関数の値によって計算量が不連続に変化することが影響していると思われる。途中経過の表示をやめると計算時間を5～10秒短縮できる。

印象と課題

HP35sのプログラム言語は、キーボードマクロに最小限の制御命令が追加されたコンパクトな仕様だが、メモリ参照命令が効果的で、素直に計算式をプログラム化できたことに感心した。計算速度もまあまあだ。課題は計算精度。近似式の検討、計算時間短縮のために省略した項や和の打ち切りなどを上手に処理できれば精度の向上が期待できる。次の機会に考えてみよう（いつになるやら）。

感想

素数の個数が実数や複素数の関数計算で求められることは、やはり不思議な感じがする。150年以上数学者の挑戦を跳ね返し続けているリーマン予想、そのややこしい式が手のひらで計算できるのを見たら、リーマン先生も面白がってくれるだろう。