少し前に、素因数分解に関するプログラムを作っている時に、素数のことを調べてみたら、リーマン予想という未解決？と言われる問題があることを知りました。

リーマン予想：リーマンゼータ関数の、自明でない零点は、全てが実数値が1/2の値を取る。

ということで、今回は、プログラムを使って、リーマンゼータ関数の零点を、プログラムで計算してみたいと思います。

リーマンゼータ関数というものは、平たく言えば、代入値の範囲を1点を除く複素数全体まで拡張した、特殊級数の和の形を取る関数です。

$${ζ(s)=\sum\limits_{n=1}^{\infty} \dfrac{1}{n^s}}$$

元の関数の形は、この形です。この形の複素数型の関数のうち、絶対に収束が見られるsの実数値が1より大きい範囲で定義されたものをゼータ関数といいます。
ベルンハルト・リーマンはこの関数にゼータ関数と名付け、解析接続という特殊な変換を用いて、定義域を複素数全体まで拡張したそうです。

プログラムの実装にあたり、上のページを参照しました。

記述によれば、ゼータ関数は、sの実部が1よりも大きければ

$${ (1-\dfrac{1}{2^s})ζ(s)=1-\dfrac{1}{2^s}+\dfrac{1}{3^s}-…}$$

を満たすそうです。一方、この級数はより広くsの実部が正の時にいつでも収束するそうです。リーマンゼータ関数は、解析接続というテクニックを使って、複素数全体を定義域とする関数に拡張されます。拡張された後のリーマンゼータ関数は、実数部が正の複素数のうち、$${s=1}$$を除く全ての複素数値で収束します。今回はこの式をベースに使い、また、シキノートさんのページの記述にある漸化式表示を参照させて頂き、リーマンゼータ関数を表現していこうと思います。

漸化式

$${ζ(s)=\sum\limits_{n=0}^{\infty} b_n}$$

$${b_n=\sum\limits_{k=0}^n a_k^n}$$

$${a_k^n=\dfrac{1}{1-2^s} \dfrac{1}{2^{n+1}} \dbinom{n}{k} \dfrac{(-1)^k}{(k+1)^s}}$$

プログラムの方を上げます。

import numpy as np
import math
import cmath
s=-1
r=0.1
vec=np.zeros((20,20),dtype=complex)
sdash=1
def enbyouga(s,r):
    slist=[0]*51
    tatebai=1
    for i in range(51):
        slist[i]=s+r*np.cos(2*np.pi/50*i)+tatebai*r*1j*np.sin(2*np.pi/50*i)
    return slist
print(1/2**(2.5+1j),1**s)
print(1-1/2*2.5+1/4*2.5*1.5/2-1/8*2.5*1.5*0.5/6-2.5*1.5*0.5*0.5/24/16)
slist=enbyouga(s,r)
for i in range(20):

    for j in range(1,21):
        vec[i,j-1]=sdash*(j+1)**(-s-i)/math.factorial(i)
    #if s==-i:
        #continue
    sdash*=-s-i
riemann=0

def riemann0(slist):
    kousuu=400
    zenka3=[0]*kousuu
    zenka2=[0]*kousuu
    riemann=0
    preriemann=0
    henkaku=0
    for slis in slist:
        zenka1=1/2/(1-2**(1-slis))
        zenka3[0]=zenka1
        riemann+=zenka3[0]
        #zenka2[0]=zenka1
        #print(zenka3,zenka2)
        for i in range(1,kousuu):
            for j in range(i):
                zenka3[j]=i/(i-j)/2*zenka3[j]
                riemann+=zenka3[j]
            zenka3[i]=-zenka3[i-1]/i*(i/(i+1))**slis
            riemann+=zenka3[i]
            #print(zenka3,zenka2)
        henkaku+=(riemann-preriemann)/riemann
        #print(slis,riemann)
        #print(riemann,'riemann')
        preriemann=riemann
        riemann=0
    return henkaku
    #print(henkaku/np.pi/2,'henkaku')
superr=2
supers=0.5+20j
minis=0
skouho=[0]*4
tatebai=1
print(riemann0([2]))
henkaku=0

slist=enbyouga(supers,superr)
for i in range(20):
    henkaku=riemann0(slist)
    #print(henkaku)
    if henkaku.imag/2/np.pi>0.8:
        minis=supers
        skouho=[supers+(0.5+tatebai*0.5j)*superr,supers+(0.5-tatebai*0.5j)*superr,supers+(-0.5-tatebai*0.5j)*superr,supers+(-0.5+tatebai*0.5j)*superr]
        for s in skouho:
            slist=enbyouga(s,superr*0.7)
            henkaku=riemann0(slist)
            if henkaku.imag/2/np.pi>0.8:
                supers=s
                superr*=0.7
                break
    
    else:
        break
    print(henkaku/2/np.pi,supers,superr)
