ジーゲルのZ関数を数値計算する - tsujimotterのノートブック

リーマン予想とかリーマンの素数公式とかの文献を調べていくと「ゼータ関数の零点を求めたいな」って気分になりますよね。

下記の Andrew Odlyzko のページに行けば、零点の生データを 100,000 個まで得ることができます。

Andrew Odlyzko: Tables of zeros of the Riemann zeta function

でもここでは、自力で計算したいと思います。

原理的に言えば、 $\Re(s)=\frac{1}{2}$ となる直線（以降、クリティカルストリップ）上のゼータ関数の値を順次計算していき、それがイコールゼロとなる点をすべて計算すれば、すべての零点が得られるはずです。

よって $t$ を実数として、 $y=\zeta\left(\frac{1}{2}+it\right)$ の値を求めていき、そのグラフを観察すればよいですね。

残念なことに、ゼータ関数をそのまま計算してしまうと関数値が複素数になってしまうため、描画や計算が面倒です。

ここでは、ジーゲルのZ関数を利用しましょう。Z関数はジーゲルという数学者が、手計算が書かれたリーマンの遺稿から「発見」した関数として知られています*1。

《リーマン・ジーゲルのZ関数》

$\displaystyle Z(t) = e^{i\vartheta(t)} \zeta\left(\frac{1}{2}+it\right)$

この関数は、

絶対値がゼータ関数のクリティカルストリップ上の値と一致する
実数関数である

という非常に都合の良い性質を持っています。

話はそれますが、今回の記事に向けた調べ物の過程で、Z関数の値を「音の振幅」に見立てて、それを出力して音声ファイルとして公開いる人を見つけましたよ。

Sound of Zeta : wrong, rogue and log

さながら「ゼータ関数の音楽」ということでしょうか。面白いことを考える人もいるもんですねえ。

Z関数の計算

Mathematica にはZ関数を計算するルーチンが搭載されているそうです。

次のブログのようにコマンドをたたけば、簡単にZ関数をグラフ描画できます。

RH祭り : wrong, rogue and log

でも、「一般の人」には Mathematica なんて高くて買えないですよね。「一般の人」である tsujimotter は、無料で使える Ruby 言語のスクリプトで Z関数を計算する方法を模索していこうと思います。

さて、いろんな文献を調べてみると、Z関数の厳密な計算は難しいことがわかります。

ただ、必ずしも厳密に計算しなくてもよければ、下記のような近似式が存在します。

《ジーゲルのZ関数の近似式》

$\displaystyle Z(t) \sim 2\sum_{n^2<\frac{t}{2\pi}} n^{-\frac{1}{2}}\cos(\vartheta(t) - t\log(n))$

tsujimotter 的には、ある程度の範囲でそれなりに近いグラフが書ければ満足なので、この近似を使って計算したいと思います。精度についてはあとで、Mathematica のグラフと見比べてみることにします。

ところで、この近似式の出どころですが、「明解ゼータ関数とリーマン予想」の第六章に載ってますので、気になる人はどうぞ。

明解ゼータ関数とリーマン予想 (KS理工学専門書)

作者: ハロルドエム・エドワーズ,鈴木治郎
出版社/メーカー: 講談社
発売日: 2012/06/21
メディア: 単行本（ソフトカバー）
購入: 2人クリック: 5回
この商品を含むブログ (1件) を見る

ほかにも次のような文献もインターネットで探せば見つかるので、この近似でよさそうだということは容易に確認できるかと思います。

参考：
http://www.math.yorku.ca/~akuznets/on_the_Riemann_Siegel_formula.pdf
http://www.phy.bris.ac.uk/people/berry_mv/the_papers/berry265.pdf
Riemann-Siegel Formula -- from Wolfram MathWorld

ここで、まだ定義が不明なのが、 $\vartheta(t)$ という関数ですね。
これまた「ジーゲルのシータ関数」という名前がついています。

これも厳密に計算するのは難しいようです。

そこで、Wolfram Mathworld の該当ページの下の方に記載されている漸近展開の式を参考にします。

Riemann-Siegel Functions -- from Wolfram MathWorld

ここに書いてある項では、残念ながら不十分であることが実験により確認できたので、次のページのリストを参考に項を付け足しました。

