こんにちは、年齢 3X 歳の peria です。

6月に @square1001 さんの「超高速！多倍長整数の計算手法」シリーズ(前編 / 後編) が投稿された際、Constants コンテスト で数年前に得点していたことを思い出しつつ、今だともう少し書けるかな、というノリだけで色々コーディングしていたら総合2位まで来たのでザックリ何したかを書いていこうと思います。

注意点としては、あくまでコンテスト用に書いたコードなので全てを公開するわけではありません。また、実際に私が選んだアルゴリズムを曖昧にする意味でも、見えていない最適化があるかもしれないという意味でも、いくつかの選択肢を書いていることが多々あると思います。また、あくまで上記シリーズに付け加える二次創作・アフターストーリー的なものなので「超高速！多倍長整数の計算手法」シリーズは読破しているものとして書いていきます。

保存形式の工夫

そもそものスタート地点であるデータ保存形式ですが、2-1-2. 大きな数はどうやって管理する？に「どちらにせよ、1 桁ごとの数字の列として管理しなければならない」と書かれていますがそんなことはありません。後の章でも書かれている通り、計算を間違わない範囲なら詰め込むことができますし、多くの桁を詰め込む方が速くなります。重要なのはどの程度までなら詰め込んで大丈夫なのか(計算アルゴリズムとその実装次第ではあります)、プログラムが不要に複雑にならないか、辺りを見据えながら決めることです。

ちなみに私は int16 に整数 4 桁ずつを詰めた 10000 進数で実装しています。途中までは int64 に 16 桁ずつ詰めていましたが後述する DRM では掛け算(詰め直しが必要になる)と足し算/引き算(詰め直しの必要がない)が交互に出てくる、という事情から変更しました。また、負の数は一切扱わないようにするので符号情報は持ちません。あくまで目的が「定数を最大 $2^{25}$ 桁計算すること」なので負の数が出そうな部分は全て呼び出し側で回避しています。y-cruncher もそうしてますし、多分問題ないはずです。
あとはメモリ管理をサボるために vector<int16> を継承し、以下のようなmethodを生やしています。

class Integer : public std::vector<int16> {
 public:
  static constexpr int64 kBase = 10000;
  static constexpr int64 kDigsPerElem = 4;
  Integer();
  Integer(int16);
  Integer(int64);
  Integer& operator=(int16);
  Integer& operator=(int64);
  Integer& operator+=(const Integer&);
  Integer& operator-=(const Integer&);  // 引く数の方が小さいことが前提
  Integer& operator*=(const Integer&);
  Integer& operator*=(int64);
  Integer& half();  // `x /= 2` を簡略化するための method
  Integer& operator>>=(int);  // bit ではなく Limb 単位で桁上げ下げする
  Integer& operator<<=(int);
};
void Print(const Integer&);  // 0 詰めして標準出力に表示する
void InvSqrt(int64 x, int64 size, Integer& out);  // sqrt(1/x) を size 桁(Limb 単位)求める
void Inverse(const Integer& a, int64 size, Integer& out);  // 1/a を以下同文 

ポイントとしては operator/= がありません。当初 operator/=(int64) は用意していましたが呼び出し元が x /= 2; しか無かったので half() に変えました。また operator/=(const Integer&) については逆数計算をする関数 Inverse を作ったので必要ありません。

足し算、引き算の工夫

入力となる数がちゃんとしている前提であれば、繰り上がりや繰り下がりが 1 より大きくなることはありません。また、短い方の入力が終わった後は繰り上がりの数しか見る必要がありませんし、その途中で繰り上がりが0になればそれ以上は計算をする必要すらありません。コピーするだけです。

筆算での掛け算の工夫

次はパフォーマンス上最重要な掛け算の工夫、なのですが中身の実装が異なるものを長さに応じていくつか使い分ける必要があるのでまずは基本的な筆算から。

私の場合は int16 に 4 桁ずつ詰めてるので 1 つずつの積は 8 桁以下ですし、畳み込みをして正規化していくことを考えても、畳み込み結果は $10^8\cdot2^{23} < 2^{63}$ と int64 の配列に入れることで素直に計算できます。

FMT を使った掛け算の工夫

FMT は掛け算を $O(N\log N)$ の計算量で計算できる、現在知られている中で最速の方法です。ただ、その事実やDFTをFFTにする流れを書いた記事は多々あるものの、実際にそれをどのように適用していくかという点まで言及したものはあまり見ません。今回はココが 1 つのメインテーマなので少し深堀りしましょう。

FMT の高速化

FMT の実装ではそこそこ大きな選択が 2 つスタート地点にあります。それは

