世界記録のための工夫点

　　　　　　　　2002/12/20　後保範 ( Ushiro Yasunori )
　　　　　　　　2003/05/12　5/16の7Fセミナーのため更新　　後保範

1. π計算の全体使用領域削減

　1兆バイトのメモリーで1.2兆桁を計算するのはやはり大変です。
メモリー効率を従来の6倍にするわけですから。ちなみに、従来と今回の比較は以下です。
(1) 従来の計算 (メモリー量/計算桁数、値が小さいほどメモリー効率が良い)
　1997年(SR2201)： 2620億バイト/515億桁 = 5
　1999年(SR8000)：　　1兆バイト/0.2兆桁 = 5 ：前回の世界記録
(2) 今回の計算
　2002年(SR8000/MPP)： 1兆バイト/1.2兆桁 = 0.8

1.1 DRM法を利用してπ計算公式を変更

　今回使用したArctan公式による計算は1983年に私と金田教授が1000万桁の計算をしてから
今回の世界記録まで記録達成に使用されていません。その理由は、Arctan公式の演算量が
O(n²)とAGM方式に比較してとても遅いと思われていたためです。ちなみに、1983年の1000万桁
は、今では信じられないO(n²)のアルゴリズムのままで計算しています。
これが、DRM法 (Divide and Rationalize Method、分割有理数化法)の使用でAGM法とほぼ同等
のO(n･(log(n))^2～3)で可能になりました。1兆桁計算するとO(n²)のアルゴリズムに比較し約1000万倍
高速になります。
　一方、AGM法の欠点は桁数全体を保存する複数のワークエリアと全桁の乗算及び平方根が毎回
(約log₂(n)回)必要な事です。このため、計算桁数の何倍ものメモリが必要です。
Arctan公式は工夫によりこれらのワークエリアが不要となります。π計算では複数のArctanを
利用しますので、この値の保存は必要ですが。
　DRM法を下記に簡単に紹介します。級数関数だけでなく連分数や基底変換にも利用可能です。
更に、加法定理等と結合する事により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
　　……

1.2 分割処理計算

　DRM法によるトーナメント通分処理は、最終的に１個の分数になるまで計算し、分子、分母
共に必要桁数以上になれば必要桁数で打ち切るのが演算量を最小に出来る。
しかし、今回はメモリの使用量をできるだけ少なくするのが大前提であったため、演算量の増大
には多少目をつぶり、通分処理は必要桁数の1/8で打ち切り、正確な分数から桁数の長い値を
計算する分割方式を採用した。このため、通分化した級数の計算は後の項から行なった。
通分処理を必要桁数の1/16で打ち切れば、メモリ面は楽であったが演算量が2倍以上になるため
ぎりぎりメモリに収まる1/8で打ち切る方式を採用した。これによる演算量の増加は1.3倍と推定
この方式の採用により、１兆桁同士の乗算と除算は不要となった。　

1.3 DISKの活用

　16進1兆桁の円周率の計算で、SR8000/MPPの128ノード(2Tバイト)が利用できればArctan公式では
ワーク領域としてのDiskの利用はArctan(1/x)の値を保存しておくだけであり、合計数Tバイトの
読み書きで済、計算全体の性能には影響しない。
一方、2002年の世界記録達成はSR8000/MPPの64ノード(1Tバイト)で行なうため、通分の分割処理の
ために大量のDiskの読み書きが必要となった。その量は数百Ｔバイトの単位で、Diskの並列I/Oの
性能が全体の性能を大きく左右する。最終的に16台のDiskの並列I/Oでカタログ値の90%を超える
1.1Gバイト/sを達成できた。

2. 多数桁乗算の使用領域削減

　FFT(高速フーリエ変換)とFMT(高速剰余変換)の関係は FMT概要と FMT原理を参照下さい。
　ここでは、FFTとFMTは一部の計算順序を除いて、結果的に同じ意味で記載しています。
　多数桁乗算は適当な桁数を１要素として、畳込み演算結果を桁合わせしてたものです。

2.1 正巡回FFTと負巡回FFT

　　FFTを利用すると多数桁の乗算が高速に計算できます。しかし、この時多くの計算領域
　が必要となります。このため、入力の2倍の桁数になる多数桁の乗算を、半分の領域で
　実行可能な正巡回FFTによる乗算と、負巡回FFTによる乗算に分けて計算します。
　　計算方法はFMTの応用の｢3. 正巡回FMTと負巡回FMTによる畳込み演算] を参照下さい。

2.2 整数FFT

　　計算方法はFMTの応用の｢2. 整数FMTと複数個のFMTの重ね合せ] を参照下さい。

2.3 2段階FFTによる多数桁乗算

　　計算方法はFMTの応用の｢4. 2段階FMTを使用した畳込み演算] を参照下さい。

2.4 DISKを利用し乗算

3. 多数桁の記憶と計算方法

　世界記録のπ計算では、8バイト整数に2進60桁(16進計算)又は10進16桁(10進計算)を記憶
　した。多数桁乗算(除算)の計算は、下位FFT(FMT)使用時は4分割して、倍精度浮動小数点
　(8バイト)に2進15桁又は10進4桁を記憶して計算した。上位FFT(FMT)使用時及び多数桁の
　加減算処理では8バイト整数に2進60桁又は10進16桁を記憶して計算した。

3.1 多数桁記憶の概要