print(henkaku/2/np.pi,supers,superr)
#print(lists)
        
        

細かい関数を複素数での計算用にcmathに切り替えた方が良かったかもしれませんが、とりあえず細かい調整はしていません。真ん中の関数riemann0がリーマンゼータ関数値の計算をする関数。無限和を計算することまではさすがにできませんので、kousuu変数で指定したnが400項くらいまでの総和を計算します。上がコーシーの偏角の原理を利用して、零点の探索を行う関数です。
superr変数を半径、supersを中心座標としてプログラムを立ち上げると、定義した範囲内に零点が有れば反応して、零点の座標をゆっくりと反復計算にて推測します。
この条件では、探索の結果座標が(0.49989
663250330824+21.022533746011124j) となり、推測する半径の広さが0.00159584532595
2238まで狭めることができます。
コーシーの偏角を使った座標の推定は、以下のサイトを参考にしています。

今回はただ無限項まで計算して合ってるかどうか見るだけでは面白くなかったので、結果が再現されるnの下限を探そうと思い立ち、計算することにしてみました。kousuu変数を自然数の下限まで動かして、零点のうち原点から2番目に近い、s=0.5+21.022iの値の変化を調べていきます。

kousuu=16
　s=0.38747615 + 20.924788 i
kousuu=15
　s=0.3331979 + 20.903045 i
kousuu=14
　s=0.2596477 + 20.881023 i
kousuu=13
　s=0.1613366 + 20.85745 i
kousuu=12
　s=0.030859 + 20.833314 i
kousuu=11
　s=-0.1422114 + 20.808947 i
kousuu=10
　s=-0.373407 + 20.784952 i
kousuu=9
　s=-0.687127 + 20.7635655 i
kousuu=8
　s=-1.1227422 + 20.7532 i
kousuu=7
　s=-1.741365 + 20.78013 i
kousuu=6
　s=-2.6316786 + 20.9036202 i
kousuu=5
　s=-3.9488732 + 21.2368 i
kousuu=4
　s=-6.13617128 + 22.0935 i

ざっとこんな感じです。kousuuが2以下の時の零点は見失いました。見つからなかったです。わかったことは、nの項数に関わらず、零点の場所と数に対して、大きな変化は見られないこと。nの上限が小さくなればなるほど、一定の規則性を持って、リーマンゼータ関数の零点とされる位置からは、実数値虚数値共に、少しずつ離れていっていること。nの上限が小さい場合の零点の実数値は0よりもちいさくなることも、可能性としてはあり得ることです。

kousuu=3(n=0,1,2)の場合の零点を手計算で求めてみる

一昨日にTwitterで書いた計算です。

$${1+\dfrac{1}{2^s} +\dfrac{1}{3^s}=0}$$となる未知数sを手計算で求めるとしたらどうすれば求まるか？

①それぞれの分数を自然対数eを底とした指数表示に直す。②$${e^{-s}}$$をまとめてひとつの文字に置く。③近似した答えを求めるならば、log2とlog3を実際の数値に置き換えて、指数を整数値の中で一番比が近い数値で表せば、多項式が出来上がり、近似解を求めることができる。

①$${1+e^{-slog2}+e^{-slog3}=0}$$
②$${x=e^{-s}}$$とすると、$${1+x^{log2}+x^{log3}=0}$$
③log2=0.3010,log3=0.4771なので、log2：log3=2：3の自然数となるように$${y=x^0.15}$$として、$${1+y^2+y^3=0}$$を解く。

