プログラミング初心者でも素数公式を組んでみたい（後編） - 音楽室と化学室と美術室とPC室の融合部屋　所謂自由室

前回の『プログラミング初心者でも素数公式を組んでみたい』

aryuaryuaryuryu.hatenablog.com

の続きです．

今回はプログラムの話メインになると思います．

基本的にC++とDxlibを用いてます．

関数実装
参考文献

関数実装

前回得られた公式

$\pi(x) = \sum_{n=1}^{\infty}\frac{\mu(n)}{n}\left({\rm li}(x^\frac{1}{n}) - \sum_{\Im(\rho)\gt 0}\left({\rm li}(x^{^\frac{\rho}{n}}) + {\rm li}(x^{\frac{1-\rho}{n}})\right) + \int_{x^\frac{1}{n}}^{\infty}\frac{1}{t(t^2-1)\log{t}}dt - \log{2}\right)$

からプログラムを組んでみたい．

li以外の積分の実装

最初は簡単なものにチャレンジしたい．

以前作成した積分関数

aryuaryuaryuryu.hatenablog.com

を用いて

$\int_{x^\frac{1}{n}}^{\infty}\frac{1}{t(t^2-1)\log{t}}dt$

部分を計算したい．

これは被積分関数部分を組めば良いので

double RP(double t, std::vector<double> v) { return 1. / (t * (t * t - 1.) * log(t)); }

という感じでよいだろう．

グラフの形は

f:id:Aryuaryuaryuryu:20200224063408p:plain — 被積分関数のグラフ

このようになっており，1で無限大に発散，10ではもう殆ど0に近い形になってる．

どうせ $x$ は1超過でないといけないし，積分は20終了で十分な気がする．

メビウス関数の実装

次に組みやすそうな関数というとメビウス関数．

4で割れたら0を返す，2で割れるか調べて，3からひたすら奇数で小さい順に割っていくだけ．

平方数で割れたら0を返し割れた数で1を返すか-1を返すかを変えるだけ．

int MobiusFunc(int Num) {

int n = Num;

int count = 0;

if (n % 4 == 0)return 0;

if (n % 2 == 0) {

n /= 2;

count++;

}

for (int i = 3; i <= n; i+=2) {

if (n % (i * i) == 0)return 0;

if (n % i == 0) {

n /= i;

count++;

}

if (n == 1)break;

}

return (count % 2 == 0) ? 1 : -1;

}

対数積分の実装

積分関数実装してるとはいえ対数積分は結構面倒そうだなと思ったからいい方法無いかなと探したら指数積分を使う方法があるみたい．

正直勉強不足のためまだ理解が少々及ばないので殆どアテにしないように．

指数積分は

${\rm Ei}(x) = PV\int_{-x}^{\infty}\frac{e^{-t}}{t}dt$

と定義されるようで，これは補正指数積分

${\rm Ein}(x) = \int_{0}^{x} \frac{1-e^{-t}}{t}dt$

を定義すると， $1-e^{-t}$ のテイラー展開が $t(1+f(t))$ と表記できるので全平面で正則になる．

また

$E_1(x) = \int_{x}^{\infty}\frac{e^{-t}}{t}dt$

という関数を定義すると ${\rm Ei}(x) = -{\rm E_1}(-x) - \pi$ になり，更に

${\rm Ein}(x) - {\rm E_1}(x) - \log{x} = \gamma$

という関係が得られ，

${\rm Ei}(x) = \gamma - {\rm Ein}(x) + \log{x}$

になるそう．

で， ${\rm Ein}(x)$ は無限和表示にすると

${\rm Ein}(x) = -\sum_{k=1}^{\infty}\frac{(-z)^k}{k k!}$

になる．

この指数積分と対数積分は ${\rm li}(x) = {\rm Ei}(\log{x})$ という関係がある．

${\rm Ein}(x)$ のテイラー展開は

std::complex<double> Ein(std::complex<double> z) {

std::complex<double> result = 0;

for (int k = 1; k <= 300; k++) {

std::complex<double> tmp = 1. / (double)k;

for (int i = 1; i <= k; i++)tmp *= -z / (double)i;

result += tmp;

}

return result;

}

と実装ができ，ここから ${\rm Ei}(x)$ は

std::complex<double> Ei(std::complex<double> z) {

return M_GAMMA + log(z) - Ein(-z);

}

と実装ができる．こうしてようやく対数積分 ${\rm li}(x)$ を

std::complex<double> li(std::complex<double> z) { return Ei(log(z)); }

と実装できるようになった．

ちなみに今回は使ってないけど

$\sum_{n=1}^{\infty}\frac{\mu (n)}{n}{\rm li}(x^{\frac{1}{n}})$

は ${\rm li}(x) = {\rm Ei}(\log{x})$ という関係を用いて

$\sum_{n=1}^{\infty}\frac{\mu (n)}{n}{\rm li}(x^{\frac{1}{n}}) = \sum_{n=1}^{\infty}\frac{\mu (n)}{n}{\rm Ei}(\frac{1}{n}\log{x})\\\\ = \sum_{n=1}^{\infty}(\gamma + \log{\left(\frac{1}{n}\log{x}\right)} + \sum_{k=1}^{\infty}\frac{(\frac{1}{n}\log{x})^k}{k\cdot k!})$

