Rubyからゼータ関数を計算して、シキノートさんのサイトに掲載されている下記の画像のようなグラフを作りたいと考えています。
まずはnumo-gslのzeta関数を参照したのですが、GSLのzeta関数は引数として実数しか取れないのではないかという印象を受けました。(これは私が間違っているかも知れません)
Rubyから複素数を引数に取るようにしてzeta関数を計算できる方法はありませんか?
(Rubyから呼び出せれば実装がRubyである必要はなく、できれば速度の出る方法がいいです)
スタック・オーバーフローはプログラマーとプログラミングに熱心な人のためのQ&Aサイトです。すぐ登録できます。
登録してこのコミュニティに参加するRubyからゼータ関数を計算して、シキノートさんのサイトに掲載されている下記の画像のようなグラフを作りたいと考えています。
まずはnumo-gslのzeta関数を参照したのですが、GSLのzeta関数は引数として実数しか取れないのではないかという印象を受けました。(これは私が間違っているかも知れません)
Rubyから複素数を引数に取るようにしてzeta関数を計算できる方法はありませんか?
(Rubyから呼び出せれば実装がRubyである必要はなく、できれば速度の出る方法がいいです)
質問欄で紹介されている リーマンのゼータ関数の数値計算コード(複素平面) には Fortran による実装コードが掲載されていますので、それを ruby から FFI で利用すると良いかも知れません。
それとは別に、ruby code だけで実装してみました。
「リーマンのゼータ関数の数値計算コード」に記載されている、ゼータ関数の級数展開を利用しています。
1.
2.
記事でも言及されていますが、それぞれ収束の速さに違いがあり、引数(複素数)の値によって計算式を切り替えています。また、収束判定は |p(k) - p(k-1)|^2 < ε
で行っています。
zeta_function.rb
##
## Riemann zeta function - sum of a p-series
##
require 'cmath'
PI = Math::PI
Q = [
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
]
QLEN = Q.count
QLEN_G = QLEN - 2
SQRT_2PI = Math.sqrt(2.0*PI)
EPSILON = 1.0e-16
# Gamma function - Rosetta Code
# https://rosettacode.org/wiki/Gamma_function#C.23
def cgamma(z)
if z.real < 0.5 then
PI/(CMath.sin(PI*z)*cgamma(1.0-z))
else
z -= 1.0
x = (1...QLEN).inject(Q[0]){|s, i| s += Q[i]/(z+i)}
t = z + QLEN_G + 0.5
SQRT_2PI*(t**(z+0.5))*CMath.exp(-t)*x
end
end
def zeta_p_series(z)
def _zeta(z)
n, prev, sum = 2, Complex(1), Complex(1)
while n < 1e+6
add = 1.0 / (n**z)
sum += add
break if (add-prev).abs2 < EPSILON
n, prev = n+1, add
end
sum
end
if z.real > 1.0 then
_zeta(z)
else
(2.0**z)*(PI**(z-1.0))*
CMath.sin(PI*z/2.0)*cgamma(1.0-z)*_zeta(1.0-z)
end
end
##
## Riemann zeta function - binomial factor
##
def combinations(n, k)
return 1 if k == 0 or k == n
(k+1..n).reduce(:*) / (1..n-k).reduce(:*)
end
def zeta_binomial(z)
n, prev, sum = 1, Complex(0), Complex(0)
loop do
add = (1..n).inject(Complex(n)){|s, k|
s += combinations(n, k)*((-1)**k)/((k+1.0)**z)
} / (2**(n+1))
sum += add
break if (add-prev).abs2 < EPSILON
n, prev = n+1, add
end
sum / (1.0-(2**(1.0-z)))
end
##
## Riemann zeta function
##
def zeta(z)
case z
when Complex(1.0) then
Float::NAN
when Complex(0.0) then
Complex(0.5)
when Complex(-1.0) then
Complex(-1.0/12.0)
else
if z.real > -1.0 and z.imag.abs < 40.0 then
zeta_binomial(z)
else
zeta_p_series(z)
end
end
end
##
## ζ(z): Re(z)[-15.0, 5.0], Im(z)[-20.0, 20.0]
##
-15.step(5, 0.1) {|r|
-20.step(20, 0.1) {|i|
z = zeta(Complex(r, i))
continue if z == Float::NAN
puts '%.1f %.1f %e %e %e' % [r, i, z.real, z.imag, z.abs]
}
puts
}
実行結果
$ lscpu | grep -E '^Model name'
Model name: Intel(R) Core(TM) i5-8500T CPU @ 2.10GHz
$ lsb_release -ir
Distributor ID: Ubuntu
Release: 19.10
$ ruby --version
ruby 2.5.5p157 (2019-03-15 revision 67260) [x86_64-linux-gnu]
$ time { ruby zeta_function.rb > plot.dat; }
real 1m22.158s
user 1m22.147s
sys 0m0.008s
zeta_binomial
メソッドで使われている while
は rubocop
の自動添削で引っかかりますので loop do ~ end
に差し替えました。