リーマンのゼータ関数の数値計算コード(複素平面)

ここではゼータ関数を数値計算で得るためのアルゴリズムとfortran90によるプログラムを載せています。

目次

  1. ゼータ関数の級数表示
  2. プログラム(fortran90)
  3. 計算データ
  4. グラフ
  5. 非自明、自明なゼロ点
  6. ゼータ関数の導関数とプログラム

ゼータ関数の計算方法


ゼータ関数は通常
・\(\rm{Re}(s)\gt 1\)に対しては直接計算を行い、\(\rm{Re}(s)\le 1\)の領域に対しては相反公式
\(
\displaystyle \zeta(s)=2^{s}\pi^{s-1}\sin\left(\frac{\pi s}{2}\right)\Gamma(1-s)\zeta(1-s)
\)

を使って求める事が推奨されます。
ガンマ関数を使うことを避けたい場合、全複素平面(\(s\ne 1\))で有効なゼータ関数の表示
\(
\displaystyle
\zeta(s)=\frac{1}{1-2^{1-s}}
\sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
\sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\)

を使うのが良いでしょう。

ゼータ関数の級数表示


ゼータ関数は
\(
\displaystyle
\zeta(s)=\sum_{n=1}^{\infty} \frac{1}{n^s} ,~~~(1\lt s)
\)

と定義されます。が、これはlogの収束であるため非常に収束が遅いです。
Euler-Knopp変換
\(\displaystyle
\sum_{n=0}^\infty a_n = \sum_{n=0}^{\infty} \frac{1}{(1+q)^{n+1}}
\sum_{j=0}^{n}{n \choose j} q^{n-j}a_j
\)

(q=1のときEuler変換)によって早く収束する級数に変換されます[4]。
より詳しい証明は[5]。

Riemann Zeta Function — from Wolfram MathWorld, Eq.(21)
によると、全領域(s=1以外)で収束するリーマンのゼータ関数の級数表示は、
\(
\displaystyle
\zeta(s)=\frac{1}{1-2^{1-s}}
\sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
\sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\)

と書けます。もう一つ形式があるのですが、そちらは収束が遅そうな形式だったので用いていません。

数値計算的には2の乗数とかはあまり計算したくありません。そこで、漸化式を考えます。
\(
\begin{align}
\zeta(s) &=\sum_{n=0}^{\infty} b_n \\
b_n &= \sum_{k=0}^{n}a_k^n \\
a_k^n &= \frac{1}{1-2^{1-s}}\frac{1}{2^{n+1}}{n \choose k}\frac{(-1)^k}{(k+1)^s}
\end{align}
\)
と\(b_n, a_k^n\)を定義します。\(a_k^n\;\; (n\ge k)\)を計算するために漸化式を考えると

\(
\begin{align}
a_k^n&=\frac{1}{2}\left(\frac{n}{n-k}\right)a_{k}^{n-1} \\
a_k^n&=\frac{k-n-1}{k}\left(\frac{k}{k+1}\right)^s a_{k-1}^n \\
a_0^0&=\frac{1}{2}\frac{1}{1-2^{1-s}}
\end{align}
\)

と導く事が出来ます。
計算手順としては、

zeta_alg_c

としていけばよいでしょう。
計算量を考えて赤の手順を少なくするようにしています。

fortran90によるプログラム


本当のゼータ関数の値との相対誤差\(\delta\)は約
\(\delta\lt 3\times 10^{-14} (Re(s)\lt 1)\),
\(\delta\lt 3\times 10^{-15} (1 \lt Re(s)) \)
負の領域の方が精度が悪いのは、相反公式から来ています。
相反公式のみを4倍精度演算で行えば
\(\delta\lt 7\times 10^{-15} (Re(s) \lt 1) \)
にまで精度が上がります。目立って精度はあがりません。

2017/01/12に大幅な更新をしました。
以前のものと比べ変更点は、
・4倍精度を使用していない
・計算速度は以前の54
になっているので、以前のものを使っている方は更新することを推奨します。
ガンマ関数に関しては複素ガンマ関数の数値計算をご参照ください。

