4

Rubyからゼータ関数を計算して、シキノートさんのサイトに掲載されている下記の画像のようなグラフを作りたいと考えています。

シキノートさんより引用

まずはnumo-gslのzeta関数を参照したのですが、GSLのzeta関数は引数として実数しか取れないのではないかという印象を受けました。(これは私が間違っているかも知れません)

Rubyから複素数を引数に取るようにしてzeta関数を計算できる方法はありませんか?
(Rubyから呼び出せれば実装がRubyである必要はなく、できれば速度の出る方法がいいです)

0

1 件の回答 1

5

質問欄で紹介されている リーマンのゼータ関数の数値計算コード(複素平面) には 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 function

3
  • このような素晴らしい解答が頂けると予想していなかったので大変嬉しく思っています。zeta_binomial メソッドで使われている whilerubocop の自動添削で引っかかりますので loop do ~ end に差し替えました。
    – kojix2
    2020年3月11日 18:56
  • metropolisさま、連絡方法が他にないのでここに書きます。上記のriemann zetaを計算するコードをGR.rbのサンプル に収載してよろしいでしょうか。metropolisさまとshikinoさまのクレジットを入れましたが、もしも不都合がございましたらご連絡ください。
    – kojix2
    2020年3月13日 1:11
  • 1
    @kojix2 はい、ご自由にどうぞ。ただ、計算精度にはご注意下さい。
    – user39889
    2020年3月13日 1:46

この質問に回答するには、ログインする必要があります。

求めていた回答ではありませんか? のタグが付いた他の質問を参照する。