Int 乗算アルゴリズム

概要

Int の乗算は、オペランドのワード数に応じて最適なアルゴリズムが自動選択される。 ユーザーが明示的に選ぶ必要はない。 ワード数が小さい領域では定数倍の軽い Basecase を使い、 大きくなるにつれて漸近的に高速なアルゴリズムに切り替わる。

1 ワード = 64 bit であり、10 進で約 19 桁に相当する。 たとえば 100 桁の数は約 6 ワードなので Basecase が使われる。

API の使い方については Int API リファレンス — 算術演算子 を参照。

アルゴリズム選択と閾値

sangi は実測ベンチマークに基づいて閾値を決定している。 以下のテーブルは $a \times b$ において、$\min(\text{words}(a), \text{words}(b))$ に対するアルゴリズムの選択を示す。

Toom-Cook-3 の区間が $140$–$199$ と狭いのは、 sangi の Toom-Cook-4 実装が Toom-Cook-3 より定数倍が小さいためである。 Toom-Cook-3 は 3 分割で 5 回の乗算を行うが、 Toom-Cook-4 は 4 分割で 7 回の乗算を行う。 分割数あたりの乗算回数では Toom-Cook-4 のほうが効率的であり、 さらに sangi では addmul_1 ベースの固定サイズ評価と submul_1 ベースの補間によって Toom-Cook-4 のオーバーヘッドが大きく削減されている。

上記の閾値は Windows MSVC ビルド時の値である。 Linux gcc/clang ビルドでは NTT 経路が相対的に遅いため、 Toom-Cook-8 の上限を $2{,}999$ ワードまで引き上げ、$\geq 3{,}000$ ワードで NTT に切り替える。

Basecase (筆算乗算)

小学校で習う筆算と同じアルゴリズムを、64 ビットワード単位で行う。 $n$ ワード $\times$ $m$ ワードの乗算は、$n \times m$ 回のワード乗算と加算からなる。

$$c = a \times b \quad \Rightarrow \quad c_k = \sum_{i+j=k} a_i \cdot b_j + \text{carry}$$

計算量: $O(nm)$（$n$, $m$ はそれぞれのワード数）

実装では MULX 命令 (BMI2) を使い、 64×64→128 ビットの乗算を 1 命令で行う。 内部ループの addmul_1 は MASM 手書きアセンブリで最適化されており、 MULX + ADCX/ADOX を用いて キャリーチェーンを並列化している。

ワード数が 22 未満の場合、分割統治のオーバーヘッドが Basecase の $O(n^2)$ を上回るため、 Basecase が最速となる。

Karatsuba 法

1962 年に Karatsuba が発見した分割統治法。 $n$ 桁の乗算を、$n/2$ 桁の乗算 3 回 と加減算に帰着する。

$a$ と $b$ をそれぞれ上位半分と下位半分に分割する:

$$a = a_1 \cdot B + a_0, \quad b = b_1 \cdot B + b_0 \quad (B = 2^{32n})$$

通常の筆算では 4 回の半分サイズの乗算が必要だが、Karatsuba は次の 3 つの積から結果を構成する:

$$z_0 = a_0 \cdot b_0, \quad z_2 = a_1 \cdot b_1, \quad z_1 = (a_0 + a_1)(b_0 + b_1) - z_0 - z_2$$

$$a \cdot b = z_2 \cdot B^2 + z_1 \cdot B + z_0$$

計算量: $O(n^{\log_2 3}) = O(n^{1.585})$

乗算回数が 4 回から 3 回に減る代わりに加減算が増えるため、 オーバーヘッドが乗算コストを上回る小さいサイズでは不利になる。 sangi では 22 ワード以上で Karatsuba に切り替わる。

Toom-Cook-3

入力を 3 分割し、5 回の部分積から結果を再構成する。 Karatsuba の一般化であり、分割数 $k$ に対して $2k-1$ 回の乗算で済む。

入力を $a = a_2 x^2 + a_1 x + a_0$, $b = b_2 x^2 + b_1 x + b_0$ ($x = B^{n/3}$) と見なし、 多項式の乗算として処理する。 5 つの評価点 ($0, 1, -1, 2, \infty$) での値を計算し、 部分積から補間で係数を復元する。

計算量: $O(n^{\log_3 5}) = O(n^{1.465})$

補間ステップでの除算と加減算のオーバーヘッドが Karatsuba より大きいため、 小さいサイズでは不利。sangi では 140 ワード以上で切り替わる。

Toom-Cook-4

入力を 4 分割し、7 回の部分積から結果を再構成する。

入力を 3 次多項式として表現し、7 つの評価点での値を計算する。 評価点は $0, 1, -1, 2, -2, 3, \infty$ を使用する。

計算量: $O(n^{\log_4 7}) = O(n^{1.404})$