function zeta(s)
  implicit none
  complex(kind(0d0)),intent(in)::s
  complex(kind(0d0))::zeta
  !
  !Checked : ifort   (success : version 16.0.2 20160204)
  !        : gfortran(success : version 4.8.4)
  !        : gfortran(success : version 4.4.7)
  !
  !http://slpr.sakura.ne.jp/qp/
  ! Author : sikino
  ! date: 2016/06/02, (yyyy/mm/dd)
  !     : 2016/06/03
  !     : 2016/06/04
  !     : 2016/06/05
  !     : 2016/12/25
  !     : 2017/01/12
  !     : 2017/01/14
  !
  integer::n,k
  integer,parameter::nmax=1000
  double precision::nq,nq1,ln1,ln2
  double precision,parameter::eps=1d-15
  double precision,parameter::pi=3.14159265358979324d0
  complex(kind(0d0))::ds,is,d,s0,tc,tt,qzeta
  complex(kind(0d0))::a(0:Nmax),b(0:nmax),cgamma
  external::cgamma
  !complex*32::zgamma
  !external::zgamma

  ds=dble(s); is=imag(s)
  if(ds.eq.1d0.and.is.eq.0d0)then
     zeta=1d+300
  elseif(ds.eq.0d0.and.is.eq.0d0)then
     zeta=-0.5d0
  else
     if(real(s).lt.1d0)then
        s0=1d0-s
     else
        s0=s
     endif
     
     tt=1d0/(1d0-2d0**(1d0-s0))    
     qzeta=0.5d0*tt

     a=cmplx(0d0,0d0,kind=8); b=cmplx(0d0,0d0,kind=8)
     b(0)=0.25d0*tt
     do n=1,nmax
        nq=n
        nq1=nq+1d0

        a(0)=b(0)
        do k=1,n-1
           a(k)=b(k)*nq/(n-k)
        enddo
        !----
        ln2=log(dble(n)/dble(n+1))        
        if(n.eq.1)then
           d=s0*ln2
        else
           d=d*ln2/ln1
        endif
        a(n)=-a(n-1)*exp(d)/nq
        ln1=ln2
        !----

        tc=cmplx(0d0,0d0,kind=8)
        do k=0,n
           tc=tc+a(k)
        enddo

        qzeta=qzeta+tc

        if(abs(qzeta).ge.1d0)then
           if(abs(tc/qzeta).le.eps)exit
        else
           if(abs(tc).le.eps)exit
        endif

        b=0.5d0*a
     enddo
     if(n.eq.nmax+1)write(6,*)" ***worning, didn't converge at zeta"

     if(real(s).lt.1d0)then
        zeta=(2d0**s)*pi**(-s0)*sin(0.5d0*pi*s)*qzeta*cgamma(s0)
        !zeta=cmplx((2q0**s)*acos(-1q0)**(-s0)*sin(0.5q0*acos(-1q0)*s)*qzeta*zgamma(s0*1q0),kind=8)
    else
        zeta=qzeta
     endif
  end if

  return
end function zeta