• 計算の中身として複素数を扱う FFT を使うのか、剰余環を扱う NTT を使うのか
• メインアルゴリズムを Cooley-Tukey のアルゴリズム にするか Stockham のアルゴリズム にするか

です。私の現在の実装では手元で色々実験をした結果として Cooley-Tukey の FFT を利用していますが、何か見逃している条件や実装次第で他の選択肢の方がより速くなる可能性ももちろんあります。

Cooley-Tukey vs Stockham

人の名前ばかりで何の違いなのか分かりにくいですね。
Cooley-Tukey アルゴリズムはバタフライ演算を行った変数にそのまま値を書き込む in-place 方式の演算をするため最終的な結果がビット逆順に並ぶ、もしくは入力値をビット逆順に並べておく必要があるアルゴリズムです。FFT アルゴリズムの説明で出てくるバタフライ演算をそのまま実装する形ですね。FFT ライブラリとして有名な FFTW はこの方式で実装されているらしいです。今回の畳み込み演算への応用だけを考えた場合、本来の FFT なら周波数データが並んでいる変換後のデータが順序よく並ぶ必要はありません。なので順変換と逆変換が別実装になる実装コストは上がりますが、並べ替えを行う処理を省いて高速化することが期待されます。
Stockham アルゴリズムは self-sorted アルゴリズムとも呼ばれるように、入出力データがちゃんと並んだままになってるアルゴリズムです。もちろん計算の途中もインデックス順に並んでいるためベクトル化しやすいメリットがあります。ただしその引き換えとしてデータと同じ大きさの作業用メモリが必要になります。スパコンのランキングに使われる FFTE は Stockham 方式で実装されています。

4 基底 と 8 基底 FMT

DFT を FFT にしていく過程では $\sum$ を 2 つに分解していく、2 基底 FFT がよく紹介されています。この分割は他にも 3 基底や 4 基底など自由に組み合わせることができますが、今回は長さ $n=2^k$ の FMT だけを対象にしているので奇数基底は用いません。逆に 4 基底、8 基底などを組み合わせることで計算が高速化することが期待されます。概要としての高速化の理由は

1. 同時アクセスするメモリが少ないので L2 キャッシュ以下のメモリにアクセスする頻度が減る
2. 回転子 $\omega^{jk}$ を掛ける回数が減るので演算が減る
3. FFT(複素数) の場合、虚数単位 $i$ や 45° の $\frac{1}{\sqrt{2}}(1+i)$ を掛ける演算が実装次第では簡略化できる

辺りがあります。現在の私の実装ではあまりパフォーマンスに差が無かったこととメンテナンス性から 4 基底 FFT を採用しています。

回転子のテーブル化

今回 FMT は畳み込み乗算を行うために導入されています。畳み込み乗算では一般的に

1. 長さ $n$ のデータ列 $\{a_j\}$ を変換して $\{A_k\}$ を得る
2. 長さ $n$ のデータ列 $\{b_j\}$ を変換して $\{B_k\}$ を得る
3. 項別乗算して長さ　$n$ のデータ列　$\{C_k\} = \{A_kB_k\}$ を得る
4. $\{C_k\}$ を逆変換して $\{c_j\}$ を得る

という過程を経ます。つまり $\omega$ が同じ値になる長さ $n$ の FMT を 3 回行うということです。そうなるとそれぞれの FMT で $\omega$ を求めていては計算の無駄になるので最初に $\{w^j\}$ のデータ列を作って再利用しましょう。
さらに、今回の計算対象では $n=2^k$ に限られていますし $n$ のパターンも少ないので作った $\{w^j\}$ のテーブルは再利用することで計算コストを減らすことができます。

実数限定 FFT の高速化

「多倍長整数の掛け算に FFT を利用する」という話の紹介内で最も省略されがちなのがこの複素数と実数の変換です。ご存知の通り FFT は複素数での計算で定義されてるものの、畳み込み乗算の入出力は実数です。実数データ $x$ を複素数 $x+0i$ と解釈して計算することももちろん可能ですが、半分のメモリ・計算が無駄に思えます。この無駄に対しては以下のような 2 つの解法が知られています。

a. 4n 乗根を用いた DWT
b. 実部と虚部をまとめた DFT

2 つの方法のそれぞれのメリット(＝もう一方にとってはデメリット)としては、計算量の違い、計算順序の違いと計算位置、があります。
まず計算量については見て明らかなように、複素数 1 つを掛けるだけの DWT の方が少なくなります。計算順序については、通常の FFT に追加する処理をどの段階で挟むかという話ですが、畳み込み処理を考えた場合

