リーマンのゼータ関数で遊び倒そう (Ruby編)

今日のテーマは「リーマンのゼータ関数」です。

リーマンのゼータ関数（以下，ゼータ関数）は，複素関数と呼ばれるタイプの関数です。複素数を変数にとって，複素数を関数値として返すので複素関数というのです。ゼータ関数は以下の式で定義されます。

$\displaystyle \zeta(s) = 1 + \frac{1}{2^s} + \frac{1}{3^s} + \frac{1}{4^s} + \cdots$

ゼータ関数は，実に魅力的な関数です。それは「オイラー積」や，それを応用した「リーマンの素数公式」を通して，数学のもっとも基本的な要素である「素数」と密接に結びついているためです。あとで紹介する「リーマン予想」という，数学史上最も難しいとされる未解決問題とも関連していて，ゼータ関数は tsujimotter にとってお気に入りの関数です。

今日は，このゼータ関数の形を「自分の力で描く」ための方法をご紹介します。いくら，ゼータ関数が魅力的といっても，自分の手で理解できないのであれば面白くありませんからね。また，単にゼータ関数を描く方法を紹介するだけでなく，先ほど挙げたリーマン予想についての理解も深まるような解説に挑戦したいと思います。ぜひ一緒にゼータ関数で遊んでみましょう！

数値計算用のプログラミング言語として，今回も Ruby を使っていますが，必要があればお好みの言語に置き換えて作ってみてください。

変数 $s$ について

本題からちょっとだけ脇道にそれますが，ゼータ関数の変数について触れておきます。

ゼータ関数を初めてみた方は，変数が $s$ となっていることにぎょっとするかもしれません。これは，高校のときに習った関数の $f(x)$ における，変数 $x$ の代わりに $s$ という記号を使っていると思ってもらえればよいです。

変数として $x$ という記号を使わない主な理由は， $x$ では実数のイメージが強いためです。ゼータ関数では， $x$ の代わりに $s$ を使うことが慣習になっています。

一般に，複素関数の変数は $z$ を使うことが多いのですが，ゼータ関数に限っては $s$ が好まれます。その理由は，ゼータ関数を発見したベルンハルト・リーマンが，最初の論文で $s$ を使ったからだと言われています。

解析接続とクリティカル・ライン

さて，まずはゼータ関数を「より扱いやすい形」に変えることから始めます。実は，上で定義したゼータ関数は，残念ながらこのままではあまり面白くありません。

なぜかというと，右辺の級数は，変数 $s$ の実部が $1$ 以下の点で発散してしまうからです。変数 $s$ の実部のことを ${\rm Re}(s)$ と表記すると，上の級数は ${\rm Re}(s) \leq 1$ において定義できないのです。

「まぁまぁ，別に ${\rm Re}(s) > 1$ で定義できれば十分じゃん。」と思うかもしれませんが，そうともいえないもっともな理由があるのです。

${\rm Re}(s) \leq 1$ の関数値にこだわる理由は，先に述べたリーマン予想と密接に関係しています。リーマン予想とは以下のような予想です。

リーマン予想：
解析接続されたゼータ関数の（非自明な）零点は，すべて ${\rm Re}(s) = \frac{1}{2}$ の直線上に存在する

「零点」とは，ゼータ関数の関数値がゼロになるような点のことです。
（「非自明な」という条件はひとまず置いておきましょう。あとで登場します。）

リーマン予想は，この零点が ${\rm Re}(s) = \frac{1}{2}$ の直線上にしか存在しない，ということを主張しています。この「 ${\rm Re}(s) = \frac{1}{2}$ の直線」のことを，別名 クリティカル・ラインと呼びます。リーマン予想が正しいとすると，ゼータ関数の零点を探そうと思ったら，このクリティカル・ラインを探せば良いというわけです。