function cgamma(z)
  ! sikinote (http://slpr.sakura.ne.jp/qp)
  ! author : sikino
  ! date   : 2017/01/09
  !          2017/01/10
 
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::cgamma
 
  integer::i,n,M
  double precision,parameter::pi=3.14159265358979324d0
  double precision,parameter::e1=0.3678794411714423216d0
  double precision,parameter::ln2pi2=0.91893853320467274d0
  double precision,parameter::sq2pi=2.5066282746310005d0
  double precision::a(1:30),AB(1:14),dz,iz
  complex(kind(0d0))::q,s,w,r,q1,q2

  dz=dble(z)
  iz=imag(z)

  if(dz.eq.nint(dz).and.iz.eq.0d0)then
     if(dz.gt.0d0)then
        s=dcmplx(1d0,0d0)
        n=nint(dz)
        do i=2,n
           s=s*i
        enddo
     else
        s=1d+300
     endif
     cgamma=s
  else
     q=z
     if(dz.lt.0d0)q=1d0-z

     if(abs(iz).lt.1.45d0)then
        ! For x=arb, |y|\lt 1.45
        n=int(dble(q))
        s=1d0
        do i=1,n
           s=s*(q-i)
        enddo
        w=q-dble(n)

        a(1) = 1d0
        a(2) = 0.57721566490153286d0
        a(3) =-0.65587807152025388d0
        a(4) =-0.42002635034095236d-1
        a(5) = 0.16653861138229149d0
        a(6) =-0.42197734555544337d-1
        a(7) =-0.96219715278769736d-2
        a(8) = 0.72189432466630995d-2
        a(9) =-0.11651675918590651d-2
        a(10)=-0.21524167411495097d-3
        a(11)= 0.12805028238811619d-3
        a(12)=-0.20134854780788239d-4
        a(13)=-0.12504934821426707d-5
        a(14)= 0.11330272319816959d-5
        a(15)=-0.20563384169776071d-6
        a(16)= 0.61160951044814158d-8
        a(17)= 0.50020076444692229d-8
        a(18)=-0.11812745704870201d-8
        a(19)= 0.10434267116911005d-9
        a(20)= 0.77822634399050713d-11
        a(21)=-0.36968056186422057d-11
        a(22)= 0.51003702874544760d-12
        a(23)=-0.20583260535665068d-13
        a(24)=-0.53481225394230180d-14
        a(25)= 0.12267786282382608d-14
        a(26)=-0.11812593016974588d-15
        a(27)= 0.11866922547516003d-17
        a(28)= 0.14123806553180318d-17
        a(29)=-0.22987456844353702d-18
        a(30)= 0.17144063219273374d-19

        M=int(11.3*abs(w)+13)
        if(M.gt.30)M=30
        r=a(M)
        do i=M-1,1,-1
           r=r*w+a(i)
        enddo
       
        cgamma=s/(r*w)
     else
        ! For x=arb, |y|\gt 1.45
        s=1d0
        if(abs(q).lt.9d0)then
           do i=0,7
              s=s*(q+i)
           enddo
           s=1d0/s
           q=q+8d0
        endif

        AB(1) = 0.83333333333333333d-1
        AB(2) =-0.27777777777777778d-2
        AB(3) = 0.79365079365079365d-3
        AB(4) =-0.59523809523809524d-3
        AB(5) = 0.84175084175084175d-3
        AB(6) =-0.19175269175269175d-2
        AB(7) = 0.64102564102564103d-2
        AB(8) =-0.29550653594771242d-1

        q1=1d0/q
        q2=q1*q1
       
        r=AB(8)
        do i=7,1,-1
           r=r*q2+AB(i)
        enddo
        cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
     endif
     if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
  endif

  return
end function cgamma

ifort version 16.0.2 20160204, 12.1.0 20110811
gfortran version 4.8.4
で動くことを確かめています。

こちらに4倍精度用のゼータ関数のコードを置いておきます。
精度は何点かで試しましたが、30~32桁一致するくらいです。

ゼータ関数のデータ


計算結果のデータも置いておきます。
(11MB)riemanndata.tar.gz
(9MB)riemanndata.zip
それぞれ、中に
“riemannzeta_high_resolution.txt”
“riemannzeta_low_resolution.txt”
が入っています。
highの方は1列目から4列目にかけて
\({\rm Re}(s):[-20\sim 20],{\rm Im}(s):[-50\sim 50],{\rm Re}(\zeta(s)),{\rm Im}(\zeta(s))\)、0.1ずつ、小数点以下13桁で出力したもので、
lowの方は1列目から4列目にかけて
\({\rm Re}(s):[-20\sim 20],{\rm Im}(s):[-50\sim 50],{\rm Re}(\zeta(s)),{\rm Im}(\zeta(s))\)、0.2ずつ、小数点以下5桁で出力したものです。

x,y,z軸を\({\rm Re}(s),{\rm Im}(s),{\rm Re}(\zeta(s))\)にとりグラフにすると、
riemann_zeta_c
というグラフが得られます。

\({\rm Re}(s)\)が負で\({\rm Im}(s)\)が小さいところで発散が起こるようですので、注意をして使ってください。
何故倍精度でおかしい問題が起こっていたのか?今回の場合はオーバーフローでもなく、アンダーフローでもありません。
原因は桁落ちです。和と差によって激しく桁落ちしてしまうのです。

また、倍精度よりも計算が早いです。なぜならば目標精度までスムーズに収束してくれるからです。

相反公式を用いたので上記問題は発生しません。

[adsense1]

ゼータ関数の実部、虚部、絶対値のグラフ

ゼータ関数の実部、虚部、絶対値のグラフは以下の通りになります。

実部\({\rm Re}(\zeta(s))\)

zeta_real_quad
カラーマップ(横軸:\({\rm Re}(s)\),縦軸\(:{\rm Im}(s)\),カラー:\({\rm Re}(\zeta(s))\))
map_riemann_real

虚数\({\rm Im}(\zeta(s))\)

zeta_imag_quad_c
カラーマップ(横軸:\({\rm Re}(s)\),縦軸\(:{\rm Im}(s)\),カラー:\({\rm Im}(\zeta(s))\))
map_riemann_imag

絶対値\(|\zeta(s)|\)

zeta_abs_quad_c
カラーマップ(横軸:\({\rm Im}(s)\),縦軸\(:{\rm Re}(s)\),カラー:\(|\zeta(s)|\))

対数プロット(横軸:\({\rm Re}(s)\),縦軸\(:{\rm Im}(s)\),カラー:\((\log{10}(|\zeta(s)|\))
map_riemann_logabs

(おまけ)ゼータ関数のニュートン写像(詳しくはリーマンのゼータ関数の複素力学系へ)

非自明なゼロ点

上記の計算方法(4倍精度のプログラム)で、リーマンの非自明なゼロ点を考えます。
「非自明なゼロ点」はゼータ関数の実数軸上には無いゼロ点のことです。
リーマン予想とは
ゼータ関数\(\zeta(s)\)の「非自明なゼロ点」が\(s=\frac{1}{2}+ix\)(\(x\)は実数)上に全てある
というものです。\(\zeta(s)=0\)が満たされるときはゼータ関数の実部も虚部もゼロになるので、簡単に言えばゼータ関数の絶対値がゼロになる\(|\zeta(s)|=0\)と言うことです。
実際に\(|\zeta(s)|\)を可視化してみますとこうなります。
quadzeta_c

自明なゼロ点

僕は他のページでは自明な点に関する図を見たことがありません。非常に小さいからでしょう。
小さいものを見るには対数スケールで見るのが良いです。実際に書かせてみます。
quadzeta_log2_c
確かに小さいですね。対数にしないと明らかに見えません。

絶対値\(|\zeta(s)|\)を上の図のように表示させるgnuplotスクリプトはこちら。

set xr[-15:10]
set yr[-40:40]
set zr[0:2]
set view 69, 73, 2, 0.3
set view  equal xy
set ticslevel 0
   
set palette defined (-9 "purple",-6 "blue",-3 "turquoise", 0 "#f5f5f5", 3 "gold", 6 "orange", 9 "red")
set colorbox horiz user origin .1,.1 size .8,.04
set pm3d at b

splot "fort.11" u 1:2:(sqrt($3**2+$4**2)) every :3 w l lw 0.1

[adsense2]

クリティカルライン上での計算


クリティカルラインとは
\(
\displaystyle s=\frac{1}{2}+ix
\)

の線です。リーマン予想は、非自明なゼロ点が全てこの上にあるか?という問いかけであり、この計算と繋がっています。

倍精度型のアルゴリズム(zeta)では、\(x\sim 400\)程度まで計算することが出来ます。
これは、

相反公式の
\(
\displaystyle \sin(\frac{\pi}{2}s)
\)

の値がオーバーフローを起こすためです。

倍精度演算では約\(10^{306}\)まで扱えるので、
\(
s=\sin^{-1}(10^{306})\frac{2}{\pi}\sim 1+i449
\)
4倍精度では約\(10^{4932}\)まで扱えるので、
\(
s=\sin^{-1}(10^{4932})\frac{2}{\pi}\sim 1+i7230
\)
となります。
しかし、実際に計算してみると倍精度では確かに\(x=450\)程度まで計算可能ですが、4倍精度では\(x=1600\)程度までしか正しい値を返してくれません。どこかで何かが起こっているはずですが、解析はしてません。

実際に試してみますと、
倍精度では

となります。

ゼータ関数の導関数とプログラム


ゼータ関数の導関数も求めましょう。
導関数はゼータ関数の式をそのまま微分すればいいです。なぜなら、正則関数として解析接続されているからです。
実際に微分しますと、
\(
\displaystyle
\zeta'(s)=\frac{1}{(2^s-2)^2}\sum_{n=0}^{\infty} \frac{1}{2^{n+1}} \sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\left[ -2^s \{(2^s-2)\ln(k+1)+\ln4\} \right]
\)

変形して以下のように
\(
\begin{align}
\zeta'(s) &=\sum_{n=0}^{\infty} b_n \\
b_n &= \sum_{k=0}^{n}a_k^n \\
a_k^n &= (-1)^k {n \choose k} \left[ \frac{-2^s (k+1)^{-s}\{(2^s-2)\ln(k+1)+\ln4\}}{(2^s-2)^2} \right]
\end{align}
\)

となります。同じように漸化式を求めてやれば、
\(
\begin{align}
a_k^n&=\frac{1}{2}\left(\frac{n}{n-k}\right)a_{k}^{n-1} \\
a_k^n&=-\frac{k-n-1}{k}\left(\frac{k}{k+1}\right)^s\frac{(2^s-2)\ln(k+1)+\ln 4}{(2^s-2)\ln(k)+\ln 4} a_{k-1}^n \\
a_0^0&=\frac{1}{2}\frac{-2^s \ln 4}{(2^s-2)^2}
\end{align}
\)

として求められます。
実際の数値計算では、2番目の漸化式は余り有能ではないかもしれません。
ゼータ関数を求めるときと同じようにやるならば、必要なのは\(a_n^n\)を\(a_{n-1}^n\)から求める作業です。漸化式を用いても直接求めても計算量はさほど変わらない気がします。

導関数の相反公式


相反公式ももちろん存在します[3]。
解析接続されているので、元の相反公式の左辺、右辺をそのまま微分することで求められます。
\(
\begin{align}
-\zeta'(1-s)&=2(2\pi)^{-s}\left\{x\cos\left(\frac{\pi s}{2}\right)+y\sin\left(\frac{\pi s}{2}\right)\right\}\Gamma(s)\zeta(s) \\
&\hspace{3em} +2(2\pi)^{-s}\cos\left(\frac{\pi s}{2}\right)\left\{\Gamma(s)\zeta(s)\right\}’ \\
& x = -\ln(2\pi), \;\; y=-\frac{\pi}{2}
\end{align}
\)

\(
\Gamma'(s)=\Gamma(s)\psi(s)
\)
ここで\(\psi(z)\)はディガンマ関数

確かめはwolfram alphaを用いて以下のコマンドで行いました。
ReplaceAll[D[Zeta[x], x], {x -> -6+2i}]

fortran90によるプログラム


fortran90によるプログラムはこんな感じになるかと思います。
相反公式で、ゼータ関数、倍精度ガンマ関数とディガンマ関数を用います。
ゼータ関数は上の、倍精度ガンマ関数とディガンマ関数はつけておきました。詳しくはディガンマ関数の数値計算をご覧ください。
ここに記述しているのは倍精度のみです。

注釈

他のページでの例えば
リーマンのゼータ関数で遊び倒そう (Ruby編)
等では、漸化式を使わないで直接計算しています。このため、余り項の数\(n\)の数を増やす事が出来ないこと、倍精度と思われるため、桁落ちの問題が回避できない問題が起こります。
実際、私も直接計算するコードをfortran90で書きましたが、項の数を増やせず収束しない事、また計算時間が2~3倍余計にかかります。

参考文献
Riemann Zeta Function — from Wolfram MathWorld, Eq.(21)
draw a horizontal gradient somewhere at the bottom of the graph. -ColorBox

[3]Tom M. Apóstol “Formulas for Higher Derivatives of the Riemann Zeta Function”, Math. of. Comp vol 44, 169, 1985, pp. 223-232.

[4]Madhava-Gregory-Leibniz級数の収束の加速(収束の加速1)

[5]リーマンゼータ関数の級数表示による解析接続 インテジャーズ


「リーマンのゼータ関数の数値計算コード(複素平面)」への8件のフィードバック

  1. ζ関数を漸近展開の方法で計算してみたのですが、
    s=1/2+ix
    でxが大きいときは通常項の方をそれに見あった項数を足し込まないとうまく収束しないです。
    大体xの1/6以上必要の様だった。
    シキノ様の計算式に適用できるかどうかは判りませんけれど、もしかしたらxの大きさにあわせて計算項数を増やす必要があるのではないでしょうか。

    1. ゼータ関数の漸近展開とは何の話でしょうか?また、何がxの1/6で、何の計算項数を増やすべきと主張していますか?

  2. 岩波「数学公式Ⅲ」のζ関数の章に漸近展開公式が載っています。
    ζ(s)=Σn^(-s)+漸近展開項+剰余項
    でこの内のΣn^(-s)の部分の項数がxの1/6よりも小さいとどんなに漸近展開項を調整しても正規の値に近づかなかったです。(6というのは2πかもしれない)
    シキノ様のプログラムではnmaxを弄ってみてはということです。

    1. なるほど。理解しましたが、私の計算方法と全く違うアルゴリズムによる計算ですので、比較はできません。nmaxの意味も異なっています。

  3. こんにちは。
    あれっ?と思ったのですが
    臨界線上のアルゴリズムにおいて、反転しても
    1/2+it→1/2-it
    で虚部の符号が反転するだけだから、別に反転公式を使う必要はないのではないでしょうか?

  4. シキノ様のアルゴリズムで計算してみたのですが、
    nmax=2000 精度:無制限の条件です
    s 相対精度
    1/2+500i 378桁
    1/2+1000i 223桁
    12+1500i 118桁
    1/2+2000i 53桁
    1/2+2500i 16桁
    1/2+3000i 2桁
    という結果になりました。
    nmax=4000
    でも試してみましたが、ほぼ同じ相対比率でシフトする様ですね。
    計算量的にはオイラーマクローリンの和公式を使ったほうが少し有利なようです。

  5. ゼータ関数 のところはn=0ではなくてn=1ではないですか?これだと1/0になってしまうような・・・

    1. その通りですね。ご指摘ありがとうございます。
      訂正いたしました。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です