Int 除算アルゴリズム

概要

Int の除算は、除数のワード数に応じて最適なアルゴリズムが自動選択される。ユーザーが明示的に選ぶ必要はない。除数が小さい場合は定数倍の軽い Schoolbook 除算を使い、大きくなるにつれて高速乗算を活用する Burnikel-Ziegler に切り替わる。被除数が除数に対して極端に大きい不均衡なケースでは Newton 逆数反復が使われる。

1 ワード = 64 bit であり、10 進で約 19 桁に相当する。

API の使い方については Int API リファレンス — 算術演算子を参照。また、IntOps を通じて floorDiv、ceilDiv、divExact といった特殊な除算バリエーションも提供されている。

アルゴリズム選択と閾値

calx は実測ベンチマークに基づいて閾値を決定している。以下のテーブルは $a \div b$ における除数 $b$ のサイズと被除数 $a$ との関係に対するアルゴリズムの選択を示す。

ケース	アルゴリズム	計算量
1 ワード $\div$ 1 ワード	ネイティブ 64 ビット除算	$O(1)$
除数 $< 64$ ワード	Schoolbook (Knuth Algorithm D)	$O(nm)$
除数 $\geq 64$ ワード	Burnikel-Ziegler	$O(M(n) \log n)$
不均衡 (被除数 $\gg$ 除数)	Newton 逆数反復	$O(M(n))$
2 冪の除数	右シフト	$O(n)$

ここで $n$ は被除数のワード数、$m$ は除数のワード数、$M(n)$ は $n$ ワード同士の乗算コストである。

乗算と異なり、除算はオペランドの非対称性が本質的である。被除数と除数のサイズ比がアルゴリズム選択に大きく影響する。

閾値はターゲットハードウェア (x86-64, AVX2) でのベンチマークに基づく。異なるアーキテクチャでは最適な閾値が変わる可能性がある。

Schoolbook 除算 (Knuth Algorithm D)

Knuth の The Art of Computer Programming Vol. 2, §4.3.1 に記載された古典的な長除法を、 64 ビットワード単位で行うアルゴリズムである。

正規化 (Normalization)

除算を開始する前に、除数の最上位ワード (MSW) の最上位ビットが 1 になるように被除数と除数の両方を左シフトする。この正規化により、試行商の推定精度が向上し、補正ステップの回数が高々 2 回に抑えられる。

$$d' = d \cdot 2^s, \quad a' = a \cdot 2^s \quad (s = \text{clz}(d_{\text{MSW}}))$$

ここで $\text{clz}$ は先頭ゼロビット数 (count leading zeros) である。

試行商の推定

各反復で、被除数の上位 2 ワードを除数の最上位ワードで割ることで試行商 $\hat{q}$ を得る:

