ゼータ関数の零点の計算
calculation of the Riemann zeta zeros
ゼータ関数の零点の値を求める。ニュートン法を使用したプログラムにより計算する。
ゼータ関数の級数による式は「級数と乗積 二項係数の級数とゼータ関数」を参照。
sn+1 = sn-ζ(sn)/ζ'(sn) , ζ'(sn)≒(ζ(sn+h)-ζ(sn))/h hは微小な値
より sn+1 = sn-ζ(sn)h/(ζ(sn+h)-ζ(sn)) = sn+h/(1-ζ(sn+h)/ζ(sn)) とする。
h=ζ(sn)/snとして、収束の判定にはζ(s)の絶対値を使用する。UBASICプログラムは次の通り。
220 point -10:word -40:M=1850:dim C(M):C(0)=1:S=1/2+14#i
250 for K=1 to M:for J=K to
1 step -1:C(J)=(C(J)+C(J-1))/2:next:next
300 ' find zeta zero
330 repeat
350 Z=fnZeta(S):H=Z/S:W=fnZeta(S+H):S+=H/(1-W/Z)
360 print
using(8,20),S
380 until abs(Z)<1/10^18
390 end
800 ' zeta function
810 fnZeta(X)
820 local J,U
830 for J=1 to M:U+=(-1)^(J-1)*C(J)/J^X:next:U/=1-2^(1-X)
880 return(U)
非自明な零点の初期値としてs=1/2+14i から開始すると、最初の零点に収束する。
0.50758974993024427495
+14.13558444989598483799#i
0.49997720386661645733 +14.13471625740151952033#i
0.49999999984206884658 +14.13472514154016738957#i
0.50000000000000000001 +14.13472514173469379043#i
0.50000000000000000000 +14.13472514173469379046#i
最初の非自明な零点の正確な値は 0.5+14.13472514173469379045725198…i である。
実数部は初期値 s= 0 +14i から開始しても0.5 になる。
0.41833496948542243845
+14.07645024653248000524#i
0.49916484106395386811
+14.13066995863602495569#i
0.50000693105447366678 +14.13472328584511145747#i
0.49999999997996781442 +14.13472514174264008280#i
0.50000000000000000000 +14.13472514173469379046#i
0.50000000000000000000 +14.13472514173469379046#i
初期値がs= 2 +14i であっても最初の零点に収束する。
-0.92707299992992668900
+14.15826732248992932600#i
-0.12175976448439554067 +14.10335229782430610623#i
0.36268298001830905457 +14.10513082285623101699#i
0.49321098475328483663 +14.13049441833890484983#i
0.49999179184087164570 +14.13469968657916528188#i
0.50000000026583803932 +14.13472514159688702844#i
0.49999999999999999997 +14.13472514173469379048#i
0.50000000000000000000 +14.13472514173469379046#i
より大きな虚数部を持つ零点を求める。s=1/2+1000i を初期値にした例。
0.58685359001856217683 +999.83619290605366949033#i
0.49133736539816212004
+999.76560627179415829529#i
0.50167136925745494820
+999.79060459365507529421#i
0.49999387802951372332
+999.79157899348567507971#i
0.50000000000613823586
+999.79157155765172943357#i
0.50000000000000000015
+999.79157155741294046318#i
0.50000000000000000000
+999.79157155741294046316#i
s=1/2 を初期値にすると、最初の自明な零点である-2 になる。
-2.40497477767311191727
-1.60002505075782479194
-1.88254687558083333765
-1.98685749797395830033
-1.99981381392326113779
-1.99999996200374428046
-1.99999999999999841717
-2.00000000000000000000
-2.00000000000000000000
s=-1+i を初期値にしても同様である。
-1.69567811770189209518
+0.77410534680183550810#i
-2.10707853708190318087 +0.41463987093489378761#i
-2.19247214989578497832 +0.01084314265021423681#i
-1.94571714778927937344 -0.00722609287991353891#i
-1.99702825784503374630 -0.00077784154117950848#i
-1.99999101082137912548 -0.00000503994776057391#i
-1.99999999993925725127 -0.00000000009933926121#i
-2.00000000000000000001 -0.00000000000000000001#i
-2.00000000000000000000
+0.00000000000000000000#i
s=-3 を初期値にすると-4 になる。精度の関係でsをあまり小さな値にすることはできない。
-4.54368465627120233580
-3.66164448434396066935
-3.98241322990050864575
-3.99989288394708271159
-3.99999999588636420686
-3.99999999999999999393
-4.00000000000000000000
ゼータ関数の零点は「四則とべき乗演算」のみで計算可能である。ここまで単純化できる理由は次の通り。「級数と乗積 二項係数の級数とゼータ関数」を参照。
・ζ(1-s)=21-sπ-ssin((1-s)π/2)Γ(s)ζ(s) の関係ではなく、二項係数の級数を使用する。
そのため三角関数やガンマ関数は不要である。
・二項係数をそのまま使用するのではなく、パスカルの三角形の性質を利用する。
そのため二項係数の関数や階乗の関数は不要である。
・ニュートン法においてζ'(s)の値ではなく、その近似値を使用する。
そのためζ(s)の導関数は不要である。
最初の零点を高精度で求める。
220 point -35:word
-140:M=535:dim C(M):C(0)=1:S=1/2+14#i
250
for K=1 to M:for J=K to 1 step -1:C(J)=(C(J)+C(J-1))/2:next:next
300 '
find zeta zero
330
repeat
350 Z=fnZeta(S):H=Z/S:W=fnZeta(S+H):S+=H/(1-W/Z)
360 print using(5,154),S
380
until abs(Z)<1/10^70
390
end
zeta function の部分は最初のプログラムと同じで、実行結果は以下の通り。
0.50758974993024427495199811334098905312233402079728265125362572531457861810
77374401590748788456535754412261057383876701315014301138827009284771818518029915
+14.13558444989598483798911284406977559617481341814039160135843717930301271212
33852167907311487490725608759323729268072110586501578696725322900023173059921159
#i
0.49997720386661645733303694446578370827900867950146832608045760632684583478
85017418880126065947058331809296002799233937576291902727854065813189418181566812
+14.13471625740151952033331883744747939073307413277058415553913333298551823710
40771542434625355197298091885820726164226270277513206788029144585662676265229626
#i
0.49999999984206884658005842136441219790784670793537049857753525266108047381
68259813117330065704648423011733190812583573480510965887637305842999328744308721
+14.13472514154016738956897475703546043662547425037199699887363640635189212688
90837827435143235317533077083699708016302282994998573646851099763328149341700909
#i
0.50000000000000000000907393723220762685587980382709200420812894033319002637
05376922514208378602467257884284147700623783586867978630992205427137365436776310
+14.13472514173469379043258791506095657754225445552059145161748380159678472729
54172259415153087165643780894347674773963818496506651197235153109253850719961904
#i
0.50000000000000000000000000000000000000019065352927224833908483303026830305
57622365918436596030212030043192191740537828042318355061823682977016844512058299
+14.13472514173469379045725198356247027078447444125566453208616403240221837646
00388943062130059334908101700162352977588842445909092160368311259112028992532532
#i
0.50000000000000000000000000000000000000000000000000000000000000000000000000
00000095399963507391630275123917927213918057100300369380206406180391575655328168
+14.13472514173469379045725198356247027078425711569924317568556746014996342980
92567312893943001902628609788235469689087617600795327097308349733328174308711496
#i
0.50000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
+14.13472514173469379045725198356247027078425711569924317568556746014996342980
92567649490103931715610127792029715487974367661426914698822545825053632394471378
#i
最後の桁は四捨五入されているが、小数点以下150桁以上が正しい数字である。
導関数を使用する方法を考える。ただし簡単のためζ(s)の零点を含む
η(s)=(1-21-s)ζ(s) を対象として sn+1 = sn-η(sn)/η'(sn) により求める。
次のプログラムにおいて Uはη(s)であり、Vはη'(s) (正確には級数を項別微分したもの)である。
220
point -10:word -40:M=1850:dim C(M):C(0)=1:S=1/2+14#i
250 for K=1 to
M:for J=K to 1 step -1:C(J)=(C(J)+C(J-1))/2:next:next
300 ' find eta zero
330 repeat
340
U=0:V=0:F=-1
350 for
J=1 to M
360 F=-F:W=F*C(J)/J^S:U+=W:V-=W*log(J)
370
next
380 S-=U/V:print using(8,20),S
390 until
abs(U)<1/10^18
初期値としてs=1/2+14i から開始すると、ゼータ関数の最初の非自明な零点に収束する。
0.51465456894550842195
+14.13390195482083560500#i
0.49982391690453977169 +14.13474103671477568026#i
0.49999997488432152761 +14.13472514570833425528#i
0.49999999999999949602 +14.13472514173469394124#i
0.50000000000000000000 +14.13472514173469379046#i
0.50000000000000000000 +14.13472514173469379046#i
実数部は初期値 s= 0 +14i から開始しても0.5 になる。
0.35368126412084391241
+14.05597738850965127118#i
0.48836413131439907972 +14.11772612182944754364#i
0.50012831065643165891 +14.13440736073379824156#i
0.50000006722951251326 +14.13472520964526432194#i
0.50000000000000024663 +14.13472514173468636525#i
0.50000000000000000000 +14.13472514173469379046#i
0.50000000000000000000 +14.13472514173469379046#i
s=1/2 を初期値にすると、ゼータ関数の最初の自明な零点である-2 になる。
-2.62950644013668291878
-1.63535389888435073795
-1.97712024818981755184
-1.99985425875635826840
-1.99999999388939167555
-1.99999999999999998926
-2.00000000000000000000
-2.00000000000000000000
初期値によっては 1-21-s の零点 1+(2kπ/ln 2)i (kは整数)に収束することがある。
例えば 0+10i の場合は 1+(2π/ln 2)i に収束する。
0.70651995720030240169
+9.69848794457433895642#i
1.08580638698311160453 +9.27078199820368513472#i
1.02120993395715894257 +9.05026681427353198184#i
0.99984006904600630475 +9.06501284739965466332#i
1.00000002382792932377 +9.06472033423899350805#i
1.00000000000000115016 +9.06472028365438654671#i
1.00000000000000000000 +9.06472028365438761926#i
1.00000000000000000000 +9.06472028365438761926#i
上記の高精度版プログラムと実行結果は以下の通り。
220 point -35:word -140:M=535:dim C(M):C(0)=1:S=1/2+14#i
250 for K=1 to M:for J=K to
1 step -1:C(J)=(C(J)+C(J-1))/2:next:next
300 ' find eta zero
330 repeat
340 U=0:V=0:F=-1
350 for J=1 to M
360
F=-F:W=F*C(J)/J^S:U+=W:V-=W*log(J)
370 next
380 S-=U/V:print
using(5,154),S
390 until abs(U)<1/10^70
0.51465456894550842195492089427504225938837617946861742682207079646776555802
54979937092851318295260911884918187092927209053701168447727859499961579416026690
+14.13390195482083560500138423426062771279219603933886812657188384528454639411
44693448984683165044814758657317681513368338968521796432216011967110349039255296
#i
途中の値を6個省略
0.50000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
+14.13472514173469379045725198356247027078425711569924317568556746014996342980
92567649490103931715610127792029715487974367661426914698822545825053632394471378
#i