Maxima で綴る数学の旅

紙と鉛筆の代わりに、数式処理システムMaxima / Macsyma を使って、数学を楽しみましょう

-数学- リーマンの素数個数関数の明示公式 (2)

f:id:jurupapa:20201012211901j:plain

 

この記事のコードをJupyter notebook形式で見るにはこちら。

 

素数個数関数に対するリーマンの明示公式を再掲します。

$$ J\left(x\right)=-\sum_{i=1}^{\infty }{\left({\it li}\left(x^{\rho_{i}^\star}\right)+{\it li}\left(x^{\rho_{i}}\right)\right)}+{\it li}\left(x\right)+\int_{x}^{\infty }{\frac{1}{t\,\left(t^2-1\right)\,\log t}\;dt}-\log 2 $$
$$ \pi\left(x\right)=\sum_{m=1}^{\left \lfloor \frac{\log x}{\log 2} \right \rfloor}{\frac{\mu\left(m\right)\,J\left(x^{\frac{1}{m}}\right)}{m}} $$

\(J(x)\)の第1項\(\sum_{i=1}^{\infty }{\left({\it li}\left(x^{\rho_{i}^\star}\right)+{\it li}\left(x^{\rho_{i}}\right)\right)}\)は複雑で難しいので、前回の記事ではこれを除いた項だけ計算してみました。

素数を数え上げる形で計算する\(\pi(x)\)と\(J(x)\)の第1項を除いて定義した\(\pi_2(x)\)をグラフに描くと、後者は滑らかなグラフですが、前者とよく重なることが分かります。別の言葉で言えば、 \(J(x)\)の第1項\(\sum_{i=1}^{\infty }{\left({\it li}\left(x^{\rho_{i}^\star}\right)+{\it li}\left(x^{\rho_{i}}\right)\right)}\)は、\(\pi_2(x)\)上のどこに素数があるのかを完全に知っているのです。

 

今回は、この\(J(x)\)の第1項をMaximaで計算してみます。前回の結果と合わせれば、\(J(x)\)全体、及び明示公式そのものをグラフに描くことができます。素数を数え上げる形で計算したガタガタの\(\pi(x)\)とガタガタを含めて驚くほどよく重なることが分かります。それを目指して、頑張りましょう。

 

\(J(x)\)の第1項の計算の前に、\(\sum_{i=1}^{\infty }{\left({\it li}\left(x^{\rho_{i}^\star}\right)+{\it li}\left(x^{\rho_{i}}\right)\right)}\)について少し解説をします。まずこの項は\(i\)を\(1\)から∞まで変化させながら総和をとっています。総和の中身のどこに\(i\)が出てくるのか、と言えば、\(x\)の指数のところにある\(\rho_i\)と\(\rho_{i}^\star\)のところです。\(\rho_i\)が何かというとこれがあの有名な「リーマンのゼータ関数の自明でない零点」、ということになります。(リーマン予想を仮定すれば)\(\rho_i\)はすべて複素数で\(\frac{1}{2}+{\it it}\)という形をしています。

\(li(x^{\rho_i})\)には対数積分関数が登場していますが、Maximaで計算する際にここにちょっと問題があり、簡単な数学を使って書き換える必要があります。\(x\)は実数ですが、\(\rho_i\)は複素数ですから、\(x^{\rho_i}\)は複素数です。ところがMaximaの対数積分関数\(expintegral\_li(z)\)は引数が複素数の時、限定された範囲でしか計算ができないのです。(対数関数が多価関数であることに関係しています)。

この分野には指数積分関数\(Ei(z)\)というものがあり、\(Ei(log(z))=li(z)\)という関係がすべての複素数\(z\)について成り立ちます。そしてこの両辺に零点の形\(z=\frac{1}{2}+{\it it}\)を代入すると(%o12)を得ることができます。

(%i12) expintegral_ei( (1/2+%i*t)*log(x) )=expintegral_li(x^(1/2+%i*t));

$$ \tag{${\it \%o}_{12}$}{\it expintegral\_ei}\left(\left(i\,t+\frac{1}{2}\right)\,\log x\right)={\it expintegral\_li}\left(x^{i\,t+\frac{1}{2}}\right) $$

この右辺(赤い線)と左辺(青い線)を複素平面上にプロットしてみると右辺の問題が浮き彫りになります。

(%i13) draw2d(nticks=1000,
parametric( realpart(expintegral_ei(expand((1/2+%i*t)*log(20.0)))),
imagpart(expintegral_ei(expand((1/2+%i*t)*log(20.0)))),t,-50,50) ,
color=red,
parametric( realpart(expintegral_li(rectform(20.0^(1/2+%i*t)))),
imagpart(expintegral_li(rectform(20.0^(1/2+%i*t)))),t,0,2.1));