したがって，リーマン予想を確かめるためには，少なくとも ${\rm Re}(s) = \frac{1}{2}$ 上のゼータ関数を計算できなければいけませんね。一方，ゼータ関数の最初の定義では， ${\rm Re}(s) > 1$ でしか定義できません。困りました。

この問題を見事に解決するのが「解析接続」というテクニックです。解析接続とは，簡単にいってしまえば「複素関数の定義域を広げるための方法」のことです。

複素関数論には「一致の定理」と呼ばれる定理があって，とある条件*1 を満たした領域内で２つの関数が同じ値をとるとしたら，両者は同じ関数とみなしてよいというのです。この定理を利用して，共通の定義域で同じ値をとりつつ，より広い範囲で定義可能な別の関数に「接続」することで，元の関数の定義域を広げることを「解析接続」と呼びます。

ゼータ関数は，解析接続によって「 $s=1$ を除く複素数平面全域」で定義できることが知られています。これが「解析接続されたゼータ関数」です。

ここまでの話を図でまとめると，こんなイメージになります。

f:id:tsujimotter:20150210010924p:plain:w500

解析接続されたゼータ関数は，リーマン予想で重要なクリティカル・ライン上できちんと定義されています。逆に，解析接続しないと，クリティカル・ライン上で値が計算できないので，面白みが半減してしまうというわけですね。

解析接続されたゼータ関数

ここまで理解すれば，解析接続されたゼータ関数の式を知りたくなってくるでしょう。

それでは，全複素数平面に解析接続できるその式をお見せしましょう！
これです！

ドン！

$\displaystyle \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}$

もちろん，この式は ${\rm Re}(s) > 1$ においては，最初の式と完全に一致します。この式で解析接続できることの説明は，当然気になるかと思いますが，この際脇においておきましょう。

どうしても気になる人は，以下の文献を読んでください。

~~http://homepage3.nifty.com/y_sugi/sp/sp13.htm~~
二項係数の級数とゼータ関数 （2016/08/16 追記：リンク先のURL が変更になったようです）

MathWorld の以下の記事の式 (21) も，見た目は若干異なりますが，まったく同じ式です。

Riemann Zeta Function -- from Wolfram MathWorld

とにかく，これもゼータ関数の解析接続された姿なのです。この式は $s=1$ 以外の全平面で定義可能です！

次の項では，この式を使った数値計算によって，ゼータ関数の姿を可視化していきたいと思います。

ゼータ関数を Ruby で書く

まず，式の形を見てください。右辺に

$\displaystyle \binom {m-1}{j-1}$

という式があるでしょう。これは，前回登場した二項係数です。いやあ，前回用意しておいてよかったですね。実は二項係数の記事は，この式を計算するためのものだったのですよ。

というわけで，前回作った "binom.rb" を同じディレクトリにおいておきます。

前回の記事：
binom.rb - 二項係数を３パターンで実行するスクリプト · GitHub

これを require で呼んであげましょう。あとは次のスクリプトの通り。

require 'Complex'
require './binom'

# 複素変数 s に対するリーマン・ゼータ関数を計算する関数
LOWER_THRESHOLD=1.0e-6  # 今回は 10^(-6) を使った。環境によって適宜変更
UPPER_BOUND=1.0e+4      # 今回は 10^(+4) を使った。環境によって適宜変更 
INFINITY=300

def zeta(s)
   outer_sum = 0
   prev = 1000000000

   for m in 1..INFINITY do
      # "m" に関する 和 を計算する
      #   = 2^{-m} \sum_{j=1}^{m} (-1)^{j-1} \binom {m-1}{j-1} \frac{j^{-s}}{1-2^{1-s}}
      inner_sum = 0
      for j in 1..m do
         c1 = ((j-1)%2==0) ? 1 : (-1)  # (-1)**(j-1) と等価。計算を簡略化した。
         c2 = binom(m-1, j-1)
         c3 = j**(-s)
         inner_sum += c1*c2*c3
      end
      inner_sum = inner_sum * (2**(-m)) / (1-2**(1-s))

      # 和に加える
      outer_sum += inner_sum

      # 差が閾値以下であれば「収束した」と判断して計算を終了する
      if (prev - inner_sum).abs < LOWER_THRESHOLD then
         break
      end

      # 大きな値は上記の収束判定がきかないため、UPPER_BOUND を超えると「打ち止め」して計算を終了する
      if outer_sum.abs > UPPER_BOUND then
         break
      end

      prev = inner_sum
   end

   outer_sum