a. 4n 乗根の DWT : 追加処理→FFT→項別乗算→FFT→追加処理
b. 組み合わせ DFT : FFT→追加処理→項別乗算→追加処理→FFT

となっています。前提条件が限られていますが、Cooley-Tukey のアルゴリズムで並べ替えを行わない選択肢を取った場合、複素数化した DFT の方では対応する複素数の位置を見つける計算が複雑に、そして実際のアクセス位置も不連続になるので遅くなると思われます。
ちなみにいずれの方法であっても、 FFT の回転子テーブル $\{w^j\}$ を自然な順序(ビット逆順ではなく)で作っているとこちらの計算にも流用して高速化することができます。

NTT の高速化

NTT の最も大きな利点は整数化の誤差が無いので 1 要素に詰め込む桁数を多く取れることです。逆に言えば同じ桁数の演算をするときに使う要素数が少なくなるので演算全体が高速になる可能性が期待されます。さらに 64bit をフルに利用できる(浮動小数点数では正確に計算ができても有効数字は 53bit)こともあるので、$2^{25}$ 桁同士の掛け算をサポートできる範囲で考えると 1 要素に 6 桁詰めた計算が可能です。4 桁ずつ詰め込むことと比較するとザッと 50% の高速化 (33% の時間削減)が期待されます。
一方で中身の計算には uint64 同士の乗算、剰余計算が必須になるので遅くなりがちです。この "uint64 同士の乗算、剰余計算" に対する有名な対処法は Montgomery 演算です。詳細な紹介は省きますが剰余演算が実質的に必要なくなり、環境によっては 7 倍高速化もありえるっぽいです。

計算アルゴリズムの工夫

DRM

DRM (分割有理数化法) は以下のような数式を効率的に計算する方法です。ググる時は Binary Splitting 法という名前で探る方がよく当たると思います。それぞれの名前は数式の導出する過程が違う程度で、今回の範囲での実装上の違いはありません。

\frac{1}{x_0}\left(y_0+\frac{z_0}{x_1}\left(y_1+\frac{z_1}{x_2}(\cdots+\frac{z_{n-2}}{x_{n-1}}y_{n-1})
\cdots\right)\right)

導入を無視して計算法だけ言うと、上記の数式から取れる 3 数列を使って

\begin{pmatrix}Z&Y\\0&X\end{pmatrix} =
\begin{pmatrix}z_0&y_0\\0&x_0\end{pmatrix}
\begin{pmatrix}z_1&y_1\\0&x_1\end{pmatrix}\cdots
\begin{pmatrix}z_{n-1}&y_{n-1}\\0&x_{n-1}\end{pmatrix}

を分割統治で計算していき、最後に $X/Y$ を計算するアルゴリズムです。

√a の計算の工夫

√a ではなく 1/√a

超高速！多倍長整数の計算手法 で $\sqrt{2}$ はニュートン法を使うと $O(N\log N)$ で計算できる、という解説がなされていました。これ自体は間違いありません。ただ、このままでは掛け算より定数倍(私の実験では約 3 倍)遅い割り算がイテレーション中にあるのでもっと遅くなることが予想されます。ではそれをどのように速くするか、というと $\sqrt{2}=2 / \sqrt{2}$ で計算します。
計算対象になるグラフの式は $y=f(x)=\dfrac{1}{x^2}-a$ (ただし $a>0$)です。これの $x$ 軸との交点(の正の方)は $x=1/\sqrt{a}$ なのでニュートン法を適用すれば $1/\sqrt{a}$ を求められることが分かります。では改めてこのグラフ上の点 $(x_k,f(x_k))$ における接線の式を求めると

y-f(x_k) = f'(x_k)(x - x_k)\\
y-\frac{1}{x_k^2}+a = \frac{-2}{x_k^3}(x - x_k)

なのでニュートン法の漸化式は

-\frac{1}{x_k^2}+a = \frac{-2}{x_k^3}(x_{k+1} - x_k)\\
-x_k+ax_k^3=-2(x_{k+1}-x_k)\\
x_{k+1} = \frac{1}{2}x_k(3-ax_k^2)

となります。この漸化式の中では割り算は /= 2 しか使わないので half() で代用できますし、他は全て掛け算になるので速くなります。今回はコンテストでの計算対象にないので実装しませんが、$a$ が多倍長整数でも割り算の 2 倍まで遅くなることはありません。

ニュートン法での有効桁数

上記の $1/\sqrt{a}$ を求める漸化式は数式としては正しいのですが、計算上の重要な点を踏まえて書き直すと

x_{k+1} = x_k + \frac{1}{2}x_k(1-ax_k^2) \tag{n1}