実装上の最適化

sangi の Toom-Cook-4 実装には以下の工夫がある:

• addmul_1 ベースの評価: 各評価点での多項式の値を、一時バッファと normalized_size の呼び出しを削減した 固定 $k$-limb の addmul_1 ベースの一括処理で計算する。
• submul_1 ベースの補間: 補間ステップの除算を submul_1 と divexact_by3、divexact_by5 等の 小定数による正確除算に帰着し、汎用除算を完全に回避する。

これらの最適化により Toom-Cook-3 に対する定数倍の優位が大きく、 $140$–$199$ ワードという狭い Toom-Cook-3 区間の原因となっている。 $200$–$749$ ワードの範囲で使われる。

Toom-Cook-6

入力を 6 分割し、11 回の部分積から結果を再構成する。 評価点は $0, \pm1, \pm2, \pm3, \pm4, \pm5, \infty$ の 12 点を使用する（$-0 = 0$ なので 11 回の乗算）。

計算量: $O(n^{\log_6 11}) = O(n^{1.338})$

Toom-Cook-4 ($O(n^{1.404})$) より漸近的に高速であり、$750$–$799$ ワードの範囲で使われる。 Toom-Cook-4 と同様に addmul_1 ベースの評価と小定数の正確除算による補間を採用している。

$800$ ワード以上では Toom-Cook-8 がさらに高速になるため、 Toom-Cook-6 の有効範囲は比較的狭い。

Toom-Cook-8

入力を 8 分割し、15 回の部分積から結果を再構成する。 評価点は $0, \pm1, \pm2, \pm3, \pm4, \pm5, \pm6, \pm7, \infty$ の 16 点を使用する ($-0 = 0$ なので 15 回の乗算)。

計算量: $O(n^{\log_8 15}) = O(n^{1.302})$

Toom-Cook-6 ($O(n^{1.338})$) より漸近的に高速で、Windows MSVC ビルドでは $800$–$1{,}249$ ワードの範囲で使われる。 Linux gcc/clang ビルドでは NTT 経路が相対的に遅いため、Toom-Cook-8 が $2{,}999$ ワードまで延長される。

sangi の Toom-Cook-8 実装には Toom-Cook-6 と同じ最適化に加え、 固定 $k$-limb の pre-pad 配置と addmul_1 評価により評価フェーズを 1 パスに圧縮する設計を採用している。 これにより当初 $n < 4{,}000$ で Toom-Cook-6 より遅かった実装が、$n \approx 800$ で Toom-Cook-6 と同等以上に達した。

不均衡乗算 (例: $\text{words}(a) \geq 2 \cdot \text{words}(b)$) では、 専用の Toom-8,4 (2:1 比に特化した deg-10 補間) が選択される場合がある。

NTT (Schönhage–Strassen)

数論変換 (Number Theoretic Transform) を用いた畳み込みベースの乗算。 1971 年に Schönhage と Strassen が発表した。 整数の乗算を「多項式の乗算 → 畳み込み定理による変換領域での点ごとの積 → 逆変換」で実現する。

通常の FFT は複素数上で動作するため丸め誤差が生じるが、 NTT は剰余環 $\mathbb{Z}/(B^F+1)$ 上で動作するため丸め誤差が一切発生しない。 これが整数乗算に NTT が適する理由である。

計算量: $O(n \log n \log \log n)$

剰余環 $\mathbb{Z}/(B^F+1)$

sangi の NTT は Fermat 数 $B^F + 1$ ($B = 2^{64}$) を法とする剰余環上で動作する。 $M = 2^l$ 個のピースに分割し、各ピースは $K$ ワード、 $F$ ワードの剰余環で演算する。 $\omega = B^{2F/M}$ を $M$ 次の原始単位根として使用し、 2 の冪乗のビットシフトだけで twiddle 因子の乗算が行える。

piecewise コストモデル

NTT の性能は変換長 $M$ と剰余環サイズ $F$ の選択に大きく依存する。 sangi は piecewise コストモデルを用いて最適なパラメータを選択する:

• 点ごと乗算コスト: $F$ のサイズに応じて使われる乗算アルゴリズム (Basecase/Karatsuba/Toom-Cook) の 実際のコストを区分的に見積もる。
• バタフライコスト: $M \cdot l$ 回のバタフライ演算 (加算・減算・ビットシフト mod $B^F+1$) のコスト。
• $l$ の初期推定 ($M \approx \sqrt{N}$) の前後で最良のパラメータを探索する。

その他の実装上の特徴

• thread_local バッファ: ワーキングメモリを thread_local で保持し、 呼び出しごとのメモリ確保・解放コストを回避する。
• F の量子化: $M > 128$ の場合、$F$ を $M/128$ の倍数に切り上げて twiddle 因子のビットシフトが整数ワード境界に乗るようにする。

