Int GCD アルゴリズム

概要

calx の GCD は Lehmer GCD + 再帰 HGCD (Half-GCD) の 2 段切り替えで実装されている。 小サイズでは Lehmer GCD、大サイズ (126 ワード以上) では HGCD (Möller 2008) に切り替える。

  • gcd(a, b) — Lehmer GCD + HGCD により最大公約数を計算する。 小〜中サイズで GMP を上回り、大サイズでも同等の性能。
  • extendedGcd(a, b, x, y) — 同じ HGCD ベースで $\gcd(a, b)$ と Bezout 係数 $x$, $y$ ($ax + by = \gcd(a, b)$) を同時に計算する。 SignedMpn 補因子追跡により高速。

計算量: Lehmer GCD は $O(n^2)$、HGCD は $O(M(n) \log n)$（$M(n)$ は乗算コスト）。

性能 (GMP mpz_gcd 比): 16K bit で 0.99x (同等)、4K bit で 1.31x。 extendedGcd も 16K bit で 0.99x。

API の使い方については Int API リファレンス — GCD / LCM を参照。

Lehmer GCD

Lehmer GCD はユークリッド互除法の改良版で、 多倍長整数の上位 1〜2 ワードのみを用いて複数ステップの商を一度に予測し、 多倍長除算の回数を大幅に削減する。

基本原理

多倍長整数 $a$, $b$ の上位数桁から近似的な商を求め、 その商が正しい限り多倍長演算なしでステップを進める。 近似が正しくなくなった時点で多倍長演算に戻る。

  1. $a$, $b$ の上位 1 ワード（64 bit）を取り出して $\hat{a}$, $\hat{b}$ とする。
  2. $\hat{a}$, $\hat{b}$ に対してユークリッド互除法を実行し、 変換行列 $\begin{pmatrix} \alpha & \beta \\ \gamma & \delta \end{pmatrix}$ を蓄積する。
  3. $\hat{b}$ が小さくなりすぎるか、近似の精度が保証できなくなったら停止する。
  4. 多倍長で $a' = \alpha \cdot a + \beta \cdot b$, $b' = \gamma \cdot a + \delta \cdot b$ を計算する。
  5. $b' = 0$ でなければ 1 に戻る。

計算量: $O(n^2)$（ユークリッドの互除法と同じオーダーだが、定数倍が改善される）。

HGCD (Half-GCD)

HGCD は分割統治法を用いて GCD の計算量を $O(M(n) \log n)$ に改善するアルゴリズムである （$M(n)$ は $n$ ワードの乗算コスト）。 calx では Möller (2008) の方式に基づき、126 ワード以上で HGCD に切り替える。

基本原理

$n$ ワードの GCD を $n/2$ ワードの HGCD 問題に帰着させ、再帰的に解く。 各再帰レベルで変換行列を求め、ボトムアップでマージする。

  1. 入力 $a$, $b$ の上位半分に対して HGCD を再帰的に呼ぶ。
  2. 得られた変換行列を $a$, $b$ 全体に適用する。
  3. 残りの部分に対してもう一度 HGCD を呼ぶ。
  4. 2 つの変換行列を合成して返す。

計算量: $O(M(n) \log n)$。NTT 乗算を使えば $O(n \log^2 n \log \log n)$。

閾値と性能

Lehmer GCD と HGCD の切り替え閾値は 126 ワード（約 2,400 桁）。 これはベンチマークにより決定された損益分岐点であり、 HGCD の再帰や行列乗算のオーバーヘッドが Lehmer の線形走査コストを下回るポイントである。

GMP mpz_gcd との性能比較（calx / GMP、1.0x で同等）:

  • 4K bit (1,233 桁): 1.31x（calx がやや遅い — Int 構築オーバーヘッド）
  • 16K bit (4,932 桁): 0.99x（同等）
  • 64K bit (19,728 桁): 1.27x

extendedGcd での HGCD

拡張 GCD (Bezout 係数の計算) でも同じ HGCD を使用する。 SignedMpn による補因子追跡と phys_swap によるバッファ整合性保証により、 ナイーブ $O(n^2)$ 比で 100 桁 4.6x、500 桁 33x、5000 桁 52x 高速。

