円周率世界記録樹立

　　　　2002/12/20 日立製作所 & 早稲田大学　後保範 ( Ushiro Yasunori )
　　　　2007/04/01　東京工芸大学　後保範 (Tokyo Polytechnic University )

1. 円周率世界記録達成の概要

　2002年11月24日13時39分に、10進表現で1兆2411億桁、16進表現で１兆307億桁の
円周率が東京大学のHITACHI SR8000/MPP で正確に求まったことを確認しました。
本結果は、東京大学情報基盤センターで12月6日10時におこなわれた記者会見で
同センターの金田教授から公表しました。また、12月8日のTBSの報道特集で
「新記録樹立！円周率計算の挑戦」と題して放映されました。この記録は、
これまでの金田教授による10進表現で約2061億の世界記録(1999年9月)を6倍更新
したものです。
世界記録達成までの過程と概要は円周率世界記録達成の概要を参照下さい。

2. 世界記録達成のための主な工夫点

　世界記録達成のための主な工夫点は円周率世界記録達成の工夫点を参照下さい。

3. 計算結果の抜粋

　10進の結果の抜粋は金田研究室の10進円周率結果の抜粋を参照して下さい。
　16進の結果の抜粋は金田研究室の16進円周率結果の抜粋を参照して下さい。

4. 結果の統計値

　10進の結果の分布は金田研究室の10進円周率結果の抜粋を参照して下さい。
　16進の結果の分布は金田研究室の16進円周率結果の抜粋を参照して下さい。

5. 世界記録達成に至る余談

　1兆桁挑戦への過程
　(1)　1998年秋：　DRM法考案。級数公式で1兆桁計算可能
　(2)　同年末：　DRM法論文とπ計算の打ち合わせ開始
　(3)　2000年7月：　DRM法π計算原理プログラム(ラマヌジャン方式)完成
　(4)　2001年6月：　1ノード版π計算プログラム(ラマヌジャン方式)完成
　(5)　同年 8月：　並列版完成(ラマヌジャン方式、128ノード版)
　(6)　同年10月：　2公式で16進5000億桁計算、3597億桁で不一致
　(7)　同年11月：　Arctan公式1ノード版完成(64ノード対応、I/O利用)
　(8)　同年12月：　3597億桁での不正原因判明、対策
　(9)　2002年7月：　64ノードArctan公式版の並列化着手
　(10)　同年10/4：　64ノードArctan版で16進1兆桁検証完
　(11)　同年11/24：　10進1.2兆桁への変換と検証完　

6. 円周率計算公式

　多数桁の円周率計算に使用される4種類の計算公式を示す。

7. 円周率計算プログラム

　FORTRANとCの5種類のπ計算プログラム を順に掲載。本プログラムは内容を理解
してもらうことを目的とし、計算時間及び使用メモリ量の効率化は目的としていない。
計算時間及び使用メモリ量の効率化を施した本格的なプログラムは金田研究室から
他のプログラムと合わせて公開する予定。

8. DRM法の概要

　有理数係数の級数関数の計算量をO(n²)からO(n・(log(n))^2～3)に削減する方式である。
ここで発生する多数桁の乗算はFFT(高速フーリエ変換)やFMT(高速剰余変換)を使用
してO(n・log(n))で計算する。
　DRM法の計算原理を下記に示す。
　　F(x) = [A₁/B₁+A₂/B₂] + [A₃/B₃+A₄/B₄] + [A₅/B₅+A₆/B₆] + [A₇/B₇+A₇/B₈] + …
　　　　　= [C₁/D₁+C₂/D₂] + [C₃/D₃+C₄/D₄] + …
　　　　　= [E₁/F₁+E₂/F₂]+ …
　　　　　= H₁/G₁… = …
　ここで、下記の式で項数を半分に桁数を2倍化するトーナメント式通分を続ける。
　　D_i=B_2i-1・B_2i, 　　C_i=A_2i-1・B_2i+A_2i・B_2i-1, 　　i=1,2,…,n
　　F_j=D_2j-1・D_2j, 　　E_j=C_2j-1・D_2j+C_2j・D_2j-1, 　　j=1,2,…,n/2
　　H_k=F_2k-1・F_2k, 　　G_k=E_2k-1・F_2k+E_2k・F_2k-1, 　　k=1,2,…,n/4
　　……

9. 多数桁計算の各種方式

　多数桁の計算方法は下記の様に多数ある。ここでは項目列挙に止める。
　(1) 筆算方式、計算量はO(n²)；桁数が短い時に使用
　(2) 実FFT方式、計算量はO(n・log(n))；高速だがメモリを多数使用
　(3) 中国剰余定理 と筆算方式の組み合わせ、O(n^3/2)；1000桁前後が高速
　(4) 整数FFT又は整数FMT、計算量はO(n・log(n))；現計算機では実FFTの方が少し速い
　(5) 同上と中国剰余定理の組み合わせ、計算量はO(n・log(n))；メモリは実FFTより少ない
　(6) 実FFTと最終数段に Karatsuba法; Karatsuba法は使用メモリ量を減らすために使用
　　　Karatsuba法は1回適用すると演算量が1.5倍(筆算方式は2倍)増加(筆算方式は2倍)する
　(7) 2段階FFT(FMT)、計算量はO(n・log(n))；計算量を増加させないでメモリ使用量を減少
　　　させる利点があり、今回の世界記録の更新に利用
　(8) 複素FMTによる直接乗算、計算量はO(n・log(n))；実FMT(FFT)を経由しないで多数桁乗算
　　　が可能、このため並列化には実FFTより有利、メモリ削減にはこれを下位FMTとして、
　　　(7)の2段階FMTが必要。世界記録達成のあとに考案
　(9) 複数個の異なる巡回FMTを利用した分割多数桁乗算、計算量はO(n・log(n))；上位FMTに
　　　この巡回FMTを適用して多数桁の乗算を分割して行う。世界記録達成のあとに考案
　　　演算量の増加なしで、メモリを少なくする究極の計算方式と思われる。実用テスト未