end

注意点を何点かあげておきましょう。

まず，１つ目の注意点は，無限級数の収束に関してです。

今回のゼータ関数では「外側の和」と「内側の和」があります。これを二重 for ループで内側から順に計算しているのですが，外側の和は無限級数になっています。
基本的には収束する級数なので，ある程度まで部分和を計算したら打ち止めにすればよいわけですが，どこで止めるかは難しい問題ですね。

上のスクリプトでは，"INFINITY = 300" として，300 までの和をとるようにしています。ただ，これだと時間がかかりすぎます。

そこで，ある程度部分和が収束して誤差が減らなくなってきたら break するようにします。実装は "LOWER_THRESHOLD" を使った if 文のところです。今回の計算では，現在の部分和と一個前に計算した部分和との差が $10^{-6}$ 以下になったら打ち切るようにしています。

しかし，それでもまだうまくいかない場合があります。たとえば，収束先の値は大きいが収束のスピードは遅い，という場合です。この状況は，複素数平面の左半分を計算しているときによく発生します。この場合，あまり値が大きくなってもしょうがないので，これまで計算してきた部分和の絶対値がある程度を超えたら（今回の場合 "UPPER_BOUND"，つまり $10^4$ を超えたら）同様に打ち切るようにしています。あまり大きな値を計算したところで，描画は出来ませんからね。

もう１つ大きな注意点は，複素数の計算です。ゼータ関数の変数 $s$ は複素数であるため，複素数を計算する方法が必要です。一番問題になるのは，以下の箇所でしょう。

j**(-s)

これは， $j$ の $(-s)$ 乗を表す計算です。今回のスクリプトではの "Complex" というライブラリを require しているので，この式は「複素数のべき乗」として計算されます。

「複素数のべき乗」って何じゃそりゃ，って思うかもしれませんが，具体的には，複素指数関数・複素対数関数を使って次のように定義されます。

$j^{-s} = e^{-s\cdot\log{j}}$

Ruby ではこのあたりの面倒なことを，ライブラリが勝手にやってくれて非常に便利です。逆に，ほかの言語で実装する際には，この部分に気をつけるといいでしょう。

クリティカル・ライン上のゼータの値

さぁ，お待たせしました。それではゼータ関数を実際に計算していきましょう。

まず気になるのは，リーマン予想に関連した ${\rm Re}(s) = \frac{1}{2}$ のクリティカル・ライン上の点ですね。すなわち，次の式を計算していきたいと思います。

$\zeta\left(\frac{1}{2}+it\right)$

$t$ は任意の実数をとれますが，無限に計算するわけにはいかないので，今回は $t=0$ から $100$ まで $0.05$ おきに計算してみたいと思います。

少し気をつけたいのは，ゼータ関数は複素関数であるため，本質的に４次元である，ということです。どういうことかというと，変数 $s$ も複素数で，それに対する $\zeta(s)$ 関数値も複素数です。変数を実部・虚部の２次元で表して，関数値も２次元で表してしまうと，計４次元が必要になるというわけです。４次元はなかなか描画しづらいですよね。

そこで少し妥協して，ゼータ関数の関数値として，複素数の代わりにその絶対値をとることにしましょう。つまり， $\left|\zeta\left(\frac{1}{2} + it\right)\right|$ を計算するということです。

計算するためのスクリプトは以下の通りです。上のゼータ関数の定義の下に続けて書いてください。