この方程式の答えは、y=-1.4655,0.2327±0.7955 iとなる。
この解のうちの、実数解は、$${e^{-0.15s}$$=-1.4655を求めることで、計算することができて、値はオイラーの公式を使い、s=-2.5479+20.9439i ± 41.8799 i を得ることができます。

この数値は、リーマンゼータ関数の零点の数値からは離れているものの、先程のプログラムのkousuuを動かして行った時の連続性をある程度推定できる数値となっております。

自明でない零点の実数値に0.5しか見つからない理由を、推論計算で推論してみる(前半)

リーマンゼータ関数の零点の位置には、ある程度の規則性がありそうなことが予想されるので、この規則性を探っていきます。
まず$${\dfrac{1}{n^s}}$$についてです。
eを底とする指数関数の形に直すと、$${e^{-s logn}}$$となります。lognをxとしてまとめて見てみると、複素数平面に映したこの式のグラフは対数螺旋(極座標表示で$${r(θ)=ae^{bθ}}$$)の形となっています。

直交座標における媒介変数表示は

$${x(θ)=ae^{bθ}cosθ、y(θ)=ae^{bθ}sinθ}$$
だそうで、元の式を

$${\dfrac{1}{e^{sx}} (s=a+bi:a,bは実数,x=logn)}$$
と見ると、オイラーの公式により

$${e^{-sx}=e^{(-a-bi)x}=e^{-ax}(cosbx-isinbx)}$$
と表示できます。
実数軸が$${e^{-ax}cos(-bx)}$$、虚数軸が$${e^{-ax}sin(-bx)}$$の表示となるので、xを介して対数螺旋の表示になっていることがわかると思います。

対数螺旋の性質としては、自己相似というものがあります。すなわち、任意の倍率で拡大したものは、適当な回転によって元の螺旋と一致します。特に、元の関数を$${r(θ)=ae^{bθ}}$$とした時の$${e^{2πb}}$$倍に拡大したものは、回転することなしに元の螺旋と一致することがわかっています。

また、中心から伸ばした半直線と螺旋は無限回交わりますが、この半直線との交点それぞれについて、隣り合う交点との原点との距離の比は一定で$${e^{2πb}}$$になります。

$${n^{-s}}$$について$${x=logn}$$を介した対数螺旋のグラフを描くと紹介しましたが、全ての$${n^{-s}=e^{(-a-bi)x}}$$上の点は、x=0の時の複素数1からの、道のり$${x=logn}$$毎にプロットした座標上を通ることになります。

さて、リーマンゼータ関数を
$${ζ(s)=\sum\limits_{n=1}^{\infty} \dfrac{1}{n^s}}$$
として紹介しましたが、自明でない零点は解析接続されたリーマンゼータ関数で見つかる値なので、今後は解析接続された後の
$${ (1-\dfrac{1}{2^s})ζ(s)=1-\dfrac{1}{2^s}+\dfrac{1}{3^s}-…}$$
の右辺の0への収束について考えていきます。

この式の右辺は次のように変形することができます。
$${1-\dfrac{1}{2^s}+\dfrac{1}{3^s}-…=\sum\limits_{n=1}^\infty \dfrac{1}{(2n+1)^s} -\sum\limits_{n=1}^\infty \dfrac{1}{(2n)^s}   (①)}$$
なので、リーマンゼータ関数の収束を考えるためには、この式の値がゼロになるかどうかを考えればいいのですが、$${n^{-s}}$$の値があるひとつの対数螺旋の軌道上に存在することから、あるひとつの解釈が成り立ちます。

ある自己相似たる対数螺旋$${e^{-sx}}$$の、拡大・縮小に対して対称な座標群の座標の総和がゼロになるとき、この座標群の要素の総和で構成されるリーマンゼータ関数はゼロに収束する。

もし、これが成立するならば、対数螺旋上の収束を調べることで、リーマンゼータ関数の収束を調べることにも繋がるということになります。
上の式についても、$${x=log(2n)}$$と、$${x=log(2n+1)}$$の両方の場合において、全ての$${e^{-sx}}$$の値の総和が0になる時は、当然ながらリーマンゼータ関数の値は0になることが推測できます。
また、①式は
$${x=log(2n)}$$と、$${x=log(2n+1)}$$
の場合のふたつの総和の差を示していますので、ふたつの場合の総和の差が0になる場合も、リーマンゼータ関数は0の値になります。また、ふたつの場合のうち、一方に-1を掛けて比較してみると、両者はほぼ実数軸に対して対称な分布を示すので、虚数部の値の総和が0になる時にも、対称的に相殺されて、リーマンゼータ関数が0になることが考えられます。

次の節では、対数螺旋の性質を利用して、sの実数値が0.5の時にリーマンゼータ関数が0になることを推論で示していきます。

自明でない零点の実数値に0.5しか見つからない理由を、推論計算で推論してみる(後半)

まず最初に、対数螺旋の性質として自己相似で、任意の倍率で拡大、縮小したものは、適当な回転によって元の関数と一致する。という性質があります。
特に、極座標表示で極座標表示で$${r(θ)=ae^{bθ}}$$となる対数螺旋であれば、$${e^{2πb}}$$倍に拡大したものは回転することなしに、元の螺旋と一致することがわかっています。
そのため、まずは、$${e^{-sx}}$$が$${x=log(2n)}$$と、$${x=log(2n+1)}$$のふたつの場合において、自己相似な座標の分布が成り立つかどうかを調べていきます。

複素数sをa+biとします。$${e^{-sx}}$$の対数螺旋としての表現は、$${e^{-ax}(cosbx-isinbx)}$$となります。対数螺旋を拡大するには、元の式にそのまま倍率の定数を掛ければいいので、自然数、n,mについてx=lognの関数を$${e^{\frac{2πa}{b}}}$$した場合にx=logmの場合で等号が成り立つmがあることを示していきます。

$${e^{\frac{2πa}{b}}×e^{-a logn}(cos(b logn)-isin(b logn))=e^{-a logm}(cos(b logm)-isin(b logm))}$$

$${e^{(-a-bi)logn+\frac{2πa}{b}}=e^{(-a-bi)logm}}$$

$${(-a-bi)logn+\dfrac{2πa}{b}=(-a-bi)logm}$$

$${log\dfrac{n}{m}=-\dfrac{2πa}{b(a+bi)}}$$

$${log\dfrac{m}{n}=-\dfrac{b(a+bi)}{2πa}}$$

$${\dfrac{m}{n}=e^{-\dfrac{b(a+bi)}{2πa}}}$$

n/mが有理数になれば、共に整数となるn、mが存在することになりますので、aとbの値により、ぴったり自己の対数螺旋と一致する倍率が存在することになります。必要条件は式変形の最後の式の右辺において、eの指数の虚数項がπの倍数であること。その上で同じ式が有理数となることです。

この条件を満たした時に、調べるnの範囲を大幅に小さく制限することが可能になります。
つまり、$${e^{-sx}、x=logn}$$   (②)の自然数nからmまでの範囲での、複素数平面全体もしくは実数軸上の値の総和が0であれば、自己相似の性質により、nもしくはmより大きい値での総和はすべて0になることがわかるからです。nもしくはmより小さい自然数での総和を調べることで、リーマンゼータ関数の値を調べることが可能になります。

ここでの理想はnもしくはmのうち、小さい方の自然数が1の時に成り立つ式になれば、他の一切の値を調べずに、収束値が0であることを言うことができる訳です。

では、ここからはひとまず、n、mの最小、最大からは離れて、nからmの範囲で、総和が0を取るかかどうかを調べてみます。

自然数n、m(n>m)が有理数p(p>1)によってm=npと表すことができ、nの場合の②式にp倍すると、mの場合の②式に一致するものと仮定します。
この時、nからmまでの②式の総和は次で表すことができます。
$${\sum\limits_{k=n}^{m} \dfrac{1}{(k)^s}=\dfrac{1}{n^s}+\dfrac{1}{(n+1)^s}+…+\dfrac{1}{m^s}(③)}$$

この式ですが、級数の総和の公式をそのまま使って求めるのはかなり難しいので、積分を使ってアプローチをすることにします。

ある実数の値を取る関数f(x)について、この関数がx=nからx=n+1の範囲で単調増加もしくは単調減少であることがわかっている場合は、次のことが言えると思います。

f(x)が単調増加の時、
$${f(n)<\int\limits_n^{n+1}f(x) dx < f(n+1)}$$

f(x)が単調減少の時、
$${f(n+1)<\int\limits_n^{n+1} f(x) dx < f(n)}$$

また、範囲を拡張して、f(x)について、x=nからx=n+q(q>0)の範囲で単調増加もしくは単調減少であることがわかっている場合は、次のことが言えると思います。

f(x)が単調増加の時、
$${f(n)<\dfrac{\int_n^{n+q}f(x) dx }{q}< f(n+q)}$$

f(x)が単調減少の時、
$${f(n+q)<\dfrac{\int_n^{n+q}f(x)dx }{q}< f(n)}$$

この関係性を使って③式を変形していきます。
③式のシグマの内部の各項の原点からの距離については、対数螺旋であることから、常にnに対して単調増加か単調減少となります。
s=a+biとして、|a|<1の時には単調減少、|a|>1の時には単調増加、|a|=1の時には常に変わりません。
よって、|a|<1のとき、
$${\int_n^{n+1} |e^{-s logx}| dx < |e^{-s logn}| < \int_{n-1}^n |e^{-slogx}| dx}$$

また、原点から実数軸の正軸方向に伸ばした半直線と、原点から対数螺旋上の一点に向けた直線の成す角はxに対して常に単調増加であり、2π以上の急激な角度の変化がある場合も、重複をそのまま数え上げれば関係性は変わらない。よって、重複をそのまま数え上げれば、
$${e^{-s logn}の偏角 < \int_n^{n+1} e^{-s logx} dxの偏角 < e^{-s log(n+1)}の偏角}$$
であることがわかります。
よって、|a|<1たる$${e^{-s logn}}$$の各項はnに対して、原点からの距離が、
$${\int_n^{n+1} |e^{-s logx}| dx<|e^{-s logn}|<\int_{n-1}^n |e^{-slogx}| dx}$$(|a|>1の場合は逆)
偏角が、
$${\int_{n-1}^n e^{-s logx} dx の偏角 < e^{-s logn}の偏角 < \int_n^{n+1} e^{-s logx} dx の偏角 }$$の領域の中にあることがわかる。特にnが十分大きければ、この領域は十分狭く、$${e^{-slogn}}$$と$${\int_n^{n+1} e^{-slogx} dx }$$はほぼ同値であると考えて良い。ということで、十分に大きな、自然数nにおいて、$${e^{-s logn}}$$を$${\int_n^{n+1} e^{-s logx} dx }$$に置き換えて総和を求めます。

nからmまで、自然数は連続する整数であるので、nからmまでの領域の積分は③式の総和に置き換えても、ほぼ同値であるとする。
$${\int_n^m e^{-s logx} dx}$$

$${=\int_n^m x^{-s } dx }$$

$${=\left[\dfrac{x^{-s+1}}{-s+1}\right]_n^m}$$

$${=\dfrac{m^{-s+1}-n^{-s+1}}{-s+1}}$$

$${=\dfrac{(np)^{-s+1}-n^{-s+1}}{-s+1}}$$

$${=\dfrac{n^{-s+1}(p^{-s+1}-1)}{-s+1}}$$
ここまで式変形ができる。
この式が常に0を取るためには、sが固定値であるため、
$${p^{-s+1}-1}$$が0を常に取れば満たすことができる。
ここで、m/nについて求めた式をpに代入して、値が0になるa,bを求める。

$${e^{-\dfrac{b(a+bi)}{2πa}×(-a-bi+1)}-1=0}$$

$${-\dfrac{b(a+bi)(-a-bi+1)}{2πa}=0}$$

$${-\dfrac{b(a^2+b^2+2abi-bi-a)}{2πa}=0}$$

$${a^2+b^2-a=0,2abi-bi=0}$$

a=1/2,b=±1/2が求まりました。ただし、この解は、pが有理数であるという条件をみたせません。そのため、この数値は零点とはなり得ません。
しかし、bの値に依らず、a=1/2で$${p^{-s+1}-1}$$が常に実数値を取ることがわかりました。
そのため、後は、$${\dfrac{n^{-s+1}}{-s+1}}$$が常に実数値を取ることができれば、十分に大きなnにおいて、③式が常に実数値を示すaとbが求まることになります。この解は実数軸に対して対称なので、大雑把に、
$${\int\limits_n^m e^{-s logx} dx}$$
のnより大きい自然数での項の総和が0になると考えられます。(④)

$${\dfrac{n^{-s+1}}{-s+1}}$$

$${=\dfrac{e^{(-a-bi+1)logn}}{-a-bi+1}}$$

$${=\dfrac{e^{(-a+1)logn}(cos(b logn)-isin(b logn))(1-a+bi)}{(1-a)^2+b^2}}$$

$${(cos(b logn)-isin(b logn))(1-a+bi)}$$の部分が全て実数になるbを求める。a=1/2とする。
$${cos(b logn)bi-isin(b logn)(1-a)=0}$$

$${cos(b logn)b=sin(b logn)/2}$$

$${tan(b logn)=2b}$$

この式を満たすnとbの組み合わせがあれば、④の条件を満たすことになります。
この条件を満たすn、bの組み合わせにおいて、自己相似の性質により、n/pからnの自然数までにおいて、③式が0になることを導くことができます。n/pが自然数になる場合に限って、n/p以降の全ての対数螺旋上の点の座標の総和が0になることが導けます。これは、n/pよりも小さい数n/p/p以降についても、同じことが言えるため、1になるまで再帰的に遡ることができます。n/p/p…を1になるまで遡った時に、ちょうどぴったりn/p/p…=1になることがわかった場合について、全ての②式の総和が0になることが導けます。これが、リーマンゼータ関数の零点であると考えられます。

(追記)
最後の結論は、最終的には間違っています。しかし、修正が厳しく、これ以上の深い答えが見つけられなかったので、このままにしておきます。