zeta 関数の非自明零点の計算

昨日求めた von Mangoldt 関数 $\Lambda(n)$ の定義から、zeta 関数の対数微分が

$\displaystyle{\zeta(s)'\over\zeta(s)}=-\sum_{n=1}^\infty {\Lambda(n)\over n^s}$

のように得られるようです。この対数微分を逆さまに見れば Newton 法での解の更新を連想させます。

$\displaystyle x_{n+1}=x_{n} - {f(x)\over f'(x)}$

$\Lambda(n)$ は 0 が多いのでこのまま行けないかと思ったのですが、やはり収束が悪くてどうにもなりませんでした。ただ検索すると Newton 法で zeta 関数の零点を求めている方がおられたので、それを参考に計算してみることにしました。

cosmos.art.coocan.jp cosmos.art.coocan.jp

zeta 関数計算のための級数

いくつか zeta 関数を計算しているネット記事などを見ますと、みな下記の同じ級数を用いているようなので、それを追うことにします。

$\displaystyle \zeta(s)={1\over1-2^{1-s}}\sum_{m=1}^{\infty}2^{-m}\sum_{j=1}^m(-1)^{j-1}\binom{m-1}{j-1}j^{-s}$

導関数は、項別微分すればよいのですが zeta そのものが出てくる項は零点で消えるはずなので零点を求めるだけならば zeta 関数に $(1 - 2^{1-s})$ を掛けたものを考えれば十分で計算は簡便化されます。上記サイトでもそのように計算しています。

以下の二つの級数を計算し、その比をとって Newton 法に持ち込みます。

$\displaystyle \sum_{m=1}^{\infty}2^{-m}\sum_{j=1}^m(-1)^{j-1}\binom{m-1}{j-1}j^{-s}$

$\displaystyle \sum_{m=1}^{\infty}2^{-m}\sum_{j=1}^m(-1)^{j-1}\binom{m-1}{j-1}j^{-s}\log(j)$

計算例

取りあえず理解のため、プログラムは何も考えずに、式の通り力ずくで進めます。

実行結果

初期値 (0, 14)

 (0.353681264120843912410326165114867,14.0559773885096512711848285131486)
 (0.488364131314399079717833284796826,14.1177261218294475436372086746189)
 (0.500128310656431658912558466227835,14.1344073607337982415586213507536)
 (0.500000067229512513257628043161170,14.1347252096452643219409359668890)
 (0.500000000000000246632265634545671,14.1347251417346863652486230667610)
 (0.500000000000000000000000000044723,14.1347251417346937904572519835665)
 (0.500000000000000000000000000000000,14.1347251417346937904572519835625)

正しく 0.5 + 14.1347251417346937904572519835624_702707842571156992i に収束しました。

初期値 (0, 101)

 (0.256069484490943666787489813775311,101.025358346443434023241782295442)
 (0.481942442980168055821830736687186,101.106742309285289934516211943197)
 (0.588559085383994993298609342005680,101.288337808967690842316625018029)
 (0.481494597190548538690508425123605,101.326114747152462162521421686812)
 (0.499363832700973876130295496302147,101.318313833220907873792593880834)
 (0.499999389309487223502748219132107,101.317852065106741980562206117842)
 (0.500000000000916858174656863132202,101.317851005734203265517874815704)
 (0.500000000000000000000015796736012,101.317851005731391228785440874372)
 (0.499999999999999999999999999999988,101.317851005731391228785447940292)

正しく 0.5 + 101.317851005731391228785447940292_30890633286638430089479992831871523i に収束しました。

www.dtc.umn.edu

Fortran プログラム

    program main
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, parameter :: kd = real128
        integer :: i
        complex(kd) :: s
     !
     !14.1347251417346937904572519835624_702707842571156992
     !
        ! zero point of zeta function by Newton's method
        s = (0.0_kd, 14.0_kd)
        do i = 1, 10
            s = s - z2(s) ! Newton's method
            print *, s
            if (abs(z2(s)) < 1e-32) exit 
        end do    
     
    contains
    
        complex(kd) function z2(s) result(res) !  {(1 - 2.0_kd**(1-s))zeta(s)} /  {(1 - 2.0_kd**(1-s))zeta(s)}'
            import, only : kd
            complex(kd), intent(in) :: s
            integer, parameter :: mmax = 250 
            integer :: m, j
            complex(kd) :: d0, e0, z0
            complex(kd) :: d1, e1, z1
            real(kd), allocatable :: bino(:)
            bino = [1.0_kd]
            d0 = (0.0_kd, 0.0_kd)
            d1 = (0.0_kd, 0.0_kd)
            do m = 1, mmax
                e0 = (0.0_kd, 0.0_kd)
                e1 = (0.0_kd, 0.0_kd)
                do j = 1, m
                    e0 = e0 - (-1)**j * bino(j) * j**(-s) 
                    e1 = e1 + (-1)**j * bino(j) * j**(-s) * log(real(j, kd))
                end do
                d0 = d0 + e0 / 2.0_kd**m
                d1 = d1 + e1 / 2.0_kd**m
                bino = [bino, 0.0_kd] + [0.0_kd, bino] ! binomial coefficients
            end do
            res = d0 / d1
        end function z2
    end program main

その他参考サイト

ja.wolframalpha.com

slpr.sakura.ne.jp

tsujimotter.hatenablog.com

おまけ　

よく見る (0.5, xi) 上の zeta 関数の絶対値のグラフ。

    program main
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, parameter :: kd = real128
        integer :: i
        complex(kd) :: s, ds, z, y1, y0, y2
        do i = 0, 1000
            write(10, *) i / 100.0, ',', abs(zeta(cmplx(0.5_kd, i / 1.0_kd, kd)))
        end do

    contains

        complex(kd) function zeta(s) result(res)
            import, only : kd
            complex(kd), intent(in) :: s
            integer, parameter :: mmax = 250 ! 100
            integer :: m, j
            complex(kd) :: c, d, e, w(mmax)
            real(kd), allocatable :: bino(:)
            bino = [1.0_kd]
            c = 1.0_kd / (1 - 2.0_kd**(1-s))
            d = (0.0_kd, 0.0_kd)
            do m = 1, mmax
                e = (0.0_kd, 0.0_kd)
                do j = 1, m
                    e = e - (-1)**j * bino(j) * j**(-s)
                end do
                d = d + e / 2.0_kd**m
                bino = [bino, 0.0_kd] + [0.0_kd, bino] ! binomial coefficients by Pascal's triangle
            end do
            res = c * d
        end function zeta

    end program main