File.open("critical-line.csv", "w") do |io|
   # クリティカル・ライン上を，実軸から離れる方向に， 0.05 おきに計算する．
   # t = 0 から t = 100 まで 
   0.step(100.0, 0.05) do |t|  

      s = 0.5 + Complex::I*t   # 1/2+it
      z = zeta(s)

      column = "#{s.real},#{s.imag},#{z.real},#{z.imag},#{z.abs},#{z.arg}"

      io.puts column
      puts column
   end
end

このスクリプトを実行すると，"critical-line.csv" という CSV ファイルが出力されます。2 列目に $t$ の値が，5 列目に対応する $\left|\zeta(\frac{1}{2}+it)\right|$ の値が現れます。

この CSV ファイルをグラフ描画ソフト（今回は gnuplot を使っています）でグラフ化するとこんな感じになります。

f:id:tsujimotter:20150209234419p:plain:w400

グラフの点が横軸に接する点が，

$\zeta\left(\frac{1}{2}+it\right) = 0$

の点，すなわち零点です！

最初の零点は，たしかに $t = 14$ 付近に存在していることが分かります。そのあとは，つかず離れずで登場します。ちゃんと零点って存在するんだな，ってことがわかりますね！

先ほどはゼータ関数の絶対値を計算しましたが，せっかくなので複素数のまま見てみたいですよね。そんな欲張りなあなたには，このような描画の方法もありますよ。

横軸に $\zeta\left(\frac{1}{2} + it\right)$ の実部，縦軸に $\zeta\left(\frac{1}{2} + it\right)$ の虚部をとります。先ほどの CSV でいうと，3列目と4列目に相当します。これを $t = 0$ から $t = 34$ まで計算して，線でつないだのがこのグラフです。

f:id:tsujimotter:20150210001822p:plain:w360

だいたい周期的にぐるぐると回っているように見えますが，図の中で５回ほど

$\zeta\left(\frac{1}{2} + it\right) = 0$

となる点を通過していることが分かります。これはもちろん零点のことです。

たしかに，先の図と見比べても， $t = 0$ から $34$ までの間に，５つ零点が存在しますから，その辺もちゃんと一致していますね。

全複素数平面上のゼータの値

ここまでは，クリティカル・ライン上に限定して関数の形状を部分的に見てきましたが，そろそろ全体像が見たいですね。

今度は虚軸方向だけでなく，実軸方向にもデータを取ってみましょう。
以下は，複素数平面の部分的な領域

$D = \{s\mid -40\leq {\rm Re}(s) \leq 40$ , $-40\leq {\rm Im}(s) \leq 40\}$

の上で，ゼータ関数を計算するスクリプトです。

File.open("commplex-plane.csv", "w") do |io|
   # 虚軸方向に 0.1 おきに [0.0, 40.0] の範囲で計算 
   (-40.0).step(40.0, 0.1) do |y|

      # 実軸方向に 0.1 おきに [-20.0, 20.0] の範囲で計算 
      # (ただし s=1 で「ゼロ割り」が発生するので、実軸方向の初期値を若干ずらしている。)
      (-40.00001).step(40.0, 0.1) do |x|  

         s = x + Complex::I*y  # x+iy
         z = zeta(s)

         column = "#{x},#{y},#{z.real},#{z.imag},#{z.abs},#{z.arg}"

         io.puts column
         puts column
      end

      io.puts ""
      puts ""
   end
end

コメントにも書きましたが， $s=1$ の点をドンピシャで通ってしまうと「ゼロ除算」が発生して止まってしまうので，その点は外すように工夫しています。

このスクリプトを実行すると "complex-plane.csv" というファイルが出力されます。
計算するときは，注意してください。収束判定があまりうまくないということもあって，計算にかなり時間がかかります。出力されるファイルサイズも 60MB を超える大きさです。

ではいよいよ，お待ちかねの全体像を見てみましょう。

