円周率世界記録の概要
2002/12/20 後 保範 ( Ushiro Yasunori )
2003/05/11 5/16の7Fセミナーのため更新 後 保範
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. 2002年世界記録達成の経過
主計算と検証計算の両方とも、全体の計算効率を考えて最初に16進で計算し、計算
の正しさが確認された16進結果を10進に基底変換することで10進の結果を得る方法を
今回は採用した。計算で使用する基本演算法として、2000年6月に情報処理学会論文誌
で発表した、適応の範囲が広く計算効率も高い、新しく考案した
DRM法を、
主計算、検証計算、基底変換のそれぞれに適応した。16進でのπ計算結果の正しさは、
独立なπに関する２種類の公式に基づくそれぞれの計算結果を比較することで確認する
こととし、16進->10進への基底変換の正しさは、16進入力と同じ結果が、16進->10進、
10進->16進の連続する２つの基底変換により得られるかどうかで確認することとした。
16進でのπ計算は主計算、検証計算共に1,030,775,439,375桁で実施した。比較結果は
最後の桁まで一致した。その理由は、16進約100万桁を１単位とした浮動小数点表示で
各Arctanの値(総て1より小さい）を計算し、整数係数を掛けて合算し、π(>1)の値を
求めたため、πの結果として求めた値より、各Arctanの値の計算は結果的に約100万桁
多く求めたことになっていたためである。16進π値の公表値は75,439,375桁桁を除いた
1兆307億桁とした。
16進->10進への基底変換は1,241,177,304,180桁で実施した。更に10進->16進への基底
変換をおこない、元の16進πの値と比較し当該桁の一致を確認した。10進π値の公表は
一時、4,180桁を除いた1兆2411億7730万桁を考えていた。最終的に億の単位の1兆2411
億桁で公表することにした。
2.1 利用公式
今回は、最近20年間ほとんどの世界記録の計算に使用されていた、
算術幾何平均法 |
(AGM法、相加相乗平均化法、サラミン-ブレント法、ガウス-ルジャンドル法とも呼ばれる）
を採用しないで、逆正接関数(Arctan)による
級数展開法 | を使用した。
円周率の多数桁計算の原理を示すソースプログラムはπ計算プログラム | を参照。
2.1.1 主計算に使用した公式
1982年に高野喜久雄発見の下記公式を使用
π/4 = 12Arctan(1/49) + 32Arctan(1/57) - 5Arctan(1/239)
+ 12Arctan(1/110443)
2.1.2 検証計算に使用した公式
1896年F.C.M. stoemer発見の下記公式を使用
π/4 = 44Arctan(1/57) + 7Arctan(1/239) - 12Arctan(1/682)
+ 24Arctan(1/12943)
2.1.3 Arctan公式の級数展開式
Arctan(1/x) = 1/x - 1/(3x3) + 1/(5x5)
- 1/(7x7) + 1/(9x9) -+ . . .
2.2 計算経過
2.2.1 16進1兆桁の主計算
開始日時： 2002/9/3 22:09
終了日時： 2002/9/25 17:00
使用時間： 1,439,973 s = 400h 0m
中断回数： 5回
磁気記憶装置への入出力量： 395,000 GB
磁気記憶装置への入出力性能： 1.1 GB/s
磁気記憶装置への入出力時間： 100 h
2.2.2 16進1兆桁の検証計算
開始日時： 2002/9/26 16:46
終了日時： 2002/10/4 10:40
使用時間： 565,468 s = 157h 4m
中断回数： 4回
2.2.3 16進の比較検証
開始日時： 2002/10/4 12:47
終了日時： 2002/10/4 13:58
使用時間： 4275 s = 1h 11m
2.2.4 16進->10進1.2兆桁への変換
(1) 初回（結果不正）
開始日時： 2002/11/10 00:06
終了日時： 2002/11/11 18:27
使用時間： 83783 s = 23h 16m
中断回数： 3回
(2) 記録達成時
開始日時： 2002/11/15 13:41
終了日時： 2002/11/16 13:01
使用時間： 84019s = 23h 20m
中断回数： 0回
2.2.5 10進->16進1兆桁への逆変換
(1) 初回（結果不正）
開始日時： 2002/11/12 01:56
終了日時： 2002/11/13 13:04
使用時間： 76809 s = 21h 20m
中断回数：１回
(2) 記録達成時
開始日時： 2002/11/23 16:07
終了日時： 2002/11/24 13:39
使用時間： 77540 s = 21h 32m
中断回数： 0回
2.2.6 再変換結果の検証
開始日時： 2002/11/24 13:48
終了日時： 2002/11/24 15:07
使用時間： 4661 s = 1h 18m
2.3 計算結果
16進でのπ計算結果の正しさは、主計算と検証計算の２種類の公式に基づくそれぞれの
計算結果を比較することで項番2.2.3で2002/10/4に確認した。
16進->10進への基底変換の正しさは、16進入力と同じ結果が、16進->10進、10進->16進
の連続する２つの基底変換により得られるかどうかで2002/11/14に確認した。
2.3.1 16進1兆桁の正確性確認結果
前記の2個の異なるArctan公式で16進数表示の1,030,775,439,375桁を計算。
両公式の計算結果は最後の桁まで一致。最後の桁まで一致する理由は、約100万桁単位の
で指数が１変わる浮動小数点表示を採用し、1以下の値を計算し定数を掛けて加える操作
で16進(3.243F…)以下の円周率を計算しているため。要するに途中計算は自動的に約100万
桁余分に計算している。
計算桁(1,030,775,439,375桁)の最後の50桁は下記の通り。
B556AACBA8 4A22B25C79 26AB3CD02F 532B7E71EE 4E10E249E5
計算宣言桁数は億より下位の桁を除き1兆307億桁とした。この最後の50桁は下記の通り。
28C7C6A1DC 75E2D2B765 2D23340DAE 5573D62B57 1BAF9CC340
2.3.2 10進変換の正確性確認結果
1,030,775,439,375桁の16進数を入力し、10進表示の1,241,177,304,180桁に変換。
検証は計算した10進数を別プログラムで16進数に逆変換して最後の桁までの一致を確認。
計算桁(1,241,177,304,180)の最後の50桁は下記の通り。
8508566667 3937845843 1454411298 9344841413 3413557641
計算宣言桁数は憶より下位の桁を除き1兆2411億桁とした。この最後の50桁は下記の通り。
5591198918 2262704528 2696896699 2856706487 3410311045
2.4 プログラム規模
プログラム作成言語はFORTRANで一部にCを利用した。
並列化のためのノード間通信ライブラリーはMPI及びリモートDMAを利用した。
プログラムはArctanによる16進π計算、計算した16進πを10進πへの変換及び逆に10進π
16進πに検証のために変換する３部に別れている。異なるArctanによる16進πの計算は
パラメータの変更でおこなった。また、このプログラムステップ数には前年におこなった
ラマヌジャン方式の16進π計算プログラムは含まれていない。ステップ数はコメント文も
含めた行数である。ステップ数がこれほど多くなった原因は、極限に近いチューニング
及びメモリーの使用量を少なくするための各種工夫のためである。
16進π計算(初回+検証) ： 5万4600行
10進数πへの変換 ： 1万4600行
16進数πへの逆変換 ： 1万行
合計 ： 7万9200行
3. 2001年世界記録チャレンジ過程
2001年に円周率1兆桁の世界記録にSR8000/MPP 128ノードを使用して挑戦しました。
倍精度浮動小数点の精度保持に関するプログラムミスで3597億余桁以降で、2種類の公式
の計算結果が異なり、2001年の挑戦を諦めました。
2001年に世界記録に挑戦した過程は2001年世界記録挑戦過程 | を参照下さい。
4. 記録挑戦への5年間の経過
4.1 DRM法の考案
1998年11月9日、京大数理解析研究所の研究集会「数値計算における前処理の研究」で
「逆数型無限級数のn桁計算の演算量を削減する前処理方式」と題して発表。
同年末に金田教授と本方式による円周率１兆桁の可能性に付いて会話。
金田教授から本方式は多数桁関数計算の基本システムとなるため、円周率の計算の準備の
前に各種応用方法を検討し、論文化した方が良いと御教授いただいた。
4.2 DRM法の論文化
1998年末から東大で、金田教授、高橋大介助手(現筑波大学講師)と一緒に応用方法の検討
及び京都大学数理解析研究所講究録及び情報処理学会論文誌の作成を開始。
ここで、連分数及び基底変換への応用をまとめた。また、入力値が出力値と同じ多数桁の
場合にも加法定理等と合せて使用する事で、多数桁関数の計算がシステムとして完成。これは、
この論文化の過程で得られた大きな成果である。このとき、英語論文も作成したが、これは
世界記録達成の結果も追加して投稿する事になった。今後この作業に着手する。
4.3 PC用π計算プログラムの開発
1999年夏頃より、ラマヌジャン公式で円周率１兆桁計算に向けた原理プログラムの作成を開始。
最初は10進でDISKをワークとして使用しないプログラムを作成。多数桁の乗算は実FFT+Karatsuba
法や整数FFTの重ね合せなど多数の試作プログラムを作成した。その結果、やはりメモリがきつい
事が判明した。そこで、16進でπを計算し、16進から10進に変換する事にした。
SR8000/MPPは倍精度浮動小数点の乗算が演算器内では、4倍精度浮動小数点として表わされる
(浮動小数点レジスターには取り出せない) ことを上手に利用すると、A,B,Cを倍精度浮動小数点として
mod(A*B,C)が高速に計算可能である。但し、FORTRANやCではこの機能が上手に利用できないため
アセンブラー化する必要がある。そのため、整数FFTの重ね合せで多数桁乗算をする方式が速度及び
メモリー使用面で有効だと考えた。しかし、SR8000でしか利用出来ない事及び一部のアセンブラー化
が必要で、PCを使用した世界記録用の円周率原理プログラムの開発には向かないと判断し断念した。
そこで、メモリを有効に利用し演算効率を落とさない方式として2段階FFTを考案した。
これで、円周率１兆桁ができる見通しはできたが、やはり現実に１兆桁を計算するプログラムの
ベースプログラムにするためには多くの困難が合った。
4.4 2001年の挑戦
2001年春からSR8000/MPPの128ノードで1兆桁計算可能な並列プログラムの作成を開始した。分担
し、東大チームが全体の並列化を分担した。金田教授は全体の纏め及びチューニングを担当された。
私(後)は１ノード用のプログラムが、並列化のベースプログラムとして完全ではないため、その修正、
並列化プログラムの調整及び並列化後のプログラムテストに多くの時間を割いた。そのため、並列化
チューニングもほとんど担当していない。私の所属する(HPC)は並列化とチューニングは専門部署なの
だが。言い換えると、１兆桁計算は生易しいものではなく、一人で多数の事が兼任できるほどの時間
的余裕はなかった。途中計算は例え、計算機が計算誤りをしてもほとんど摘出出来る仕組みを入れた
はずであったが。結果的には2公式での計算結果が3597億桁以降で不一致となり、計算機利用の都合
で記録への挑戦は来年度へ持ち越しとなった。
4.5 2002年の挑戦
昨年の様にSR8000/MPPの128ノードを専有して利用するのではなく、64ノードを使用することにした。
このため、時間は2〜3倍必要となるがArctan公式を利用する事に変更した。Arctan公式の１ノード
プログラムは昨年末に完成していた。6月までは、16進10進変換及び10進16進の逆変換のPC用プログラム
の作成をおこなった。本格的にプロジェクトとしてArctan公式の並列化に取り組んだのは7月からである。
もちろん、多数桁乗算等の並列計算の基本部品はそのまま使用できるが、64ノード計算でメモリ利用量
が半減したため、Diskの読み書きが昨年版より桁違いに増えた。Diskの並列I/O（16台並列)の性能を最大
限引き出すための工夫(プログラム上だけでなく、OSやDiskコントローラの最適化も含めて)をした。
その結果、16台のDiskで最高性能の90%を超す1.1GB/sの読み書きが出来る並列I/Oが達成できた。
これでもDiskのI/O時間は、初回の16進1兆桁の計算時間400時間のうち100時間に相当する。
何しろ、全体ではTバイト(1012)のさらに上の単位のFバイト(1015)のI/Oなので。
今更ながら、計算機もDiskもよく耐えてくれたと感心しています。とても普通のサーバでは無理でしょう。
F1マシンのレースを10人のドライバーで2週近く連続でおこなっているわけですから。ピットには時々入る
にしてもなかなか耐えられるF1マシンはないでしょう。
何しろ、チューニングによる性能向上はプログラム全体で4倍を超えました。最初から連続アドレス化
や並列オプションはほとんどもれなく入れおり、一般のプログラムよりは相当注意深く作成していたもの
からの向上比率ですから凄いです。何しろ、初期計算の途中でも中断時点でチューニングをしていました。
FORTRANのシフト関数などは、仕様的に遅いので日立のコンパイラー部隊が専用のルーチンを作成して
います。ノード内のチューニングは本プロジェクトのオーナーでリーダの金田教授がほとんどされました。