$$\hat{q} = \left\lfloor \frac{a'_{j+n} \cdot B + a'_{j+n-1}}{d'_{n-1}} \right\rfloor \quad (B = 2^{64})$$

正規化の効果により、$\hat{q}$ は真の商 $q_j$ に対して $q_j \leq \hat{q} \leq q_j + 2$ を満たす。したがって高々 2 回の補正減算で正しい商が得られる。

計算量: $O(nm)$ ($n$ = 被除数のワード数、$m$ = 除数のワード数)

除数が 64 ワード未満の場合、分割統治のオーバーヘッドが Schoolbook のコストを上回るため、 Schoolbook が最速となる。

関連記事: 多倍長整数の基本演算

Burnikel-Ziegler 除算

1998 年に Burnikel と Ziegler が発表した、高速乗算を活用する除算アルゴリズムである。除算を乗算に帰着することで、大きなオペランドに対して Schoolbook 除算より高速に動作する。

基本的なアイデア

$2n$ ワードの被除数を $n$ ワードの除数で割る問題を、 $3n/2$ ワード規模の部分問題に分割して再帰的に解く。各再帰ステップでは Schoolbook 除算ではなく、乗算 (Karatsuba / Toom-Cook / NTT) を使って商と剰余の検証を行う。

具体的には、被除数 $a$ を上位ブロック $a_1$ と下位ブロック $a_0$ に分割し、まず $a_1$ を除数 $d$ で割って部分商 $q_1$ と剰余 $r_1$ を得る。次に $r_1$ と $a_0$ を連結したものを $d$ で割って $q_0$ と $r_0$ を得る。最終的な商は $q = q_1 \cdot B^k + q_0$、剰余は $r = r_0$ である。

計算量: $O(M(n) \log n)$ ($M(n)$ は $n$ ワード同士の乗算コスト)

除数が 64 ワード以上の場合に使用される。高速乗算アルゴリズム (Karatsuba / Toom-Cook / NTT) を活用できるため、大きな除数に対して Schoolbook 除算を大幅に上回る。

関連記事: 多倍長整数の高度な演算

Newton 逆数反復

被除数が除数に対して極端に大きい不均衡なケースに特化したアルゴリズムである。除数 $d$ の近似逆数 $1/d$ を Newton 法で求め、被除数との乗算で商を得る。

Newton 反復

$f(x) = 1/x - d$ の零点を求める Newton 法は以下の反復式になる:

$$x_{k+1} = 2x_k - d \cdot x_k^2$$

この反復は二次収束する。すなわち、各反復で有効精度が 2 倍になる。初期値は除数の上位数ワードから低精度の逆数として求め、目標精度に達するまで反復を繰り返す。

商の計算

十分な精度の逆数 $x \approx 1/d$ が得られたら、商を次のように計算する:

$$q \approx a \times x$$

その後、$r = a - q \cdot d$ で剰余を求め、$0 \leq r < d$ となるよう必要に応じて商を補正する。

計算量: $O(M(n))$ — 乗算 1 回と同等のコスト

実装は GMP の mpn_ni_invertappr (mu_div_qr) に準拠している。逆数の計算自体にコストがかかるため、被除数と除数のサイズ比が十分に大きい場合にのみ使用される (動的閾値)。

関連記事: 多倍長整数の高度な演算

除算バリエーション

calx は標準の切り捨て除算 (operator/) に加えて、 IntOps を通じて以下の除算バリエーションを提供している。 API の詳細については Int API リファレンスを参照。

floorDiv — 床除算

$$\text{floorDiv}(a, b) = \lfloor a / b \rfloor$$

Python の // 演算子と同等である。 C++ の operator/ はゼロ方向への切り捨て (truncation toward zero) だが、 floorDiv は常に負の無限大方向へ丸める。負のオペランドを含む場合に結果が異なる:

$$(-7) / 2 = -3 \quad (\text{C++ 切り捨て}), \quad \text{floorDiv}(-7, 2) = -4 \quad (\text{床除算})$$

ceilDiv — 天井除算

$$\text{ceilDiv}(a, b) = \lceil a / b \rceil$$

正の無限大方向へ丸める除算。メモリ確保サイズの計算など、「少なくとも何個必要か」を求める場面で有用である。

$$\text{ceilDiv}(10, 3) = 4, \quad \text{ceilDiv}(9, 3) = 3$$

divExact — 厳密除算

被除数が除数で割り切れることが事前に分かっている場合に使用する。剰余の検証をスキップするため、通常の除算より高速に動作する。

割り切れない入力を渡した場合の結果は未定義である。

設計判断

Schoolbook / Burnikel-Ziegler 閾値が 64 ワードの理由

Burnikel-Ziegler は乗算を内部的に使用するため、乗算のオーバーヘッドが存在する。除数が小さい場合は Schoolbook の $O(nm)$ のほうが定数倍が軽い。 x86-64 環境でのベンチマーク結果として、64 ワード (約 1,200 桁) がクロスオーバーポイントであった。

Newton 逆数反復を不均衡ケースのみに限定する理由

Newton 逆数反復は $O(M(n))$ と漸近的に最速だが、逆数の計算にはそれ自体が乗算数回分のコストを要する。被除数と除数が同程度のサイズの場合、逆数計算のオーバーヘッドが除算本体のコストを上回る。逆数の恩恵が活きるのは、

被除数が除数に対して著しく大きい (逆数計算のコストが除算全体に対して相対的に小さくなる)
同じ除数で繰り返し除算する (逆数を再利用できる)

といったケースである。calx では被除数と除数のサイズ比を動的に判定し、 Newton が有利な場合にのみ切り替える。

参考文献

Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997. §4.3.1.
Burnikel, C. and Ziegler, J. (1998). "Fast Recursive Division". Max-Planck-Institut für Informatik Research Report, MPI-I-98-1-022.
GMP Development Team. GNU Multiple Precision Arithmetic Library. https://gmplib.org/