　　SR8000/MPPでの計算では、多数桁の記憶はノード内で桁数が小のときはA-Typeデータで、
　桁数が大のときはB-Typeデータの形で記憶し、多数桁をノードにまたがり記憶するときは
　C-Typeデータの形で記憶した。いずれのTypeのデータも、固定小数点の形式と浮動小数点
　の2種類を用意し、必要に応じて利用した。

3.2 A-Typeデータ

　(1) 固定小数点
　　　8バイト整数(I*8)に2進60桁又は10進16桁を詰めて記憶した。記憶要素数(桁数)nは別に
　　記憶した。下記のようにn要素の多数桁をEのn-1次の多項式fで示すと、fの係数 a_kをa₁
　　から順にa_nまで記憶する。
　　　　f = a₁ E^n-1 + …… + a_n-1 E + a_n
　　　　　　E = 2⁶⁰ or 10¹⁶ ,　　| a_k | < E/2 , 　　 k = 1, 2, …… , n　　　　　　- - - (1)
　(2) 浮動小数点
　　　8バイト整数(I*8)の指数部とn個のI*8要素の仮数部で構成する。
　　小数点は仮数部の先頭にあり、指数の単位は(1)式のEに対応する。

3.3 B-Typeデータ

　　本Typeのデータの乗算は2段階FFT(FMT)を利用して実施される。2段階FFTの上位FFTの
　要素数をNとし、１要素は8バイト整数(I*8)n個で構成する。各要素の先頭の１個(I*8)はダミー
　として常にゼロにしておく。
　(1) 固定小数点
　　　A-Typeデータと同様に8バイト整数(I*8,ワード)に2進60桁又は10進16桁を詰めて記憶した。
　　下記のようにN要素(N・nワード)の多数桁をEのN・(n-1)-1次の多項式gで示すと、gの係数b_k
　　をb₁から順にb_N・(n-1)まで記憶する。
　　　　g = b₁ E^N・(n-1)-1 + …… + b_n E^{(N-1)・(n-1)-1} …… + b_N・(n-1)-1 E + b_N・(n-1)
　　　　　　E = 2⁶⁰ or 10¹⁶ ,　　| b_k | < E/2 , 　　 k = 1, 2, …… , N・(n-1)　　　- - - (2)
　(2) 浮動小数点
　　　8バイト整数(I*8)の指数部とN・n個のI*8要素（ワード)の仮数部で構成する。
　　小数点は仮数部の先頭にあり、指数の単位は(2)式のE^n-1に対応する。即ち、指数部の1
　　は仮数部の1要素(nワード、2進60(n-1)桁、10進16(n-1)桁)に対応する。
　　　仮数部は固定小数点と同様に１要素は8バイト整数(I*8)n個で、全体としてN要素で構成
　　する。仮数部の各要素の先頭の１個(I*8)はダミーとして常にゼロにしておく。

3.4 C-Typeデータ

　　本TypeのデータはB-Typeデータと基本的には同じであるが、複数ノードで１つの値を記憶
　する。ノード間のデータの配置方法はブロック・サイクリック分割である。
　(1) 固定小数点
　　　構成方法はB-Typeデータと同じで、要素数をNとし、１要素は8バイト整数(I*8)n個で構成
　　する。各要素の先頭の１個(I*8)はダミーとして常にゼロにしておく。ブロック・サイクリック
　　分割の単位は１要素、つまり8バイト整数(I*8)n個である。
　(2) 浮動小数点
　　　構成方法はB-Typeデータの浮動小数点と同様である。仮数部のブロック・サイクリック
　　分割ほ方法はC-Typeデータの固定小数点と同じである。一方、指数部は各ノード毎に保持
　　する。ただし、有効な指数部は先頭ノードのものだけである。

4. 多数桁計算での有効ビットの確保方法

4.1 2段階FFTの利用

　　入力A,Bが5000億桁で出力Cが倍の1兆桁(共に10進換算)の多数桁乗算で、入力値A,Bは保存しない
　ケースで3種類の計算方式のメモリ使用量を以下に比較する。
　(1) 2段階FMTによる多数桁乗算(整数8バイトに2進60桁詰めが可能)
　　　FFT(A),FFT(B)及び計算結果領域：　1.0兆バイト
　　　AとBの正巡回乗算結果保存領域：　0.25兆バイト
　　　下位FFT及び並列化通信領域　　：　0.05兆バイト
　　　　　　　　合　　　計　　　　　　　　　　：　1.3兆バイト
　(2) 実FFTによる多数桁乗算(実数8バイトに2進9桁詰めが限界)
　　　FFT(A),FFT(B)及び計算結果領域：　6.0兆バイト
　　　FFT用ワーク及び計算用テーブル：　3.0兆バイト
　　　　　　　　合　　　計　　　　　　　　　　：　9.0兆バイト
　(3) 実FFTとKaratsuba法の組み合わせで多数桁乗算
　　この方法はKaratsuba法の適用回数により、使用メモリ量と演算量が異なる。
　　Karatsuba法の適用回数を増やすと使用メモリ量は減少するが、演算量は増加する。
　　　Karatsuba法適用　FFTメモリ　ワークメモリ　合計メモリ量　演算量比
　　　　　１回適用　　　4.5兆バイト + 1.0兆バイト = 5.5兆バイト　　 1.5倍
　　　　　２回適用　　　2.3兆バイト + 1.5兆バイト = 3.8兆バイト　　 2.3倍
　　　　　３回適用　　　1.2兆バイト + 1.7兆バイト = 2.9兆バイト　　 3.4倍
　　　　　４回適用　　　0.6兆バイト + 1.8兆バイト = 2.4兆バイト　　 5.1倍