という感じにバラすことができる．

途中の関係式は ${\rm Ei}(-x) = \gamma + \log{x} - {\rm Ein}(x) = \gamma + \log{x} + \sum_{k=1}^{\infty}\frac{(-z)^k}{k k!}$

によるもの．

から

$\sum_{n=1}^{\infty}\frac{\mu (n)}{n}{\rm li}(x^{\frac{1}{n}}) = \sum_{n=1}^{\infty}\frac{\mu(n)}{n}\gamma - \sum_{n=1}^{\infty}\log{n} + \sum_{n=1}^{\infty}\frac{\mu(n)}{n}\log{\log{x}} + \sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{\mu(n)}{n}\frac{(\log{x})^k}{n^k k\cdot k!}$

ここで

$\frac{1}{\zeta(s)} = \sum_{n=1}^{\infty}\frac{\mu(n)}{n^s}$ から得られる

ｋ $\sum_{n=1}^{\infty}\frac{\mu(n)}{n} = 0$ 及び， $\sum_{n=1}^{\infty}\frac{\mu(n)\log{n}}{n} = 1$ より

$\sum_{n=1}^{\infty}\frac{\mu (n)}{n}{\rm li}(x^{\frac{1}{n}}) = 1 + \sum_{k=1}^{\infty}\frac{(\log{x})^k}{k\cdot k!\zeta(k+1)}$

になる．

ゼータ関数の非自明な零点

ゼータ関数の非自明な零点はプログラムで計算させようとはしたけどあまり精度がよくならなかったし計算時間もバカにならないので予めできあがっているものを用意させてもらった．

３分クッキングばりの省略の仕方で申し訳ない．

非自明な零点部分の対数積分

正直ここがすごく苦労した．対数積分も指数積分も理解がままなっていないので，実装が非常につらく無事に死亡したので，先駆者のtujimotterさんはどうやってるのかなととりあえずざっくり調べてみた．

H. Riesel and G. Göhl, "Some Calculations Related to Riemann's Prime Number Formula", Mathematics of Computation, Vol.24, No.112, pp. 969-983 (1970).

という論文に書いてある式で実数のまま計算ができるようになるということで，よく見つけたなぁという感じでありがたく参考にさせていただく．

式の導出は

${\rm li}(x) = \int_{0}^{x}\frac{1}{\log{t}}dt$

なので， $e^{u+iv}$ を代入して $t = e^{z}$ と置換すると， $v\neq 0$ もしくは $u\lt 0$ のとき

${\rm li}(e^{u+iv}) = \int_{-\infty + iv}^{u+iv}\frac{e^z}{z}dz$

になる．

これは部分積分すると

${\rm li}(e^{u+iv}) = \int_{-\infty + iv}^{u+iv}\frac{e^z}{z}dz = \left[\frac{e^z}{z}\right]_{-\infty + iv}^{u+iv} + \int_{-\infty + iv}^{u+iv} \frac{e^z}{z^2}dz \\\\ = \frac{e^{u+iv}}{u+iv} + \int_{-\infty + iv}^{u+iv} \frac{e^z}{z^2}dz$

になって，

$\frac{e^{u+iv}}{u+iv}$ と $\int_{-\infty + iv}^{u+iv} \frac{e^z}{z^2}dz$ の比の絶対値を考えると， $|u+iv|\to \infty$ と吹き飛ばすと

$\left|\frac{\int_{-\infty + iv}^{u+iv} \frac{e^z}{z^2}dz}{\frac{e^{u+iv}}{u+iv}}\right|\to 0$

になるので，積分部分はほぼ無視をして

${\rm li}(e^{u+iv}) \sim \frac{e^{u+iv}}{u+iv}$

と近似ができるそう．

そこで

${\rm li}(x^{\rho}) = {\rm li}(e^{x\rho}) \sim \frac{e^{\rho\log{x}}}{\rho \log{x}}$

として，もしもリーマン予想が成り立つ，つまり $\rho = \frac{1}{2} + ia$ と表せられるなら

${\rm li}(x^{\rho}) \sim \frac{e^{(\frac{1}{2})\log{x}}}{(\frac{1}{2} + ia)\log{x}} = \frac{\sqrt{x}e^{ia\log{x}}}{(\frac{1}{2} + ia)\log{x}}$

と書くことができる．

これを有理化（？）して

${\rm li}(x^{\rho}) \sim \frac{\sqrt{x}e^{ia\log{x}} (\frac{1}{2} - ia)}{|\rho|^2 \log{x}}$

展開する

${\rm li}(x^{\rho}) \sim \sqrt{x}\frac{\frac{1}{2}\cos{(a\log{x})} + a\sin{(a\log{x})}}{|\rho|^2\log{x}} + i\sqrt{x}\frac{-a\cos{(a\log{x})} + \sin{(a\log{x})}}{|\rho|^2\log{x}}$

になる．

一方 $1-\rho = \frac{1}{2} - ia$ になるので

${\rm li}(x^{1 - \rho}) \sim \sqrt{x}\frac{\frac{1}{2}\cos{(a\log{x})} + a\sin{(a\log{x})}}{|\rho|^2\log{x}} - i\sqrt{x}\frac{-a\cos{(a\log{x})} + a\sin{(a\log{x})}}{|\rho|^2\log{x}}$