今回は，実軸・虚軸・ゼータ関数の絶対値の３次元でグラフを書きたいと思います。といっても，そのまま三次元グラフを書くとよくわからないので，３軸目はカラーマップで表現します。データは，先ほど出力された CSV の 1 列目, 2 列名, 5 列目をとってきます。以下がその図です。

f:id:tsujimotter:20150211001216p:plain:w500

いやあ，ほれぼれするような図ですねぇ。デスクトップの壁紙にしたいぐらい。

見方を簡単に説明すると，ゼータ関数の絶対値が $0$ に近いところは黒っぽくなり， $4$ より大きければ白色になっています。 $s = 1$ の点は，ゼータ関数が発散するので白色になっており，縦軸よりやや右寄りに黒い点が存在しますが，この点は「一直線上に」存在しているように見えます。

ここで，リーマン予想を思い出しましょう。

リーマン予想は，実部が $\frac{1}{2}$ の直線上，すなわちクリティカル・ライン上に「非自明な零点」が存在する，というものでした。なるほど，今さっき「一直線上に」と言った部分は，まさにクリティカル・ラインといえそうですね。

クリティカル・ライン上の点を除くと，「黒い領域」つまり「零点の候補」は，中央左よりに１箇所だけ存在します。これをいったん無視して考えると，それ以外には，零点の候補となる黒い領域が存在しないようです。

つまり，これによってリーマン予想を図的に確認できたわけです。もちろん，今回計算したのは一部だけですので，他の部分に零点がある可能性は否定できません。でも，これでイメージはついたでしょう。もっとほかの領域が気になれば，数値計算の範囲を広げればよいのです。

最後に，先ほど説明を省いた，中央左よりの黒い部分について触れておきましょう。

リーマン予想の定義に「非自明な零点」という言葉があったかと思います。「非自明な」と言うからには「自明な」零点が存在するはずですよね。この中央左よりの黒い部分に存在する零点が「自明な零点」です。

実は， $s$ が「負の偶数」，つまり $s=-2n$ のとき，ゼータ関数 $\zeta(-2n)$ は，

$\zeta(-2n) = 0$

となることが知られています。この事実は，比較的簡単に示せるそうなのですが，紙面が足りないので今日は止めておきしょう。とにかく，場所が明らかな零点なので「自明な零点」なのですね。

それ以外の零点が「非自明な零点」です。こちらは，クリティカル・ライン上にしか存在せず，また登場するタイミングを予期できるような理論は今のところ存在しません。面白いですね！

おわりに

いかがでしたか。ゼータ関数をいろんな方法で可視化してみましたが，イメージは沸いてきましたでしょうか。

今回の記事を通して，

解析接続されたゼータ関数の（非自明な）零点は，すべて ${\rm Re}(s) = \frac{1}{2}$ の直線上に存在する

という「リーマン予想」の意味を直感的に納得してもらえたら幸いです。

こういうのって，自分で計算して描いてみて，はじめて分かるものだと思うのです。よくあるリーマン予想やゼータ関数の解説書では，今日描いた図とほぼ同じような図がたくさん出てきます。でもそれって，所詮ゼータ関数の一部を切り取ったものに過ぎません。図を描いた本人は理解しているのかもしれないけど，読んでいる側は描いてある部分以外のところはよくわからない。別の視点から見ないと納得できない人もいる。その意味で，本当に理解するためには，自分の好きな切り口で描いてみるのがいちばんだと思うのです。そういうわけで，今回の記事では「自分で描ける」という点に着目しました。

