63654


Rubyからゼータ関数を計算して、[シキノート](http://slpr.sakura.ne.jp/qp/code-riemann-
zeta/)さんのサイトに掲載されている下記の画像のようなグラフを作りたいと考えています。

[![シキノートさんより引用](https://i.stack.imgur.com/3xQzO.jpg)](https://i.stack.imgur.com/3xQzO.jpg)

まずは[numo-gslのzeta関数](https://ruby-
numo.github.io/gsl/yard/Numo/GSL/Sf.html#zeta-
class_method)を参照したのですが、GSLのzeta関数は引数として実数しか取れないのではないかという印象を受けました。（これは私が間違っているかも知れません）

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


CC BY-SA 4.0
2020-03-06T05:41:28.640
63624
2020-03-12T03:10:48.883
32772
question
ruby narray
Rubyでゼータ関数を計算する方法を探しています


質問欄で紹介されている [リーマンのゼータ関数の数値計算コード(複素平面)](http://slpr.sakura.ne.jp/qp/code-
riemann-zeta/) には Fortran による実装コードが掲載されていますので、それを ruby から FFI で利用すると良いかも知れません。

それとは別に、ruby code だけで実装してみました。

「リーマンのゼータ関数の数値計算コード」に記載されている、ゼータ関数の級数展開を利用しています。

1.  
  
![](https://latex.codecogs.com/svg.latex?%5CLarge&space;%5Czeta\(s\)=%5Csum_%7Bn=1%7D%5E%7B%5Cinfty%7D%20%5Cfrac%7B1%7D%7Bn%5Es%7D~~\(1%20<%20s\))

2.  
  
![](https://latex.codecogs.com/svg.latex?%5CLarge&space;%5Czeta\(s\)=%5Cfrac%7B1%7D%7B1-2%5E%7B1-s%7D%7D%5Csum_%7Bn=0%7D%5E%7B%5Cinfty%7D%20%5Cfrac%7B1%7D%7B2%5E%7Bn+1%7D%7D%5Csum_%7Bk=0%7D%5E%7Bn%7D%20%7Bn%20%5Cchoose%20k%7D%20%5Cfrac%7B\(-1\)%5Ek%7D%7B\(k+1\)%5Es%7D)

記事でも言及されていますが、それぞれ収束の速さに違いがあり、引数(複素数)の値によって計算式を切り替えています。また、収束判定は `|p(k) -
p(k-1)|^2 < ε` で行っています。

**zeta_function.rb**

[code]

    ##
    ## 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
    }
    
[/code]

**実行結果**

[code]

    $ 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
    
[/code]

  
  
[![Zeta
function](https://i.stack.imgur.com/y1RRz.png)](https://i.stack.imgur.com/y1RRz.png)


CC BY-SA 4.0
2020-03-07T17:27:28.487
63654
2020-03-12T03:10:48.883
2020-03-12T03:10:48.883
32772
63624
answer