今回はマンゴールト関数の総和関数をフーリエ変換(正確にはメリン変換)すると周波数成分としてゼータ関数の非自明な零点が現れることを見てみます。
この記事を書くにあたってThe Mobius function Blogの記事Fourier transform of the von Mangoldt function with first term equal to a harmonic number及びそこに登場するMathematicaのコードを参考にさせて頂きました。このブログはそのタイトルからして渋いのですが、中身も面白いです。しかしその話はまた別の機会に。
まずマンゴールト関数を定義し(%o1〜%o3)、%o4でメリン変換の定義を確認。
(%i1) mangoldt_c(N):=block([fct],fct:ev(ifactors(N),factors_only:true),
if length(fct)=1 then log(fct[1]) else 0)$
(%i2) (matchdeclare(n,integerp),tellsimpafter(mangoldt(n),mangoldt_c(n)))$
(%i3) texput(nounify(mangoldt), "\\Lambda")$
(%i4) h(t)=integrate(x^(-1/2+%i*t)*f(x),x,0,inf);
$$ \tag{%o4} h\left(t\right)=\int_{0}^{\infty }{x^{i\,t-\frac{1}{2}}\,f\left(x\right)\;dx} $$
この一般的な定義にマンゴールト関数の総和関数を代入した物がこちらです。
(%i5) integrate(x^(-1/2+%i*t)*f(x),x,0,inf), f(x):=sum(mangoldt(n),n,1,floor(x))-x;
$$ \tag{%o5} \int_{0}^{\infty }{\left(\sum_{n=1}^{\left \lfloor x \right \rfloor}{\Lambda\left(n\right)}-x\right)\,x^{i\,t-\frac{1}{2}}\;dx} $$
すこし計算を進めて、計算が高速になるように準備します。
(%i6) assume(x>0);
$$ \tag{%o6} \left[ x>0 \right] $$
(%i7) rectform(x^(-1/2+%i*t));
$$ \tag{%o7} \frac{i\,\sin \left(t\,\log x\right)}{\sqrt{x}}+\frac{\cos \left(t\,\log x\right)}{\sqrt{x}} $$
(%i8) define(img(x,t),imagpart(%));
$$ \tag{%o8} \mathrm{img}\left(x , t\right):=\frac{\sin \left(t\,\log x\right)}{\sqrt{x}} $$
上記の関数img(x,t)は後の計算で使います。
ここから、上記のメリン変換を離散化した物を計算するプログラムとなります。先ほど紹介したブログ記事のMathematicaのプログラムをMaximaで動くように書き換えた物です。
(%i9) kill(f)$
(%i10) scale:5000$
(%i11) f:make_array(flonum,scale+1)$
(%i12) f[1]:0$
まずscaleまでの範囲でマンゴールト関数の総和関数を計算してしまいます。
(%i13) for i:2 thru scale do f[i]:float(mangoldt(i))+f[i-1]$
以下、xを対数軸で考えるため、ちょっと特殊な準備をしています。
(%i14) xres:.002$
(%i15) xmax:floor(log(scale)/xres);
$$ \tag{%o15} 4258 $$
(%i16) xlist:make_array(flonum,xmax)$
(%i17) for i:0 thru xmax-1 do xlist[i]:exp(xres*i)$
(%i18) fx:make_array(flonum,xmax)$
対数軸に沿ってマンゴールト関数の総和関数を評価して、fx[i]という配列に格納します。
(%i19) for i:0 thru xmax-1 do fx[i]:f[floor(xlist[i])]-xlist[i]$
以下配列の大きさを取り出す関数、2つの配列を次元の大きなベクトルと思って内積をとる関数を定義します。
(%i20) array_size(a1):=subvar(subvar(arrayinfo(a1),3),1)$
(%i21) inner_product(a1,a2):=block([len,total], /* a1 and a2 must be arrays */
len:min(array_size(a1),array_size(a2)),
total:0,
for i:0 thru len-1 do total:total+a1[i]*a2[i],
total)$
(%i22) tmax:100$
(%i23) tres:0.14$
関数f1()が実際に離散的にメリン変換を行う関数です。img(x,t)を計算し、fx[i]との内積をとります。それをどんどん繰り返します。
(%i24) f1():=block([maxt,t,xexp,errlist],
errlist:,
maxt:floor(tmax/tres),
xexp:make_array(flonum,xmax),
for i:0 thru maxt do (
if mod(i,100)=0 then print(i,maxt),
t:i*tres,
for j:0 thru xmax-1 do xexp[j]:float(img(xlist[j],t)),
errlist:append(errlist,t,inner_product(xexp,fx))),
errlist)$
以下、定義した関数をコンパイルします。
(%i25) compile(array_size,inner_product,f1,imf)$
離散化したメリン変換の結果はerrlistに入ります。下のグラフを見ると所々で下に大きくひげが伸びています。このひげの位置のx座標がまさにゼータ関数の非自明な零点と一致しているのです。
(%i26) errlist:f1()$
%i27〜%i29でリストzzにゼータ関数の非自明な零点を入れています。
(%i27) load("zetazeros1000.mac")$
(%i28) zz:$
(%i29) for i:1 while img_rho[i]<tmax do zz:append(zz,[[img_rho[i],1]])$
以下、zzとerrlistを同じグラフに書いてみます。ひげの位置とzzの位置(赤い三角)が正確に対応していることが分かります。
(%i30) plot2d([[discrete,errlist],[discrete,zz]],[y,-400,100],[style,line,points])$
ちなみに原点に近い方から3つの赤い三角のx座標は以下の通り、ゼータ関数の非自明な零点の最初の3つです。
(%i31) [img_rho[1],img_rho[2],img_rho[3]];
$$ \tag{%o31} \left[ 14.13472514173469 , 21.02203963877156 , \\ 25.01085758014569 \right] $$
まとめると、その1ではゼータ関数の非自明な零点から素数の位置(マンゴールト関数)が分かることを見ました。今回であるその2では素数の位置(マンゴールト関数の総和関数)からゼータ関数の非自明な零点の位置が分かることを見ました。