関連記事: Binary GCD (Stein の算法)  |  多倍長整数の高度な演算

拡張ユークリッド互除法

拡張ユークリッド互除法 (Extended Euclidean Algorithm) は、 $\gcd(a, b)$ を求めると同時に、次の等式を満たす整数 $x$, $y$ (Bezout 係数) を計算する:

$$ax + by = \gcd(a, b)$$

この等式は Bezout の等式と呼ばれ、 $\gcd(a, b)$ が $a$ と $b$ の整数線形結合で表現できることを示す。

HGCD ベースの拡張 GCD

calx の拡張 GCD は gcd() と同じ Lehmer + HGCD アルゴリズムを使用し、 各ステップで Bezout 係数 (補因子) を同時に追跡する。 SignedMpn (符号付き多倍長配列) で補因子を管理し、 phys_swap によるバッファ整合性保証で効率的に実装されている。

アルゴリズム

初期値 $(r_0, x_0, y_0) = (a, 1, 0)$, $(r_1, x_1, y_1) = (b, 0, 1)$ として、 HGCD の各ステップで変換行列と同時に $x$, $y$ を更新する:

$$q_i = \lfloor r_{i-1} / r_i \rfloor$$

$$r_{i+1} = r_{i-1} - q_i \cdot r_i$$

$$x_{i+1} = x_{i-1} - q_i \cdot x_i$$

$$y_{i+1} = y_{i-1} - q_i \cdot y_i$$

最後の非零の $r_i$ が $\gcd(a, b)$ であり、対応する $x_i$, $y_i$ が Bezout 係数である。

応用

  • モジュラー逆元: $\gcd(a, m) = 1$ のとき、$ax \equiv 1 \pmod{m}$ を満たす $x$ は 拡張 GCD で求められる。inverseMod() は内部でこれを使用する。
  • 線形ディオファントス方程式: $ax + by = c$ の整数解は、$\gcd(a, b) \mid c$ のとき拡張 GCD を用いて構成できる。

計算量: $O(M(n) \log n)$（HGCD ベース）。 GMP mpz_gcdext 比では 16K bit で 0.99x（同等）。 ナイーブ拡張ユークリッド互除法 $O(n^2)$ 比では 5000 桁で 52x 高速。

LCM の計算

最小公倍数 (LCM: Least Common Multiple) は GCD を用いて次のように計算する:

$$\text{lcm}(a, b) = \frac{|a| \cdot |b|}{\gcd(a, b)}$$

ただし、実装では中間結果のオーバーフローを避けるため、 除算を先に行う:

$$\text{lcm}(a, b) = |a| \times \left( \frac{|b|}{\gcd(a, b)} \right)$$

$\gcd(a, b)$ は $b$ を割り切ることが保証されているため、 $|b| / \gcd(a, b)$ は整数除算で正確に計算できる。 先に除算することで、$|a| \times |b|$ という不必要に大きな中間結果の生成を回避し、 メモリ使用量と計算時間を削減する。

設計判断

Lehmer + HGCD の 2 段構成を選択した理由

GCD アルゴリズムにはいくつかの選択肢がある:

  • ユークリッドの互除法: $O(n^2)$、ただし毎回多倍長除算が必要で遅い。
  • Binary GCD: $O(n^2)$、シフトと減算のみで動作。ユークリッドより 30–50% 高速だが、大サイズでは HGCD に劣る。
  • Lehmer GCD: $O(n^2)$、上位ワードで複数ステップを一括予測。小〜中サイズで最速。
  • HGCD: $O(M(n) \log n)$、分割統治法。大サイズで漸近的に最速。

calx では Lehmer (小サイズ) + HGCD (大サイズ) の 2 段構成を採用した。 Binary GCD はシンプルだが Lehmer に劣り、ユークリッドの互除法は多倍長除算のコストが高い。 Lehmer + HGCD の組み合わせにより、 全サイズ範囲で GMP と同等以上の性能を達成している。

参考文献

  • Möller, N. (2008). "On Schönhage's algorithm and subquadratic integer GCD computation". Mathematics of Computation, 77(261), 589–607.
  • Stein, J. (1967). "Computational problems associated with Racah algebra". Journal of Computational Physics, 1(3), 397–405.
  • Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997. §4.5.2.