calx の多倍長整数乗算はどのアルゴリズムを使っていますか？

calx はオペランドのワード数に応じて6段階のアルゴリズムを自動選択します。22ワード未満は Basecase (筆算)、22–139 は Karatsuba、140–199 は Toom-Cook-3、200–599 は Toom-Cook-4、600–1,049 は Toom-Cook-6、1,050ワード以上は NTT (Schönhage–Strassen) を使用します。

NTT (Schönhage–Strassen) 乗算とは何ですか？

数論変換 (Number Theoretic Transform) を用いた畳み込みベースの乗算アルゴリズムです。剰余環 Z/(B^F+1) 上で動作するため浮動小数点 FFT と異なり丸め誤差が一切発生せず、計算量は O(n log n log log n) です。calx では 1,400ワード (約26,600桁) 以上で使われます。

Karatsuba 法の計算量はいくつですか？

Karatsuba 法の計算量は O(n^1.585) です。入力を半分に分割し、4回の乗算を3回に削減する分割統治法です。1962年に Karatsuba が発見しました。

自乗 (squaring) は通常の乗算より速いですか？

はい。自乗では両オペランドが同一であるため、交差項 a_i·a_j と a_j·a_i を1回だけ計算して2倍することで、Basecase では乗算回数が約半分になり約30%高速化されます。Karatsuba や NTT にも自乗専用の最適化変種があります。

calx の乗算アルゴリズムの閾値はどのように決定されていますか？

AMD Ryzen Threadripper PRO 5995WX (Zen 3) 上で、予熱期間 (warm-up) と冷却期間 (cool-down) を含むベンチマークにより決定しています。将来的には利用者の環境で自動チューニングを実行し、閾値を最適化できる仕組みを提供する予定です。

Int 乗算アルゴリズム

概要

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

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

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

アルゴリズム選択と閾値

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

ワード数	アルゴリズム	計算量
$< 22$	Basecase (筆算)	$O(n^2)$
$22$–$139$	Karatsuba	$O(n^{1.585})$
$140$–$199$	Toom-Cook-3	$O(n^{1.465})$
$200$–$599$	Toom-Cook-4	$O(n^{1.404})$
$600$–$1\,049$	Toom-Cook-6	$O(n^{1.338})$
$\geq 1\,050$	NTT (Schönhage–Strassen)	$O(n \log n \log \log n)$

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

閾値は AMD Ryzen Threadripper PRO 5995WX (Zen 3) 上で、予熱期間 (warm-up) と冷却期間 (cool-down) を含むベンチマークにより決定している。 CPU キャッシュの状態やクロック周波数の安定化を考慮した計測であるが、異なるマイクロアーキテクチャ (Intel, ARM 等) では最適な閾値が変わる可能性がある。将来的には、利用者の環境で自動チューニングを実行し、閾値を最適化できる仕組みを提供する予定である。

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 回に減る代わりに加減算が増えるため、オーバーヘッドが乗算コストを上回る小さいサイズでは不利になる。 calx では 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 より大きいため、小さいサイズでは不利。calx では 140 ワード以上で切り替わる。

関連記事: 高速乗算アルゴリズム — Toom-Cook 法

Toom-Cook-4

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

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

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

実装上の最適化

calx の 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$–$599$ ワードの範囲で使われる。

関連記事: 高速乗算アルゴリズム — Toom-Cook 法

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})$) より漸近的に高速であり、$600$–$1{,}049$ ワードの範囲で使われる。 Toom-Cook-4 と同様に addmul_1 ベースの評価と小定数の正確除算による補間を採用している。

$1{,}050$ ワード以上では NTT が Toom-Cook-6 を上回るため、 Toom-Cook-6 の有効範囲は比較的狭い。

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)$

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

piecewise コストモデル

NTT の性能は変換長 $M$ と剰余環サイズ $F$ の選択に大きく依存する。 calx は 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 因子のビットシフトが整数ワード境界に乗るようにする。

$1{,}050$ ワード以上 (約 $43{,}000$ 桁) で NTT に切り替わる。 $1{,}050$–$1{,}249$ ワードでは Double FFT を試行し、 $1{,}250$ ワード以上では Prime NTT を直接使用する。数十万桁以上では他のアルゴリズムを大幅に凌駕する。

関連記事: 高速乗算アルゴリズム — FFT 乗算

不均衡乗算

オペランドのサイズが大きく異なる場合 (たとえば 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{ ワードのブロック})$$

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

関連記事: 不均衡 Toom-Cook 乗算

自乗の最適化

自乗 ($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 つの部分積がすべて自乗になるため、再帰的に対称性を活用できる。

自乗の閾値

自乗は一般の乗算より高速であるため、アルゴリズム切り替えの閾値が異なる。 calx では SQR_TOOMCOOK3_THRESHOLD = 200 であり、一般の乗算 (140) より Karatsuba の優位区間が広い。 Toom-Cook-4 専用の自乗は実装されていないため、 Toom-Cook-3 (sqr) から直接 NTT (sqr) に切り替わる。 $140$–$199$ ワードの自乗で Karatsuba (sqr) が使われ続けるのは、 Basecase 自乗の対称性削減 ($n^2/2$ 回の乗算) が Karatsuba の再帰にも波及するため、この区間では Toom-Cook-3 (sqr) より Karatsuba (sqr) のほうが依然として高速なためである。

ワード数	一般乗算	自乗
$< 22$	Basecase	Basecase (sqr)
$22$–$139$	Karatsuba	Karatsuba (sqr)
$140$–$199$	Toom-Cook-3	Karatsuba (sqr)
$200$–$2{,}599$	Toom-Cook-4	Karatsuba (sqr)
$200$–$2\,999$	Toom-Cook-4	Toom-Cook-3 (sqr) 以上
$\geq 3\,000$	NTT	NTT (sqr)

設計判断

6 段階選択の理由

GMP は Toom-Cook-6,7,8 やさらに高次の分割も実装しているが、 calx は 5 段階の分割統治 (Basecase, Karatsuba, Toom-3, Toom-4, Toom-6) + NTT の 6 段階としている。 Toom-Cook-5 は補間行列が複雑な割に Toom-Cook-6 に対する利点が薄いためスキップしている。

Toom-Cook-6 は $600$–$1{,}049$ ワードの区間で Toom-Cook-4 より 5〜10% 高速
NTT の閾値 ($1{,}050$ ワード $\approx$ $43{,}000$ 桁) を超えると NTT が圧倒的に速い
Toom-Cook-7,8 は有効区間がさらに狭くなるため実装の効果が薄い

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.