f:id:jurupapa:20201010155010p:plain

青い線は線\(\frac{1}{2}+{\it it}\)が無限に伸びるのを、渦巻きに巻きつけています。これが本来の\(li(x^{\rho_i})\)のあるべき姿なのです。一方そのまま組込のexpintegral_li()を呼び出して描いた赤い線は、青い線と重なっている部分(ここは正確です)はわずかで、その範囲を超えると値が元に戻ってしまうことが分かります。これではうまくありません。(%o12)の左辺を使うことで、対数積分関数を複素数の範囲でも計算できるようになります。

第1項の計算の最後の工夫は、\(x^{\rho_{i}^{\star}}\)の処理です。ちなみに\(\rho_{i}^{\star}\)は\(\rho_{i}\)の複素共役を表します。グラフでわかるように、\(li(x^{\frac{1}{2}+{\it it}})\)が\(t\)の関数として複素平面上で実軸に対象であることなどを用いると、\(li(x^{\rho_i})\)と\(li(x^{\rho_i^{\star}})\)は複素共役であることが分かり、第1項は次の式で計算できることが分かります。

(%i15) Li_power(x,t):=2*realpart(expintegral_ei(expand( (1/2+%i*t)*log(x) ) ) );

$$ \tag{${\it \%o}_{15}$}{\it Li\_power}\left(x , t\right):=2\,{\it realpart}\left({\it expintegral\_ei}\left({\it expand}\left(\left(\frac{1}{2}+i\,t\right)\,\log x\right)\right)\right) $$

さて、ゼータ関数の非自明零点の計算なのですが、これは他の方がやった計算結果をダウンロードすることにしました。ここでは数論データベースLMFDB [1] から最初の200個を使わせて頂いています。
(%i15) /*
The LMFDB Collaboration, The L-functions and Modular Forms Database,
home page of the Zeros of zeta(s),
https://www.lmfdb.org/zeros/zeta/?limit=200&N=1, 2020 , [Online; accessed 10 October 2020].
*/
img_rho:[
14.1347251417346937904572519835625,
21.0220396387715549926284795938969,
25.0108575801456887632137909925628,

途中省略。合計200この零点をダウンロードしました。この計算のipynbファイルはすべての零点を含んでいます。

392.2450833395190967490151841709930,
393.4277438444340259366989529201288,
395.5828700109937209708777113231417,
396.3818542225921869319994544917305
]$

 

これですべての準備が整いました。以下に\(J3(x)\)として第1項の計算に先ほど定義した\(Li\_power(x,t)\)を使って定義します。また\(J3(x)\)を呼び出す\(\pi_3(x)\)も定義します。

(%i16) J3(x):=expintegral_li(x)-sum(2*realpart(Li_power(x,img_rho[i])),i,1,200)-log(2.0)+romberg(1/(t*(t^2-1)*log(t)),t,x,20);

$$ \tag{${\it \%o}_{16}$}J_{3}\left(x\right):={\it expintegral\_li}\left(x\right)-{\it sum}\left(2\,{\it realpart}\left({\it Li\_power}\left(x , {\it img\_rho}_{i}\right)\right) , i , 1 , 200\right)-\log 2.0+{\it romberg}\left(\frac{1}{t\,\left(t^2-1\right)\,\log t} , t , x , 20\right) $$
(%i17) pi3(x):=ev(sum(moebius(m)/m*J3(x^(1/m)),m,1,floor(log(x)/log(2))),numer);

$$ \tag{${\it \%o}_{17}$}\pi_{3}\left(x\right):={\it ev}\left({\it sum}\left(\frac{\mu\left(m\right)}{m}\,J_{3}\left(x^{\frac{1}{m}}\right) , m , 1 , \left \lfloor \frac{\log x}{\log 2} \right \rfloor\right) , {\it numer}\right) $$

 

素数を数え上げる形で定義された\(\pi(x)\)と上記で定義されたリーマンの明示公式通りの関数\(\pi_3(x)\)をグラフに描画してみます。
(%i19) plot2d([pi3(x),pi(x)],[x,3,200]);

f:id:jurupapa:20201010155406p:plain

見事に重なっています。もちろんよく見ると\(x=200\)に近づくと重なりが甘くなっていることが分かります。もっと多数の非自明零点を使えばこの重なりはどんどん一致するはずです。

 

[1] The LMFDB Collaboration, The L-functions and Modular Forms Database,
home page of the Zeros of zeta(s),
https://www.lmfdb.org/zeros/zeta/?limit=200&N=1, 2020 , [Online; accessed 10 October 2020].

 

今回もこちらの本を参考にしました 。特に第21章「誤差項」です。