Googleの円周率31兆4000億桁までの計算について、きちんと解説

Facebook にシェア
Pocket
LINEで送る
このエントリーを Google ブックマーク に追加

さて、問題です(デーデン)

円周率とは何でしょう?

割と答えられない人が多いと思う。わかっている人でも「直径を1とするときの円の周囲の長さ」という人も多い。が、厳密には「円の直径と円周の長さの比率」というのが正しい。

とはいえ、円周率が「円の直径と円周の長さの比率」として使われる事のほうがいまとなっては少なくなっている。円周率は数学や物理学、広くは工学や経済学など幅広い分野で出現してくる不思議な数。もはや、「円の直径と円周の長さの比率」以外の意味があるのではないかとすら思えてくる数だ。

歴史的には古く、紀元前より円周率という考え方ではないものの、小学校で習う円の半径と面積の関係であったり、中学校以降で習う球の半径と表面積と体積の関係などが知られていた。昔の人、すごい。

円周率の厳密な値

現代数学では、円周率は厳密な値を求めることに対してあまり興味を持っていない。しかし、過去には様々な人達が円周率の厳密な値を求めることに注力していた。例えば、紀元前の数学者アルキメデスは、次のような方法で円周率を求めようとした。

円を円に内接する正多角形と円に外接する正多角形ではさむことで、円周は内接する多角形の辺の長さの和よりは大きく、外接する円の辺の長さよりは小さいという事実を利用して求めていく。正多角形の辺の数が増えれば増えるほど円に近づくため、徐々に円周に近づいていくというわけだ。この手法によって、円周率は3.14に近い数値であることが明らかとなる。ちなみにこの手法はアルキメデスの時代より2000年のあとに、ドイツの数学者・ガウスによって定式化される「はさみうちの原理」の元となっている。その後も様々な人によって3.14に近かったり遠かったりの円周率の値が導き出されることになる。西洋ばかりでなく中国の学者たちも円周率を計算し、有名なところでは祖冲之が3.1415929まで求めたことが知られている。かなりあとになるが、日本人も村松茂清が3.1415926まで求めている。またそれまでの数学的に言うといわゆる幾何学(図形を扱う分野)ではなく、例えばフランスの数学者・ド・モアブルは解析的に確率論を用いた手法で円周率を求めたりした。

厳密な値を求めることから少しそれて、円周率が無理数であることは知られている。無理数とは、ご存知だろうか。詳しくはSoftware Design3月号あたりを参照いただきたいが、無理数とは「分数で表せない数」である(個人的には実数のうち有理数でない数、と言いたい)。円周率が無理数であることは古代の哲学者(といってしまっていいのか疑問の余地はあるが)アリストテレスによって予想されていた。それを1761年にドイツの数学者・ランベルトによって証明される。つまり円周率は2000年くらい、よくわからない数だったのだ。ちなみに「「10桁で終了」 円周率ついに割り切れる」というニュースを真に受けている人とあったことがある。その後、様々な方法で円周率が無理数であることの証明がされている。無理数である以上、小数点以下の最後の数というのはないのである。

1900年代に入ると、コンピュータの登場で躍進する。天才として名高いフォン・ノイマンによるノイマン型コンピュータの開発により、2037桁まで計算された。その後はアルゴリズムの進歩はありつつも、コンピュータの進化によって今回のような莫大な桁数までの円周率が計算されている。

円周率はどうやって求める?

それほど膨大な円周率をどのように求めるか。最初は上で図示したとおり、図形によって計算していたが、徐々に解析的な手法での計算方法が編み出され、アルゴリズム化された。円周率が出てくる代表的な関係式でいうと、ライプニッツの定理(後にライプニッツよりその300年前にインドの数学者・マーダヴァによって考案されたことがわかった)と呼ばれる

\displaystyle \frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} + \dots = \sum_{n=0}^{\infty} \frac{(-1)^{n}}{2n+1}

やゼータ関数(バーゼル問題とも呼ばれる)

\displaystyle \zeta(2) = \frac{1}{1^2} + \frac{1}{2^2} + \frac{1}{3^2} + \dots = \frac{pi^2}{6}

などが知られているが、こういった式は円周率を求めるのには不向きである。円周率を求めるには、収束速度が鍵となる。収束速度とは簡単に言えば、一気にある値に定まる速度である。次の図を見てみよう。


青いグラフと赤いグラフだと、右に行くに従って赤いグラフのほうがはやく0に近づいていることがわかる。円周率を求める際にも、できる限り早く円周率に近づく数式を使ったほうが効率的である。

いくつか有意義な円周率の求め方を紹介する。まずインドの数学者・ラマヌジャンが発見したラマヌジャン法である。ラマヌジャンは次のような式を考案した。

\displaystyle \frac{4}{\pi} = \sum_{n=0}^{\infty} \frac{(-1)^{n}(4n)!(1123+21460)n}{882^{(2n+1)}(4^{n n!})^4}

Oops, ページを閉じないでほしい!

このような数式を求められたのはラマヌジャンくらいのもので、神が教えてくれない限り思いつくわけがない数式をラマヌジャンは発見している。しかもラマヌジャン法は精度が高く、また収束速度も早いことが知られている。

