raspberry piをクラスタ計算機として使ってみるために円周率計算のプログラムを作成していたところ、数学に現れる定数の計算速度を競うコンテストサイト（？）Constantを発見しました。 せっかくだし任意精度計算のためのプログラムから作ってみるか、ということで色々やってみたことを書きます。とはいっても円周率計算が目的なので、実装した機能はかなり限定的です。あくまでもやってみたことの記録であり、解説記事ではありません。また、必ずしも適切な方針をとっているとは限らないのでご注意ください。
任意精度整数（多倍長整数）
データの持ち方
16ビット整数を要素にもつVectorとして管理しています。16ビット整数型の範囲は-32,768～32,767ですが、表示のしやすさを考慮して10,000進数（基数10,000）として扱います。また、負の数は使わない（負の値が出ないように運用でカバーする）予定なので、負の数は実装しません。
演算
円周率や2の平方根計算には整数同士の割り算が登場しないため、今回は足し算、引き算、掛け算のみを対象として考えることにしました。
足し算・引き算
桁ごと（Vec<i16>の要素ごと）に足し算・引き算を行い、最後に繰り上げ・繰り下げを行えば良いです。
掛け算
掛け算を実現するアルゴリズムには選択肢が複数あります。代表的な例として、愚直な筆算（桁数を$N$として$O(N^2)$）, 分割統治法的に計算するカラツバ法（$O(N^{1.59})$）, カラツバ法よりも細かく分割することで計算量を削減したToom-Cook法（$O(N^{1.47})$）, そして$O(N\log N)$で計算可能な高速フーリエ変換（FFT）があります。FFTを用いた掛け算の仲間に、数論変換（NTT）と高速剰余変換（FMT）があります。
NTTは整数の世界で計算が完結するため数値誤差が発生しないですが、計算する桁数が大きくなるとオーバーフロー的現象が発生します。複数の$p$に対する計算結果からGarnerのアルゴリズムを用いて正しい計算結果を復元することでオーバーフロー的現象を回避することができます。
FFTは浮動小数計算を含むのでどうしても数値誤差の影響はありますが、NTTのように桁溢れ対策をしなくてもかなり大きな数まで計算できます。FFTの高速な実装が公開されているので、そういった点でも実装しやすいかと思います。
当初はNTTを用いて掛け算を実装していましたが、速度や桁溢れの点からFFTに乗り換えました。FFTを用いた掛け算の擬似コードは以下の通りです。ifftはfftの逆変換（$x=IFFT(FFT(x))$）なので、FFTとIFFTの違いは1のn乗根として$e^{\frac{2\pi i}{n}}$を使うか$e^{-\frac{2\pi i}{n}}$を使うか、また計算後に値を$n$で割るかどうか、です。
fn mul(a: &[i16], b: &[i16]) -> Vec<i16> {
let m = a.len() + b.len() - 1;
let mut n = 1;
while n < m {
n <<= 1;
}
let mut ar_ = vec![0.0; n]; // aの実部
let mut br_ = vec![0.0; n];
let mut ai_ = vec![0.0; n];
let mut bi_ = vec![0.0; n];
// a, bの値をar_, br_に書き込む
fft(ar_, ai_); // NTTやFMTに置き換え可能
fft(br_, bi_); // NTTやFMTに置き換え可能
let mut cr_ = vec![0.0; n];
let mut ci_ = vec![0.0; n];
// cr_[i], ci_[i]にar_[i], br_[i], ai_[i], bi_[i]の積を代入する。
// ここでいう積はA=ar_+i ai_, B=br_+i bi_とした時の複素数積である。
ifft(cr_, ci_); // NTTやFMTに置き換え可能
let mut c = vec![0_i16; n];
// c[i] = cr_[i]に繰り上げ処理+整数丸めした値
return c;
}実数列a, bを複素数列に変換する方法にも選択の自由度があります。例えば実数列の偶数部分は複素数の実部、奇数部分は虚部と対応させることで配列の長さを1/2にすることができます。しかし、私の環境では特に工夫しない場合がもっとも速く動作したので擬似コードでは実数列をそのまま実部に対応させる方法を記載しています。
表示
今回は基数が10,000なので、各桁の値を0埋めしながら表示すれば十分です。
任意精度浮動小数点
データの持ち方
多倍長整数$\times$基数$^{EXP}$として任意精度浮動小数点を表現しています。例えば$1.1212121212$であれば、多倍長整数部分は[1200, 1212, 1212, 1], 指数はEXP=-3となります。多倍長整数部分は、16ビット整数を要素にもつVectorとして実装しています。小数点以下の桁数を固定していないため、何桁まででも表現できます。一方で省メモリ化のための工夫は一切行っていないので改善の余地はあると思います。
演算
足し算・引き算
今回の実装では小数点以下の桁数が異なる数同士の演算、整数部分の桁数が異なる数同士の演算が発生します。そこで、小数点以下の桁数を揃えた状態で計算を行います。それ以外の計算手順は多倍長整数の場合と同様です。計算結果を小数点以下の桁数が長い方に揃えるため、指数部分は小さい方に揃えます。
掛け算
多倍長整数部分の掛け算と基数$^{EXP}$部分の掛け算に分割して考えます。基数$^{EXP}$部分の掛け算はEXPの足し算なので特にいうことはありません。多倍長整数部分の掛け算はFFTを用いた実装を使い回します。
逆数・平方根
逆数・平方根の計算にはニュートン法を使用します。割り算は逆数を用いて表現します。例えば$A$の逆数を計算したい場合、$f(x)=A-\frac{1}{x}$に対して$f(x)=0$を満たす$x$を求めれば良いことになります。平方根については、平方根の逆数（$\frac{1}{\sqrt{A}}$）を求める方が演算量が少ないです。逆数と同様に考えると、$f(x)=A-\frac{1}{x^2} = 0$の解を求めることで$\sqrt{A}=Ax$を求めることができます。
表示
小数点の位置だけ気をつければ、あとは多倍長整数の場合と同様です。
参考文献
- 超高速！多倍長整数の計算手法【前編：大きな数の四則計算を圧倒的な速度で！】
- 超高速！多倍長整数の計算手法【後編：N! の計算から円周率 100 万桁の挑戦まで】
- 高速フーリエ変換・数論変換を改めて理解しようという話
- 【４基底の Stockham FFT】
- 円周率.jp
- 超高速！！多倍長整数の計算手法
- FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法
この記事は役に立ちましたか？
もし参考になりましたら、下記のボタンで教えてください。
コメント