ゼータ関数を描いてみる

最近ゼータ関数芸人の辻君が熱い(http://tsujimotter.hatenablog.com/entry/riemann-zeta-function )。リーマン予想については知っていても、その零点が素数階段に関係していることは知らなかった。素数の音楽(マーカス・デュ・ソートイ著,冨永 星 翻訳 http://amzn.to/1QnCAvM )を読んでみたが、翻訳も素晴らしくお勧め。ゼータ関数には少しばかり思い入れがあるので、ipython notebookでちょっと試してみる。このページはダウンロードして試せます。@edy555

ipython notebookでグラフを描く準備

In [2]:
%matplotlib inline
%config InlineBackend.figure_format='svg'
import numpy as np
import pylab as pl

Pythonでの複素数の取り扱いを確認

In [3]:
1+1j
Out[3]:
(1+1j)
In [4]:
(1+1j)**2
Out[4]:
2j
In [5]:
2**2j
Out[5]:
(0.18345697474330172+0.9830277404112437j)
In [6]:
abs(1+1j)
Out[6]:
1.4142135623730951

まずはナイーブな計算。 $\displaystyle \zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^s}$

In [7]:
def zeta0(s):
    x = 0
    for i in range(1,1000000):
        x += 1 / (i ** s)
    return x
In [8]:
zeta0(0.5+14.134j)
Out[8]:
(35.44160211530015+61.18274604837292j)

当然収束しない

scipy.specialにゼータ関数があるけど、複素数では使えないようだ

In [26]:
#import scipy.special
#scipy.special.zeta(0.5+14.134j, 0)

解析接続した関数で計算する

辻君のページから引用(http://tsujimotter.hatenablog.com/entry/riemann-zeta-function )

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

Pythonでは計算に必要なコンビネーションはscipy.miscにある

In [11]:
import scipy.misc as scm
scm.comb(1000000, 10)
Out[11]:
2.7556079168595421e+53

数値が小さい間はexact=Trueを付けると速い

In [12]:
scm.comb(10,2,exact=True)
Out[12]:
45

辻君のrubyのコードを素直に変換。(実はrangeにハマって悩んだ。終端がrubyはinclusive, pythonはexclusive)

In [13]:
from scipy.misc import comb
LOWER_THRESHOLD=1.0e-6
UPPER_BOUND=1.0e+4
INFINITY=300

def zeta(s):
   outer_sum = 0.0
   prev = np.inf
   if s == 1:
      return np.inf
   for m in range(1,INFINITY):
      inner_sum = 0.0
      for j in range(1, m+1):
         k = comb(m-1, j-1, exact=True) * j**(-s)
         inner_sum += k if j%2 else -k
      inner_sum = inner_sum * 2**(-m)
      outer_sum += inner_sum
      if abs(prev - inner_sum) < LOWER_THRESHOLD:
         break
      if abs(outer_sum) > UPPER_BOUND:
         break
      prev = inner_sum
   return outer_sum / (1-(2**(1-s)))

最初の零点付近で試す

In [14]:
zeta(0.5+14.134j)
Out[14]:
(9.031587952272404e-05-0.000568140412905536j)

クリティカルライン上の区間で計算してみる。numpy用に関数を変換してから計算。ついでに計算時間も計測

In [15]:
x = np.linspace(0,50,1000)
%time y = np.vectorize(zeta)(0.5 + x*1j)

CPU times: user 6.38 s, sys: 83.9 ms, total: 6.47 s
Wall time: 6.46 s

グラフにしてみる

In [16]:
#x = np.linspace(0,30,1000)
#%time y = np.vectorize(zeta)(0.8 + x*1j)
pl.grid()
pl.plot(x, np.abs(y))
pl.figure()
pl.grid()
pl.plot(np.real(y), np.imag(y))
Out[16]:
[<matplotlib.lines.Line2D at 0x11151d050>]

面で計算してみる

In [22]:
re = np.linspace(-5, 5, 51)
im = np.linspace(-40, 40, 101)
IM,RE = np.meshgrid(im, re)
%time Z = np.vectorize(zeta)(RE + IM*1j)

CPU times: user 20.4 s, sys: 218 ms, total: 20.6 s
Wall time: 20.5 s
In [27]:
pl.imshow(np.abs(Z), vmin=0, vmax=5, aspect=0.4)
Out[27]:
<matplotlib.image.AxesImage at 0x11780afd0>

3Dでプロットしてみる。inlineを外すと別ウィンドウでぐりぐり回せる

In [25]:
#%matplotlib 
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
fig = pl.figure(figsize=(10,6))
ax = Axes3D(fig)
ZZ = np.abs(Z)
ZZ[ZZ > 10] = 10
ax.plot_wireframe(RE, IM, ZZ)
#ax.plot_wireframe(RE, IM, np.log(np.abs(Z)))
#ax.plot_surface(RE, IM, np.log(np.abs(Z)), rstride=1, cstride=1, cmap='hot')
Out[25]:
<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x117668790>

引き続きGPUで計算させてみる予定