モンテカルロサンプリング法という方法もある。これは数式というよりもアルゴリズムで、正方形の中に正円を書きランダムに点を打って、円の中に点が入る確率を求めるというものである。プログラミングをしてもそれほど難しくなく実装できる。

また、つい最近まで使われていたかなり精度も高く収束速度も早いアルゴリズムで、ガウス=ルジャンドルのアルゴリズムがある。詳しくはWikipediaを参照してもらえるとわかるが、ラマヌジャン法よりもとてもシンプルに求めることができる。実際に、2009年までの円周率の最高桁数、約2兆6000億桁はガウス=ルジャンドルのアルゴリズムで求められた。

では今回はどのような方法で求められたのだろう。今回はアメリカの数学者・チュドノフスキー兄弟によって考案されたチュドノフスキーのアルゴリズムが使用された。嫌にならないでほしい。チュドノフスキーのアルゴリズムは次のとおりである。

\displaystyle \frac{1}{\pi} = 12 \sum_{n=0}^{\infty} \frac{(-1)^{n} (6n)! (545140134n + 13591409)}{(3k)!(k!)^{3} (640320)^{\frac{3k+3}{2}}}

記事を打っているだけでも嫌になってしまったが、今回の結果はこの数式が使われた。Wikipediaに書いてあったため、一応Pythonでのコードの例も記載する。

from decimal import Decimal as Dec, getcontext as gc

def PI(maxK=70, prec=1008, disp=1007): # parameter defaults chosen to gain 1000+ digits within a few seconds
    gc().prec = prec
    K, M, L, X, S = 6, 1, 13591409, 1, 13591409
    for k in range(1, maxK+1):
        M = (K**3 - 16*K) * M // k**3 
        L += 545140134
        X *= -262537412640768000
        S += Dec(M * L) / X
        K += 12
    pi = 426880 * Dec(10005).sqrt() / S
    pi = Dec(str(pi)[:disp]) # drop few digits of precision for accuracy
    print("PI(maxK=%d iterations, gc().prec=%d, disp=%d digits) =\n%s" % (maxK, prec, disp, pi))
    return pi

Pi = PI()
print("\nFor greater precision and more digits (takes a few extra seconds) - Try")
print("Pi = PI(317,4501,4500)") 
print("Pi = PI(353,5022,5020)")

 プログラムに起こすと少しは簡単に見えるかもしれない。単純にループして計算しているだけだ。単純に……。

円周率を求めるには、アルゴリズムの選定はかなり重要だ。しかし現実的にはむしろ、コンピュータの性能がものをいうのである。Googleは演算に特化したスーパーコンピュータを持っている。これが今回の結果につながったことは街がないと言える。ただし、すごいコンピュータを使えばいいというわけでもない。コンピュータ内でのCPUやメモリの効率化や分散方法など、様々な条件を整えなければならない。数学だけでなく、技術力も試される世界だ。

計算結果って正しいの?

今回、円周率は約31兆桁計算された。つまり、これまで人類が知ることができなかった桁の円周率がわかったということである。しかし、これは本当に正しいのか?正しいと誰が判断したのか?円周率の正しさを証明するために使われるのはまた数学だ。

今回の検証では、2つの数式によって確認された。1つはベラードの公式、もう一つはベイリー=ボルウェイン=プルーフの公式(BBP法) だ。ベラードの公式は、BBP法よりも高速に計算できることが知られている。

ベラードの数式は次のとおりである。もうここまで読んだのだから付き-合ってほしい。

\displaystyle \pi = \frac{1}{2^{6}}\sum_{n=0}^{\infty} (\frac{(-1)^{n}}{2^{10n}}-\frac{2^{5}}{4n+1}-\frac{1}{4n+3}+\frac{2^{8}}{10n+1}-\frac{2^{6}}{10n+3}-\frac{2^{2}}{10n+5}-\frac{2^{2}}{10n+7}+\frac{1}{10n+9})

また、BBP法は次のとおりである。

\displaystyle \pi = \sum_{n=0}^{\infty}[\frac{1}{16^{n}}(    \frac{4}{8k-1}    - \frac{2}{8k+4}    - \frac{1}{8k+5}    - \frac{1}{8k+6})]

よく打ち込んだと褒めてほしい(褒めなくて良い)。これらの数式は高速で円周率を計算できる。「ではこれを使えばいいのではないか?」と思うだろうが、そうは簡単には行かない。ベラードの公式とBBP法は両方とも、任意の桁数目の円周率を高速に求めることができる数式なのである。そのため、すべての桁数を求めることはできるが、その1桁ずつ求めていく必要があるため時間がかかってしまう。そのため、他のアルゴリズムで計算された数式の確認にのみ用いられるのである。

最後に

数学的にはそれほど新しい技術が使われていないが、Googleの世界最高峰の人材と技術が成し遂げた結果である。漠然と「円周率を求めただけじゃん」と思ってしまうが、上でも言ったとおり、数学的には具体的な円周率を知ること自体にそれほど興味はない。円周率を求めることは、コンピュータの性能を計る手段であり、それが今までの結果を圧倒する計算スピードとなったのである。

もし今後量子コンピュータが完成すれば、おそらく円周率の桁数は兆とか京とかそのレベルではない計算ができるようになるだろう。「2位ではだめなんですか?」とか「三角関数なんて誰が使うんですか?」などと政治家が言っている国では、とても追いつけるものではない。

Related posts

コメントを残す