Maxima で綴る数学の旅

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

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

f:id:jurupapa:20201010190514j:plain

(サバのトマト缶パスタ)

 

このシリーズでは素数個数関数に関するリーマンの明示公式をMaximaで計算して、そのグラフをお見せしたいと思います。

 

\(x\)を正の実数として、\(x\)以下の素数の個数を表す関数を\(\pi(x)\)と書くことにします。例えば\(\pi(10)=4\)です。この素数個数関数\(\pi(x)\)の解析的な計算式(リーマンの明示公式)を与えたのがリーマンです。リーマンの明示公式は近似式ではなく、正確な素数の個数を与えていることに驚かされます。

この辺のお話は、

という本が非常に詳しく、楽しく読めます。

またリーマンの明示公式に関連する話題はtsujimotterさんのブログ「tsujimotterのノートブック」にいくつも記事があり、導出方法については、

に詳しく載っています(ダービーシャーの本では一部省略されています)。

 

この記事で再度、この話題を取り上げた理由は、この公式は複雑でややこしく見えるが、実際に計算できること、その際にMaximaを使うと、式をMaximaでそのまま入力し、少しの数学を使うことでほぼそのまま計算できること、などを知っていただきたいからです。

では早速始めましょう。まず(%i1)で、数値計算を中心に行うこと、メビウス関数を綺麗に表示すること、\(\rho\)という変数が複素数であることをMaximaに知らせます。

(%i1) keepfloat:true$
assume(x>1)$
texput(moebius,"\\mu")$
declare(rho,complex)$

いよいよ、リーマンの明示公式とのご対面です。リーマンはまず\(J(x)\)という関数を定義し、次にそれを使って\(\pi(x)\)を定義しました。

(%i4) J(x)=li(x)-sum(li(x^rho[i])+li(x^(conjugate(rho[i]))),i,1,inf)-log(2)+integrate(1/(t*(t^2-1)*log(t)),t,x,inf);

$$ \tag{${\it \%o}_{4}$}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 $$
(%i5) pi(x)=sum(moebius(m)/m*J(x^(1/m)),m,1,floor(log(x)/log(2)));

$$ \tag{${\it \%o}_{5}$}\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}} $$

この2つの式で定義される\(\pi(x)\)を計算すると、\(x\)以下の素数の個数が計算される、というのですから驚きです。式自体、とても複雑です。ちょっと説明しますが、まあ分からなくても良いです。

(%o4)を見てください。\(J(x)\)の第2項\(li(x)\)は対数積分関数と呼ばれる関数です。Maximaでは\(expintegral\_li(x)\)として実装されています。第3項の積分は複雑そうに見えますが、\(x\)の値が決まれば数値積分で近似値を得ることができます。第4項\(log 2\)は2の自然対数で、大体\(0.6931\)くらいです。第1項は、、、かなり複雑なのでこのシリーズの第2回で説明することにします。ただし1つだけ覚えておいてください。第1項に現れている\(\rho_i\)はリーマンのゼータ関数の自明でない零点(無限個あります)を表しているのです。

 (%o5)に\(\pi(x)\)の明示公式の定義があります。この中で\(\mu(x)\)という関数はメビウス関数と呼ばれる有名な関数です。\(J(x)\)は(%o4)で定義されたものです。メビウス関数と\(J(x)\)を掛け算してから総和をとっていることがわかります。

 

ところで、素数個数関数のグラフをかけ、と言われたら、普通のプログラマは以下のように実装すると思います。要は1から順番に素数の個数を数えていくわけです。ただしちょうど関数がジャンプする(不連続になる)素数のところではジャンプ前と後の中間の値を取るように調整します。Maximaで書くとこんな感じです。

(%i6) pi1[n]:=if n<2 then 0 elseif primep(n) then pi1[n-1]+1 else pi1[n-1]$

(%i7) pi(x):=if integerp(x) and primep(x) then pi1[x]-1/2 else pi1[floor(x)]$

この\(\pi(x)\)をグラフに書くとこうなります。

(%i8) plot2d(pi(x),[x,0,100])$

f:id:jurupapa:20201010154813p:plain

素数の個数を数えているのですから当然なのですが、xが素数の時に1ジャンプし、それ以外のところでは平です。

ここで、この数え上げによる定義とリーマンの明示公式が大局的にどのくらい一致するのかを確認しておきます。先ほどの\(J(x)\)のややこしい第1項を除いた関数を\(J2(x)\)として定義します。対数積分関数をMaximaで実装されている\(expintegral\_li(x)\)に置き換えます。最後に積分を数値積分コマンドromberg()に置き換え、積分範囲を無限から10に置き換えます(ここは被積分関数のグラフの概形から10で十分なのです)。

(%i9) J2(x):=expintegral_li(x)-log(2)+romberg(1/(t*(t^2-1)*log(t)),t,x,10);

$$ \tag{${\it \%o}_{9}$}J_{2}\left(x\right):={\it expintegral\_li}\left(x\right)-\log 2+{\it romberg}\left(\frac{1}{t\,\left(t^2-1\right)\,\log t} , t , x , 10\right) $$

\(\pi_2(x)\)は\(J2(x)\)を呼び出すように変更しただけです。

(%i10) pi2(x):=ev(sum(moebius(m)/m*J2(x^(1/m)),m,1,floor(log(x)/log(2))),numer);

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

素数を数え上げる関数\(\pi(x)\)と解析的に定義した\(\pi_2(x)\)を同じグラフにプロットしてみます。

(%i??) plot2d([pi(x),pi2(x)],[x,2,1000])$

f:id:jurupapa:20201010154853p:plain

ガタガタしている青い線が\(\pi(x)\)なのですが、その上を、滑らかに赤い\(\pi_2(x)\)が重なっていることが分かります。素数定理をご存知の方なら、この\(\pi_2(x)\)の方が遥かに正確に素数個数関数を大局的に近似していることがわかると思います。

\(\pi_2(x)\)は\(J(x)\)関数の第1項を除いて計算したものです。別の言い方をすると、\(J(x)\)関数の第1項を計算に加えると、滑らかな\(\pi_2(x)\)が、\(x\)が素数のところでジャンプするようになるのです。

もう少しいえば、滑らかな\(\pi_2(x)\)の「どこに素数があるのか」を\(J(x)\)の第1項は正確に知っているのです。そしてその第1項こそが、ゼータ関数の零点からの寄与なのです。