Maxima で綴る数学の旅

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

-数学- マンゴールトの明示公式

 

f:id:jurupapa:20130224013415j:plain

 

マンゴールト関数の総和関数の明示公式(マンゴールトの明示公式)のお話です。次の本の第3章にマンゴールトの明示公式の導出方法などが載っています。

明解 ゼータ関数とリーマン予想 (KS理工学専門書)

明解 ゼータ関数とリーマン予想 (KS理工学専門書)

前回マンゴールト関数とその総和関数を定義して、そのグラフを書いてみました。

(%i1) mangoldt(N):=if integerp(N)

  then block([fct],fct:ev(ifactors(N),factors_only:true),

                   if length(fct)=1 then log(fct[1]) else 0)

  else 'mangoldt(N);

$$ \tag{%o1} \mathrm{mangoldt}\left(N\right):=\mathbf{if}\;\mathrm{integerp}\left(N\right)\;\mathbf{then}\;\mathbf{block}\;\left(\left[ \mathrm{fct} \right]  , \\ \mathrm{fct}:\mathrm{ev}\left(\mathrm{ifactors}\left(N\right) , \mathrm{factors\underline{\quad}only}:\mathbf{true}\right) , \\ \mathbf{if}\;\mathrm{length}\left(\mathrm{fct}\right)=1\;\mathbf{then}\;\log \mathrm{fct}_{1}\;\mathbf{else}\;0\right) \\ \;\mathbf{else}\;\Lambda\left(N\right) $$

(%i2) psi0(x):=sum(mangoldt(n),n,1,floor(x));

$$ \tag{%o2} \mathrm{psi0}\left(x\right):=\mathrm{sum}\left(\mathrm{mangoldt}\left(n\right) , n , 1 , \left \lfloor x \right \rfloor\right) $$

その際にマンゴールト関数の総和関数は、素数個数関数 \(\pi(x)\) と同じような素数に関するステップ関数であり、研究ではよく使われる、と記しました。

実際、この関数もゼータの零点を使って明示公式を与えることができます。この公式はマンゴールトの明示公式と呼ばれています。

(%i3) load("zetazeros1000.mac")$

以下がそのマンゴールトの明示公式です。リーマンの素数個数明示公式よりもはるかに簡単な姿をしています。

(%i4) psi1(x,k)=x-log(2*%pi)-1/2*log(1-x^(-2))-sum(1/rho[i]*x^rho[i]+1/conjugate(rho[i])*x^conjugate(rho[i]),i,1,k);

$$ \tag{%o4} \mathrm{psi1}\left(x , k\right)=-\sum_{i=1}^{k}{\left(\frac{x^{\rho_{i}^\star}}{\rho_{i}^\star}+\frac{x^{\rho_{i}}}{\rho_{i}}\right)}+x-\frac{\log \left(1-\frac{1}{x^2}\right)}{2} \\ -\log \left(2\,\pi\right) $$

この関数、このままだと計算に長い時間がかかります。そこで式変形を行って高速に計算出来る形にしてみます。問題はゼータの零点を使っているsum()の部分です。以下、式変形のための準備です。

(%i5) (matchdeclare(agr1,true),tellsimpafter(atan2(0,arg1),0))$

(%i6) declare(rho,complex,x,real,a,real)$

(%i7) assume(x>0);

$$ \tag{%o7} \left[ x>0 \right]  $$

では早速始めましょう。

(%i8) x^rho/rho+x^conjugate(rho)/conjugate(rho);

$$ \tag{%o8} \frac{x^{\rho^\star}}{\rho^\star}+\frac{x^{\rho}}{\rho} $$

ここで使うゼータの零点は最初の数百個で、それらの実部は1/2であることが証明されています。なので、以下のような代入を行います。

(%i9) %,rho:1/2+%i*a;

$$ \tag{%o9} \frac{x^{i\,a+\frac{1}{2}}}{i\,a+\frac{1}{2}}+\frac{x^{\frac{1}{2}-i\,a}}{\frac{1}{2}-i\,a} $$

簡単にしてみます。

(%i10) %,rectform,ratsimp;

$$ \tag{%o10} \frac{\sqrt{x}\,\left(8\,a\,\sin \left(a\,\log x\right)+4\,\cos \left(a\,\log x\right)\right)}{4\,a^2+1} $$

得られた式に関数名を付けます。

(%i11) define(periodic(x,a),%);

$$ \tag{%o11} \mathrm{periodic}\left(x , a\right):=\frac{\sqrt{x}\,\left(8\,a\,\sin \left(a\,\log x\right)+4\,\cos \left(a\,\log x\right)\right)}{4\,a^2+1} $$

定義した\( \mathrm{periodic}\left(x , a\right) \) を使ってマンゴールトの明示公式を書き直します。

(%i12) psi2(x,k):=float(x-log(2*%pi)-1/2*log(1-x^(-2))-sum(periodic(x,img_rho[i]),i,1,k));

$$ \tag{%o12} \mathrm{psi2}\left(x , k\right):=\mathrm{float}\left(x-\log \left(2\,\pi\right)+\frac{-1}{2}\,\log \left(1-x^ {- 2 }\right) \\ -\mathrm{sum}\left(\mathrm{periodic}\left(x , \mathrm{img\underline{\quad}rho}_{i}\right) , i , 1 , k\right)\right) $$

コンパイルします。

(%i13) compile(periodic,psi2)$

(%i14) wxplot2d(['psi0(x),'psi2(x,50)],[x,2,100]);

$$ \tag{%o14}  $$

f:id:jurupapa:20130224013338j:plain

ゼータの零点を50個利用した赤い線は階段の上を綺麗に走っています。

(%i15) wxplot2d(['psi0(x),'psi2(x,200)],[x,2,100]);

$$ \tag{%o15}  $$

 

f:id:jurupapa:20130224013415j:plain

 

ゼータの零点を200個利用してみると、この範囲では誤差がほとんど無く、赤い線(明示公式)が階段に完全に重なっていることが分かります。