もう１つ，「描くためのツールは誰でも使えるものである」という点は重要だと思っています。ゼータ関数を描こうと思ったら，多くの人は Mathematica という商用ソフトを使いたいと思うでしょう。実際，Mathematica にはプリセットでゼータ関数を計算するモジュールが入っています。一方で，こうした商用ソフトは一般に高価です。数学科の学生でもないのに，こうしたソフトに手を出す人は多くないと思います。購入できたとしても，周りに聞ける人がいなければ，始めるのはハードルが高いかもしれません。今回紹介した Ruby であれば，誰もが無料で使うことが出来て，なおかつ，インターネットでも簡単に調べること出来て独学ではじめることができます。また，複素数が計算しやすいとか，桁落ちという概念が存在しないとか，そういうメリットもあったりします。

結局のところ，tsujimotter が目指すのは「高等数学の大衆化」なのだと思います。つまり，誰でも高度な数学で遊べるようにしたい，という壮大な目標です。「ゼータ関数」や「リーマン予想」といった専門的な数学に触れることは，一部の数学者だけの特権ではないと tsujimotter は思っています。非専門家であっても，興味を持って取り組めば理解できないものではないし，実際やってみるととても面白いです。「分野外だから」というだけで見向きもしないのは，あまりにももったいないです。「誰でも」と言わないまでも，少なくとも私と同じぐらい数学に興味を持った素人が，気軽に高等数学に触れて遊ぶことができる，そんなお手伝いをできればと思っています。

ソースコード

というわけで，みなさんが遊べるようにソースコードを Github に上げておきました。MITライセンスなので，著作権表示さえすれば，無許可で自由に利用が可能です。ぜひ使って遊んでみてくださいね。

GitHub - junpeitsuji/ruby_zeta: s=1 を除く全平面で解析接続された「リーマンのゼータ関数」の関数値を Ruby で計算する

"zeta.rb" がメインのスクリプトです。ターミナルから以下のように実行してください。

ruby zeta.rb

今回登場した図は，すべて Gnuplot で描画しています。上の Github リポジトリに入っている "plot-" で始まるファイルが，Gnuplot で描画するためのスクリプトです。Ruby で CSV を出力したあと，Gnuplot から以下のようにスクリプトを呼び出すと，今回の記事で使ったものと同じ図を png として生成してくれます。

gnuplot plot-******.gp

それでは，今日はこの辺で。

参考文献

今日紹介した図のいくつかは，この本にも載っています。もちろん，書き方までは載っていませんが。笑

素数に憑かれた人たち ~リーマン予想への挑戦~

作者:John Derbyshire,松浦俊輔,ジョン・ダービーシャー
発売日: 2004/08/26
メディア: 単行本

計算した零点の値が，本当に正しいかどうか確認したい方のために，以下のサイトを紹介します。こちらは，Andrew Odlyzko という計算器科学の専門家のサイトで，ゼータ関数の零点を 100,000 個載せたリストが公開されています。言わば，零点のデータベースです。

Andrew Odlyzko: Tables of zeros of the Riemann zeta function
http://www.dtc.umn.edu/~odlyzko/zeta_tables/index.html

数論の専門家，ブライアン・コンリーが書いたリーマン予想についての解説記事です。今回の調べものをしてて見つけました。英語です。

J. Brian Conrey, "The Riemann Hypothesis"
http://www.ams.org/notices/200303/fea-conrey-web.pdf

ゼータ関数の応用「リーマンの素数公式」についてまとめたスライドです。

「明日話したくなる「素数」のお話」第1回プログラマのための数学勉強会 #maths4pg from Junpei Tsuji

二項係数を求める関数の作り方 (Ruby編) - tsujimotterのノートブック

リーマンの素数公式を可視化する - tsujimotterのノートブック

リーマンの素数公式の可視化アプリがパワーアップしました - tsujimotterのノートブック

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

ゼータ関数のオイラー積 - tsujimotterのノートブック

追記：

@edy555 さんが python でゼータ関数でかくコードを作ってくれました。"ipython notebook" というフレームワークを使って Web 上でグラフも公開してくれています。ipython notebook 便利そうですねぇ。
http://nbviewer.ipython.org/gist/edy555/9988efd898d11fba2540