となります。ニュートン法では 1 回のイテレーションで正確な桁数が約 2 倍に増えることが知られています。つまり $x_k$ が $n$ 桁正しい状態だとすると、各計算ステップで得られる/必要な桁数は以下の通りです。ここで入力値になる $a$ は多倍長数ではないものとします。

最後のステップでは $x_k$ の最下位桁と $\frac12x_k(1-ax_k^2)$ の最上位桁がほぼ同じ位置にあるため有効桁数 $n$ 桁の数同士の和の有効桁数が $2n$ 桁になります＆ほとんど足し算的作業は必要ありません。
細かな解析については各自計算してみて下さい。$x_k$ での誤差 $1/\sqrt{a}-x_k \leqq B^{-n}$ みたいに置いてから計算すると分かるかもしれません。また、今回は平方根の計算を例にしましたが、$a$ が多倍長数になる場合や逆数を求める割り算用ルーチンでも同じような形で計算ができます。

公式の選択

さてここまでの流れで演算そのものは高速化できましたが、数式の選択を間違うとまたパフォーマンスが落ちます。この章では私が値を出すのに使った、もしくは使っていた数式を紹介しておきます。DRM を使う場合、基本的には特定の桁数を計算するのに使う項数を見積もって少なくすむ数式を選ぶと良いです。

\begin{eqnarray}
\sqrt{2} &=& \sqrt{2}\\
\phi &=& \dfrac{1+\sqrt{5}}{2}\\
\end{eqnarray}

まずは簡単な 2 定数。これらは定義式通りの計算で、プログラムもほぼニュートン法を実装するだけの InvSqrt() でパフォーマンスが決まります。除算 Inverse() を実装しなくても計算できるので最初のターゲットに向いています。

\begin{eqnarray}
e &=& \left(1+\sum_{k=1}^{\infty}\dfrac{1}{2^kk}\right)^2 \tag{e1} \\
&=& \sum_{k=0}^{\infty}\dfrac{(3k)^2+1}{(3k)!} \tag{e2}
\end{eqnarray}

次は整数演算だけの DRM で計算できる $e$ です。√を使わないので DRM と Inverse() でパフォーマンスが決まる上に DRM もかんたんになります。2 式ありますが両方とも DRM の $z_k=1$ となるので $z_k$ に関わる計算をすっ飛ばせます。
ちなみに (e1) の式は $\exp(x)$ の Taylor 展開に $x=1/2$ を代入したものです。最後に 2 乗するコストはかかるものの $x=1$ を代入するより収束が速い(はず)です。 (e2) の式は Wikipedia から拝借しました。

\zeta(3) = \frac{1}{64}\sum_{k=0}^{\infty}\frac{(-1)^k(k!)^{10}(205k^2+250k+77)}{((2k+1)!)^5}

次は整数演算だけの DRM で計算できる定数その 2、 $\zeta(3)$ です。DRM の基本形を地道に計算します。式は Wikipedia から。

\begin{eqnarray}
\ln(2) &=& \dfrac{3}{4} + \dfrac{1}{8}\sum_{k=1}^{\infty}
 \begin{pmatrix}2n\\n\end{pmatrix} \dfrac{(-1)^k (5k+1)}{16^k k(k+\frac12)} \tag{l1}\\
&=& 18\coth^{-1}(26) - 2\coth^{-1}(4801) + 8\coth^{-1}(8749) \tag{l2}\\
&&\coth^{-1}(x) = \sum_{k=0}^{\infty} \frac{1}{(2k+1)x^{2k+1}}\\
\end{eqnarray}

(l1) 式は Wolfram から、(l2) 式は y-cruncherのページ から拝借しました。(l2) は DRM を複数回走らせる必要がある分少しコーディングの手間が増えます。また、その出力結果をどのように合わせて最終結果を得るかもまたパフォーマンスに影響します。

\pi = \left( \frac{12}{\sqrt{C}} \sum_{k=0}^{\infty} \frac{(-1)^k(6k)!(A+Bk)}{(3k)!(k!)^3C^k} \right)^{-1}\\
A=13591409,\ B=545140134,\ C=640320^3

円周率の公式は言わずと知れた Chudnovsky の公式です。これは DRM と√計算とが組み合わさった、ある意味総合問題になってます。ちなみに Salamin-Brent のアルゴリズムに比べて3倍くらい効率よく計算できるはずです。

おわりに

冒頭部分にも書いたとおり今回はコンテスト用に書いたコードを元にした解説なので全てを公開しているわけではありません。ココに書いていない工夫、問題点だけ挙げて解決策を書いていない問題、などありますがご了承下さい。コード全公開するとコンテストそのものがダメになってしまうので…。