Windows MSVC ビルドでは $1{,}250$ ワード以上 (約 $24{,}000$ 桁) で NTT に切り替わる。 Linux gcc/clang ビルドでは Toom-Cook-8 が $2{,}999$ まで延長されるため、$3{,}000$ ワード以上で NTT に切り替わる。 数十万桁以上では他のアルゴリズムを大幅に凌駕する。

不均衡乗算

オペランドのサイズが大きく異なる場合 (たとえば 1,000 ワード × 100 ワード) は、 大きいほうを小さいほうと同じサイズのブロックに分割し、 各ブロックの積を足し合わせる反復ブロック乗算を使う。

$$a \cdot b = \sum_{i=0}^{k-1} a_i \cdot b \cdot B^{i \cdot m} \quad (a_i \text{ は } m \text{ ワードのブロック})$$

sangi では $\text{words}(a) \geq 2 \cdot \text{words}(b)$ のとき不均衡乗算に切り替わる。 各 $a_i \cdot b$ は均衡な乗算として最適なアルゴリズムが適用される。 これにより、不均衡な入力でも効率的に処理できる。

自乗の最適化

自乗 ($a \times a$) は一般の乗算とは異なり、両オペランドが同一であるという対称性を利用できる。

Basecase 自乗

$a_i \cdot a_j$ と $a_j \cdot a_i$ は同じ値であるため、$i \neq j$ の項は 1 回だけ計算して 2 倍すればよい。 これにより乗算回数が約 $n^2/2$ に削減され、約 30% の高速化が得られる。

$$a^2 = \sum_i a_i^2 \cdot B^{2i} + 2 \sum_{i < j} a_i \cdot a_j \cdot B^{i+j}$$

高速アルゴリズムの自乗

Karatsuba・Toom-Cook・NTT にも自乗専用の変種がある。 Karatsuba 自乗では再帰呼び出しの 3 つの部分積がすべて自乗になるため、 再帰的に対称性を活用できる。

自乗の閾値

自乗は一般の乗算より高速であるため、アルゴリズム切り替えの閾値が異なる。 sangi では Basecase (sqr) の対称性削減 ($n^2/2$ 回の乗算) が効くため、 Karatsuba (sqr) の閾値が一般乗算 (22) より高い $74$ ワードに設定されている。 以降の Toom-Cook-3 (sqr) と Toom-Cook-4 (sqr) にも専用実装があり、 それぞれ $160$、$400$ ワードで切り替わる。

閾値は Windows MSVC ビルド時の値。 Linux gcc/clang ビルドでは Toom-Cook-4 (sqr) が $2{,}999$ まで延長され、$\geq 3{,}000$ ワードで NTT (sqr) に切り替わる。

設計判断

7 段階選択の理由

GMP は Toom-Cook-7 や非対称分割を多用しているが、 sangi は 6 段階の分割統治 (Basecase, Karatsuba, Toom-3, Toom-4, Toom-6, Toom-8) + NTT の 7 段階としている。 Toom-Cook-5 と Toom-Cook-7 は補間行列が複雑な割に Toom-Cook-6/8 に対する利点が薄いためスキップしている。

• Toom-Cook-6 は $750$–$799$ ワードの区間で Toom-Cook-4 より高速 (主に Toom-Cook-8 の再帰経路で効く)
• Toom-Cook-8 は $800$ ワード以上で Toom-Cook-6 より速く、固定 $k$-limb pre-pad と addmul_1 評価で評価フェーズを 1 パスに圧縮している
• NTT の閾値 ($1{,}250$ ワード $\approx$ $24{,}000$ 桁、Windows MSVC) を超えると NTT が圧倒的に速い
• Linux gcc/clang ビルドでは NTT 経路が相対的に遅いため、Toom-Cook-8 を $2{,}999$ ワードまで延長している

FFT ではなく NTT を使う理由

浮動小数点 FFT を用いた乗算は丸め誤差の管理が複雑であり、 桁数が増えると誤差が蓄積して正確性の保証が困難になる。 NTT は剰余環 $\mathbb{Z}/(B^F+1)$ 上で動作するため丸め誤差が発生せず、 結果の正確性が数学的に保証される。 多倍長整数の乗算にはこの性質が不可欠である。

参考文献

• Karatsuba, A. and Ofman, Yu. (1962). "Multiplication of Multidigit Numbers on Automata". Doklady Akademii Nauk SSSR, 145, 293–294.
• Toom, A. L. (1963). "The Complexity of a Scheme of Functional Elements Simulating the Multiplication of Integers". Doklady Akademii Nauk SSSR, 150 (3), 496–498.
• Schönhage, A. and Strassen, V. (1971). "Schnelle Multiplikation großer Zahlen". Computing, 7, 281–292.
• Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997.