最近ゼータ関数芸人の辻君が熱い(http://tsujimotter.hatenablog.com/entry/riemann-zeta-function )。リーマン予想については知っていても、その零点が素数階段に関係していることは知らなかった。素数の音楽(マーカス・デュ・ソートイ著,冨永 星 翻訳 http://amzn.to/1QnCAvM )を読んでみたが、翻訳も素晴らしくお勧め。ゼータ関数には少しばかり思い入れがあるので、ipython notebookでちょっと試してみる。このページはダウンロードして試せます。@edy555
ipython notebookでグラフを描く準備
%matplotlib inline
%config InlineBackend.figure_format='svg'
import numpy as np
import pylab as pl
Pythonでの複素数の取り扱いを確認
1+1j
(1+1j)
(1+1j)**2
2j
2**2j
(0.18345697474330172+0.9830277404112437j)
abs(1+1j)
1.4142135623730951
まずはナイーブな計算。 $\displaystyle \zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^s}$
def zeta0(s):
x = 0
for i in range(1,1000000):
x += 1 / (i ** s)
return x
zeta0(0.5+14.134j)
(35.44160211530015+61.18274604837292j)
当然収束しない
scipy.specialにゼータ関数があるけど、複素数では使えないようだ
#import scipy.special
#scipy.special.zeta(0.5+14.134j, 0)
Pythonでは計算に必要なコンビネーションはscipy.miscにある
import scipy.misc as scm
scm.comb(1000000, 10)
2.7556079168595421e+53
数値が小さい間はexact=Trueを付けると速い
scm.comb(10,2,exact=True)
45
辻君のrubyのコードを素直に変換。(実はrangeにハマって悩んだ。終端がrubyはinclusive, pythonはexclusive)
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)))
最初の零点付近で試す
zeta(0.5+14.134j)
(9.031587952272404e-05-0.000568140412905536j)
クリティカルライン上の区間で計算してみる。numpy用に関数を変換してから計算。ついでに計算時間も計測
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
グラフにしてみる
#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))
[<matplotlib.lines.Line2D at 0x11151d050>]
面で計算してみる
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
pl.imshow(np.abs(Z), vmin=0, vmax=5, aspect=0.4)
<matplotlib.image.AxesImage at 0x11780afd0>
3Dでプロットしてみる。inlineを外すと別ウィンドウでぐりぐり回せる
#%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')
<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x117668790>
引き続きGPUで計算させてみる予定