$\vartheta(t)$ の漸近展開における各項の分子のリスト http://oeis.org/A036282
$\vartheta(t)$ の漸近展開における各項の分母のリスト http://oeis.org/A114721

それが次の式です。

《ジーゲルの $\vartheta$ 関数の漸近展開》

$\displaystyle \begin{eqnarray} \vartheta(t) \sim \frac{t\log\left(\frac{t}{2\pi}\right)}{2} &-& \frac{t}{2} - \frac{\pi}{8} + \frac{1}{48t} + \frac{7}{5760t^3} \\ &+& \frac{31}{80640t^5} + \frac{127}{430080t^7} + \frac{511}{1216512t^9} + \cdots \end{eqnarray}$

以上の式を Ruby で計算するスクリプトはこちらです。

# ジーゲルのシータ関数
def siegelTheta t
	theta = 0.5*t*Math.log(t/(2*Math::PI)) - 0.5*t - Math::PI/8.0
	theta += 1/(48*t)
	theta += 7/(5760*t*t*t)
	theta += 31/(80640*t*t*t*t*t) 
	theta += 127/(430080*t*t*t*t*t*t*t)
	theta += 511/(1216512*t*t*t*t*t*t*t*t*t)

	return theta
end

# ジーゲルのZ関数
def siegelZ t
	sum = 0
	th = (t/(2*Math::PI))

	n = 1
	while n*n < th do
		sum += 2*Math.cos(siegelTheta(t) - t*Math.log(n))/Math.sqrt(n)
		n+=1
	end	
	
	return sum
end

# 0.01 刻みで [0, 100) の範囲を出力
intervals = 0.01
10000.times do |i|
	t = i*intervals + 0
	puts "#{t} #{siegelZ(t)}"
end

計算できたら、出力データを Gnuplot で描画してみましょう。

結果

得られたグラフは次の通りです。

f:id:tsujimotter:20140630175746p:plain

若干Mathematicaの結果と違うようです・・・。

とくに $t=0$ 付近の値が、まったく計算できていませんね。

tsujimotterの推測ですが、この誤差は $\vartheta(t)$ の近似精度が原因なのではないかと思います。

$\vartheta(t)$ は、 $\cos$ の中の「位相部分」を担当しているため、これが少しでもずれると $Z(t)$ の関数値が大きく変わってしまうのではないかと思われます。

零点の位置に関しては正確に調べていないですが、 $t=14$ 付近の点以降は、Andrewのものとおよそ一致しているように見えます。

あとがき

Z関数の近似式を改めて見返してみると、ゼータ関数が「三角関数」にみえてきますね。

ゼータ関数の専門家であられる黒川先生が、とある本で次のようにおっしゃっていたことを思い出しました。

「ゼータ関数を貶めようという気はないんですが、ゼータ関数は三角関数の一種だと思っているわけです。」

なるほど納得ですね。

とある本とは、次の本のことです。私もまだ買ってはいないのですが、ぱらっと見ただけでもなかなか面白そうな本でしたよ。

21世紀の新しい数学 ~絶対数学、リーマン予想、そしてこれからの数学~ (知の扉)

作者: 黒川信重,小島寛之
出版社/メーカー: 技術評論社
発売日: 2013/07/23
メディア: 単行本（ソフトカバー）
この商品を含むブログ (13件) を見る

*1:ここで、誤解のないようにジーゲルの仕事について補足しておかなければなりません。まず前提として、リーマンは完璧主義の性格から論文を書くための途中計算を残すことはほとんどありませんでした。彼の論文は簡潔すぎて、途中の過程がすっぽり抜けているのです。そのため、リーマンの真意を理解するためには彼の手計算のメモが必要不可欠でした。リーマンの死後、手計算のメモは彼の書斎に残っていたはずでした。しかしながら、リーマンの家政婦は、なんとそのメモを焼却処分してしまったのです。ジーゲルは、その焼け残った手計算メモのかけらをかき集め、つないで、足りない部分は補足してなんとか先の公式を発見したのです。ただでさえ、適当にかかれた手計算のメモの、しかも大半は焼失した中から、計算式を導き出す仕事は、ジーゲルの涙ぐましい努力なしには実現されなかったでしょう。元々はリーマンが発見したものだったわけですが、現在では「ジーゲル」の名前も残っているというわけです。