このエントリでは、Riemannのζ関数の自明でない零点と素数の分布がどう関係し、非自明な零点は素数分布に対してどんな機能を持っているか、について述べる。つまり、一般人の一番の疑問点である、Riemannのζ関数と素数分布関数Pi関数を直接見てみたい。そしてこのことはRiemann予想が成立しているかどうかに関係なく成り立つ。今回のエントリについては、Stan Wagon "Mathematica in Action のChap 25に多くを負っていて僕自身のオリジナルはほとんどない。Riemann予想で遊んでいるとよくわかるのだが、自分の手に負えない問題は巨人の肩に乗らないと何にもできないのだ。
ちなみにRiemannのζ関数の自明でない零点とか非自明な零点とかいう言い方について一言。 ζ関数の自明な零点とは、ζ関数のre s<0の領域においては、s=0, -2, -4, -6...という点でζ関数が零になり、それ以外に零点が存在しないということである。それに対して、自明でない零点だとか非自明な零点とは、それ以外の零点のことで、この非自明な零点について述べたものがRiemann予想なのだ：「ζ関数の自明でない零点の実数部はすべて1/2である」。
さて、このエントリの目的のためには、Pi関数とζ関数を結ぶRiemann-von Mangoldtの明示公式と呼ばれるものが必要だ。というか、この明示公式を出しちゃえばもう終わりである。だって、Riemann-von Mangoldtの明示公式は、Pi関数とζ関数の零点を直接繋げるものであるからだ。
この公式にあるμはメビウス関数、J関数のliとは対数積分であり、ζはζ関数である。J関数の第二項の和はζ関数の自明でない零点について和をとるということであり、これこそがζ関数の自明でない零点と素数分布関数Piを結ぶものなのだ。
この公式はRiemannにより与えられ、von Mongoldtにより証明された。Riemann-von Mongoldtの明示公式の導出と証明はEdwards, H. M.の"Riemann's Zeta Function" Chap3にあるのでここでは議論しない、なんていってしまうと、このエントリの存在意義自体がなくなってしまうのだが、そう言いたくなるのは理由がある。明示公式の導出と証明はそれなりに長く、関数論の熟達した計算技術を求められる。僕の数学能力ではこの公式を導出するだけでも死ぬかと思ったくらいなので、人に説明することなど困難だ。無理に説明するとなると、EdwardsのChap3をコピるだけになってしまうので、意味がなくなってしまうという次第だ。
そこで今回の目的は、このRiemann-von Mangoldtの明示公式を既知のモノとしてそれを巨人の肩とし、これが実際にPi関数を計算するものであるかどうかというのを、Mathematicaを用いて数値計算をしてみるというように小さく設定する。Riemann-von Mangoldtの明示公式で遊ぶことで、素数の分布を表すPi関数とζ関数の非自明な零点の関係を直感的につかもうというのが主旨である。
さて、まずJ関数をPi関数に直接展開した式から始める。
このままでは数値計算を行うのはなかなか難しい。そこで、まず第一項について、 Riesel, H.の "Prime Numbers and Computer Methods for Factorization" の結果を利用して、対数積分をζ関数の級数として表した表示式に変換する。
また第二項については、実はMathematicaの対数積分の実装LogIntegralは複素数部分ではうまく働かないのだ。そこで、複素数でも働くエクスポネンシャル積分関数eiで書き換える。またζ関数の自明でない零点rhoは虚数~rhoとペアで現れる、つまり複素共役もζ関数の零点になっている、という事実があり、これを利用すると実数部分を考えるだけでよい。したがって、実数をパートを取る操作reを行って最後に二倍すればよいのだ。
更に第三項、第四項には、Riesel, H. & Goehl, G., p969-983, Mathematics of Computation, Vol24, 1970により、nについて154項までの和をとると、次のようにarctanの関数として近似ができることを用いる。
上記の近似式を各項に適用した結果、Pi関数は次のような近似式として数値計算に適した形で表示することができる。
ここまでくれば、Mathematicaのコードに直すのは比較的容易である。まず、第一項を関数として表現すると次のようになる。ここでは、無限項の和をとることができないので500までの和をとることにする。
RiemannR[x_] :=
1 + (Floor[N[Log[x]]^Range[500]]).Map[1/(#! # N[Zeta[# + 1]]) &,
Range[500]];
RiemannR[0] := 0;
Attributes[RiemannR] =\ {Listable\};
次に第二項についてζ関数の零点について和を取る前の式を、xとζ零点のインデックスを引数とする関数として表現する。ここではnについて154項までの和をとるのだが、これは先の第三項、第四項の計算においてRiesel&Goehlの近似式を用いたからだ。実際のコードにおいて和をとる操作は次のように行う。μ/nの項とeiの項の級数をそれぞれベクトルとして計算しておいてから、その二つの内積をとることにする。メビウス関数μにはMathematicaの組込関数MoebiusMuを用い、エクスポネンシャル積分関数eiには組込関数ExpIntegralEiを用いる。
P[x_, k_] := -2*
Table[MoebiusMu[n]/n, {n, 1, 154}].Table[
Re[ExpIntegralEi[(ZetaZero[k] Log[N[x]])/n]], {n, 1, 154}];
Attributes[P] = {Listable};
この第二項には問題のζ関数の非自明な零点が含まれている。これがまさにζ関数と素数分布のPi関数を繋ぐ部分だ。これをどのように処理するかというと、MathematicaにはなんとRiemannのζ関数の自明でない零点を与えてくれる関数が組込関数ZetaZeroとしてあるのだ。だから簡単にこれを利用すればよいだけである。しかし、バージョンによってはこの組込関数を用いることができないのだが、心配はいらない。Riemannのζ関数の零点計算の権威である Andrew Odlyzkoが 零点テーブル を用意してくれているのだ。だからたとえ組込関数としてなかったとしてもこのテーブルを利用すればいいのだ。つまり、ここでも巨人の肩に乗るのだ。
ここで忘れないうちに第一項、第三項の和をとっておく。ここで計算を行うための定義域は、仮に10から100としておく。
domain = Range[10., 100, 90/450]; RiemannTerms = RiemannR[domain] + ArcTan[Pi / Log[domain]] / Pi];
いよいよメインの第二項の計算を行う。先ほど定義した関数Pに対するζ関数の非自明な零点の和をとればよい。明示公式では無限個の非自明な零点の和をとることになっているのだが、そいつは無理な相談なので、とりあえず50個までの零点の和をとることにする。この50個の和を取る計算にしても、Xeonのかなり早いデュアルプロセッサでさえ25分くらい必要だ。
PeriodicTerm = Sum[P[domain, i], {i, 1, 50}];
最後に、すべての結果を足して実際のPi関数と比較する。Pi関数はPrimePiという組込関数をそのまま利用する。
ListPlot[{Transpose[{domain, RiemannTerms + PeriodicTerm}],
Transpose[{domain, PrimePi[domain]}]}, Joined -> True,
PlotRange -> {{5, 100}, {0, 25}}, AxesOrigin -> {10, 4}]
これを表示すると次のようになる。
このように、それなりの精度でPi関数がζの非自明な零点を用いた近似計算として表現されていることがわかる。この計算過程をみていくと、ζの非自明な零点が無限個与えられれば、素数の分布が表現できるという、Riemann-von Mangoldtの明示公式の意義がよくわかる。上記の計算ではζ関数の非自明な零点は最初の500個しかとらなかったけれど、これを多くすればするほど近似の精度が上がるのだ。そして、それは単なる近似というには及ばない。Littlewoodの計算によれば、この零点が絡む誤差項こそが主要な成分になる部分が出てくるのだ。
重要なことは、この明示公式自体はRiemann予想が真か偽であるかに関係なく成立するということだ。しかし、もしもRiemann予想が真-非自明な零点がクリティカルライン上にしかでてこない-ならば、素数の分布もRiemann-von Mangoldtの明示公式を通じてある一定の制限を受ける、ということろがとても重要な素数に関するナレッジになるのだ。