Riemann zeta function の非自明零点
先日、「ビジュアル リーマン予想入門」という本を眺めていたのですが、色々やさしく分かりやすく書いてあるので、その結果の一部を再現してみることにします。
この本の最後の章で、von Mangoldt の明示公式を元にして、いくつかの素数から zeta 関数の 上の零点でピークが立つ関数が紹介されています。
（総和の最後の数Cが素数の時（階段の段に当たると）、値を半分にするのかもしれませんが、素数を避けることで考えないことにしますｗ）
実行結果
import matplotlib.pyplot as plt x, y = np.loadtxt('fort.10',delimiter=',', unpack=True) plt.figure(figsize=(16,8)) plt.grid() plt.ylim(-1, 2) plt.plot(x, y) plt.savefig('zero.png')
ピークの立っているところが、虚軸上の零点位置になります。
例 1/2+14i 近傍の最初の零点
14.13 と 14.14 の間に零点あり
14.11000 , 1.665981 14.12000 , 1.670522 14.13000 , 1.673154 14.14000 , 1.673850 14.15000 , 1.672614 14.16000 , 1.669438
Fortran ソース
はじめに素数を求め、それを元に von Mangoldt function を求めるのに必要な数を求め、それを用いて各値に対して $\Psi_C(x)$ を求める総和を取って計算しています。
program Riemann implicit none integer, parameter :: ntab = 10000 integer, allocatable :: kp(:) integer :: lambda(ntab) = 0, ip(ntab) ! ! sieve of Eratosthenes prime: block integer :: itab(ntab) = 0 integer :: i forall(i = 2:ntab) itab(i) = i do i = 2, int(sqrt(real(ntab))) itab(i*i::i) = 0 end do kp = pack(itab, itab /= 0) end block prime ! ! von Mangoldt function Lambda = Log(lambda) Mangoldt: block integer :: i, j, nj do i = 1, size(kp) nj = log(real(ntab)) / log(real(kp(i))) forall(j = 1:nj) lambda(kp(i)**j) = kp(i) end do end block Mangoldt ! ! zero point of zeta zeta0: block integer, parameter :: nmax = 10000 real , parameter :: xmax = 100.0 integer :: i, j real :: x(0:nmax), y(0:nmax), t, u x = [(xmax / nmax * i, i = 0, nmax)] forall(i = 1:ntab) ip(i) = i forall(i = 0:nmax) y(i) = sum(-log(real(lambda)) / ip * cos(log(real(ip)) * x(i)), lambda /= 0) do i = 0, nmax write(10, *) x(i), ',', y(i) end do end block zeta0 end program Riemann
