第5章: 高速除算アルゴリズム

5.1 学校式除算 (Knuth Algorithm D)

$n$ ワードの被除数 $A = (a_{n-1}, \ldots, a_0)$ を $m$ ワードの除数 $B = (b_{m-1}, \ldots, b_0)$ で割り、 商 $Q$ と余り $R$ ($A = QB + R$, $0 \le R < B$) を求める。 Knuth の The Art of Computer Programming Vol. 2 §4.3.1 に記述される Algorithm D は、すべての高速除算アルゴリズムの土台である。

計算量は $O(nm)$ 回のワード乗算である。$n = m$ のバランス除算では $O(n^2)$、 これは学校で習う筆算と本質的に同じ計算手順である。

sangi では MpnOps.hpp の div_basecase / sbpi1_div_qr がこの役割を担う。基本構造は次の通りである。

// MpnOps.hpp の div_basecase (概略)
inline void div_basecase(uint64_t* q, uint64_t* a, size_t an,
                         const uint64_t* b, size_t bn) {
    const uint64_t d1 = b[bn - 1];   // MSB = 1 (正規化済み)
    const uint64_t d0 = b[bn - 2];
    // 事前計算: dinv = ⌊(β² − 1) / d1⌋ − β  (Möller-Granlund)
    // 反復: 各位置 j で 3-by-2 商推測 + addmul_1 による減算 + 補正
}

5.2 商の推測 — 2-by-1 と 3-by-2

Algorithm D の内ループは「$A$ の上位 3 ワードと $B$ の上位 1〜2 ワードから $\hat{q}$ を推測する」ことに尽きる。 計算時間のほとんどはこのワード単位の除算に費やされるため、この微小な操作の最適化が全体性能を左右する。

2-by-1 の推測

最も単純な推測は上位 2 ワードからの商である:

これは 1 回のハードウェア除算 (DIV 命令) で計算できる。x86-64 の DIV r64 は 128-by-64 bit 除算を約 $20\text{-}30$ サイクルで実行する。

3-by-2 の推測 (Möller-Granlund)

より精度の高い推測として、上位 3 ワードを用いる方法がある。 Möller & Granlund (2011) は、事前計算した逆数 $\text{dinv}$ を使って、 ハードウェア除算を完全に除去した 3-by-2 推測を提示した:

この MULX + 加算で $\hat{q}$ の推測が $2\text{-}3$ サイクルで済む。 事前計算した $\text{dinv}$ は除数が変わらない限り再利用できるため、長い除算では圧倒的に有利である。

// sangi: 3-by-2 商推測 (MpnOps.hpp より)
inline void udiv_qrnnd_3by2(uint64_t& q, uint64_t& r1, uint64_t& r0,
                            uint64_t u2, uint64_t u1, uint64_t u0,
                            uint64_t d1, uint64_t d0, uint64_t dinv) {
    uint64_t p_hi;
    uint64_t p_lo = _umul128(u2, dinv, &p_hi);
    uint64_t q0 = p_lo + u1;
    uint64_t q1 = p_hi + u2 + (q0 < p_lo);
    // 以下: 商の条件付き補正 (多くとも 2 回)
}

5.3 Burnikel-Ziegler 再帰除算

Knuth D は $O(n^2)$ のワード乗算を要する。 サイズが数百 ワードに達すると、乗算を内部で高速化しても全体が $O(n^2)$ に縛られる。 Burnikel & Ziegler (1998) は、分割統治により除算を $O(M(n) \log n)$ に落とした。 ここで $M(n)$ は $n$ ワード乗算のコスト (Karatsuba 以降の高速乗算が使える)。

Div-2n-by-n パターン

被除数 $A$ が $2n$ ワード、除数 $B$ が $n$ ワードで、$A / B$ の商が $n$ ワードに収まる場合を考える。 $B$ と $A$ をそれぞれ $n/2$ ワード単位で分割する:

このとき、商を上下 2 つに分けて順に再帰的に求めることができる:

1. 上半分の商 $Q_1$: $A$ の上位 $3n/2$ ワードを $B$ で割る (再帰 Div-$3n$/2-by-$n/2$)。
2. 下半分の商 $Q_0$: 更新された中間余りと $A$ の下半分から再帰で求める。
3. 結合: $Q = Q_1 \beta^{n/2} + Q_0$。

再帰の各ステップで、サイズが半分の除算と、その結果を使った 1 回の $n/2 \times n$ 乗算が発生する。乗算コストを $M(n)$ とすると:

Master 定理により $D(n) = O(M(n) \log n)$ となる。 $M(n) = O(n^{1.465})$ (Toom-3) や $O(n \log n)$ (FFT) の場合、除算は乗算より $\log n$ 倍遅いだけである。

sangi の実装

// MpnOps.hpp: dcpi1_div_qr_n (Divide-and-Conquer, pre-inverted)
inline uint64_t dcpi1_div_qr_n(uint64_t* qp, uint64_t* np,
                               const uint64_t* dp, size_t n,
                               uint64_t dinv, uint64_t* scratch) {
    if (n < DC_DIV_QR_THRESHOLD) {
        return sbpi1_div_qr(qp, np, 2 * n, dp, n, dinv); // 基底: 学校式
    }
    size_t lo = n / 2, hi = n - lo;
    // 上半分: 2hi-by-hi 再帰
    uint64_t qh = dcpi1_div_qr_n(qp + lo, np + 2 * lo, dp + lo, hi, dinv, scratch);
    // 補正: Q1 × B_low の差し引き、下半分の再帰
    // ...
}

sangi では 2 つの閾値を区別している (値は 2026-05-08 時点):

• BZ_THRESHOLD = 24: 最上位ディスパッチの切替点。除数が 24 ワード未満なら Knuth D 直接、24 ワード以上なら BZ 再帰に入る。
• DC_DIV_QR_THRESHOLD = 20: BZ 再帰内部の基底ケース閾値。再帰でサイズが 20 ワード未満まで縮んだら sbpi1_div_qr (学校式) に落とす。

入口閾値 24 と基底閾値 20 の差は小さいが、両者を独立に持つことで「サイズが 24-31 ワードあたりで BZ に入った直後、最初の再帰分割でちょうど 20 ワード台に落ちて基底ケースに切替わる」自然な遷移になる。なお BZ_THRESHOLD は当初 64 だったが、2026-05-08 に sbpi1 (学校式) ループの高速化を反映して 24 に下げられた。

5.4 Newton 逆数反復

除算 $A / B$ を 2 段階に分ける:

1. $B$ の逆数近似 $X \approx 1/B$ を高速に求める。
2. $A \cdot X$ で商を計算し、誤差を補正する。

このアプローチでは、除算を 乗算に還元できる。 乗算 $M(n)$ が FFT 等で $O(n \log n)$ なら、除算も同じオーダーになる。 Newton 反復はこの戦略の要である。

Newton 反復式

$f(x) = 1/x - B$ の零点 $x = 1/B$ に対して Newton 法を適用する:

$x_k$ が $1/B$ に近いとき $B x_k \approx 1$、$2 - B x_k \approx 1$ で、 $x_{k+1}$ の誤差は $x_k$ の誤差の 2 乗に比例する (二次収束)。 精度 $b$ ビットが必要なら、初期近似 $x_0$ (たとえば倍精度浮動小数の 48 ビット) から始めて $\lceil \log_2(b / 48) \rceil$ 回の反復で十分。

sangi の Float 除算での実装

sangi では、Float の除数が 64 ワード (約 4096 ビット) 以上で Newton に切り替わる。 初期近似は倍精度による逆数、反復中は精度を倍々に拡張する。

// Float.cpp: Newton 除算の初期近似
int64_t b_bl = static_cast<int64_t>(b_work.mantissa_.bitLength());
int64_t exp_adj = b_work.exponent_ + b_bl - 1;  // b_norm ∈ [1, 2)
b_work.exponent_ = -(b_bl - 1);
double b_d = b_work.toDouble();                  // 53 bit 近似
Float x(1.0 / b_d);
x.exponent_ -= exp_adj;                          // 元のスケールに戻す
x.effective_bits_ = 64;                          // 達成精度の記録

// Newton 反復本体 (b_work は正規化済みの除数)
int correct_bits = 48;
while (correct_bits < target_bits) {
    Float bx = b_work * x;                       // ≈ 1
    Float residual = one - bx;                   // ≈ 2^{-correct_bits}
    x = x + x * residual;                        // x_{k+1} = x_k (2 − B x_k)
    correct_bits = std::min(2 * correct_bits, target_bits);
}

精度追跡 (effective_bits)

Newton 反復の精度管理は微妙である。$k$ 回目の 目標 精度 $p_k$ と、 実際に達成された 効果的 精度 $\text{effective\_bits}$ は必ずしも一致しない。 反復で用いる丸めの方向、入力誤差の伝播、加減算での桁落ちが原因である。

sangi の実装では、各 Float が effective_bits_ フィールドを持ち、 反復毎に min(2 * correct_bits, step_bits) で更新する。 最終商の計算でも誤差を踏まえた 1 回の補正を行う:

これは「最終の $M(P, P)$ 乗算を一度節約しつつ誤差を 2 乗で縮める」 fused 補正である。

5.5 Exact Division

「余りが $0$ になることがあらかじめ分かっている」除算を特別扱いできる。 たとえば多項式除算での $A = QB$ の解消、二項係数 $\binom{n}{k} = n! / (k!(n-k)!)$ の整数性、 Toom-Cook 補間で出てくる固定除数 (/3, /5, /11 など) — いずれも余りが $0$ の除算である。

Jebelean の Exact Division

余りを切り捨てない普通の除算では商を上位から求める。 しかし余りが $0$ なら、商を 下位から求めても 結果は同じである。 Jebelean (1993) のアイデアは、$B$ の逆元 $B^{-1} \bmod \beta$ を $p$-進的に使うことである:

これを $n$ 回繰り返せば商が得られる。各ステップは 1 ワード乗算 + addmul_1 でハードウェア ネイティブ。Möller-Granlund の逆数推測も不要。結果として:

用途

• Toom-Cook 補間: 評価値を Vandermonde 逆行列で戻すとき、固定の小整数で割る (第8章 参照)。 Toom-8 では /11, /17, /31 といった小素数除算が 10〜20 回入る。
• 二項係数: $\binom{n}{k}$ を逐次更新 $\binom{n}{k} = \binom{n}{k-1} \cdot (n-k+1) / k$ で 求めると、各ステップは exact division。
• Binary Splitting: 有理級数の結合で共通因子の除去 (第10章 参照)。

5.6 特殊定数の除算

2 の冪での除算

$A / 2^k$ は単なる右シフトである。 ワード境界に揃っていれば ($k$ が 64 の倍数) ポインタの付け替えだけで済む。 揃っていなくても各ワードを SHRD 命令で 2 ワードまたぎにシフトすればよい。 コストは $O(n)$ の線形時間だが、乗除算と違って分岐も補正もない。

10 の冪での除算

10 進出力や有効桁管理では $A / 10^k$ が頻繁に必要になる。 $10 = 2 \cdot 5$ なので、$A / 10^k = (A / 2^k) / 5^k$。右シフト後に 固定除数 $5^k$ で exact division する。

小定数での除算

$d \in \{3, 5, 7, \ldots\}$ のような小定数での除算は、 Möller-Granlund の divide by invariant integer 法で実装する:

$\mu$ は事前計算できるので、実行時は 1 回の MULH (上位半分の乗算) とシフト。 最大で 1 回の補正が必要だが、それを含めても 1 ワード除算の数倍速い。

5.7 アルゴリズム選択の閾値 — sangi の設計

高速アルゴリズムには、基本的に「定数項が大きいが漸近的に優位」という性質がある。 小さい入力では素朴なアルゴリズムの方が速いため、サイズごとにディスパッチする必要がある。 sangi の除算ディスパッチは次の通り。

なぜ Int では $\mu$-div を使い Float では Newton を使うのか — Int の商は整数のため 余り $R$ を同時に計算しなければならず、Barrett 還元ベースの $\mu$-div が向いている。 一方 Float は商だけを必要とし、精度も離散的でなく連続的に管理できるため Newton が自然である。

5.8 実測性能

sangi の除算を代表的なサイズで実測し、サイズ帯別の特性を観察する。 測定環境は Zen 3 (Ryzen 5 5600X), MSVC Release x64, single thread。

サイズ帯別の実装特性

除算の速度は 除算そのものではなく、内部で呼ばれる 乗算の速度で大部分が決まる。 Newton ベースの $\mu$-div では逆数計算に不均衡乗算 (被乗数 $\gg$ 乗数) が多数発生する。

• NTT 域 ($n \ge 1200$ ワード): 5-smooth NTT が $O(n \log n)$ の漸近優位を発揮する領域。
• Karatsuba 域 ($22 \le n < 140$ ワード): 古典的 Karatsuba による $O(n^{1.585})$ が最速。
• Toom-6/Toom-8 域 (600-1200 ワード): 高次 Toom-Cook が最速。Toom-8 の pre-pad 評価が効く。
• invert_approx の 2:1 不均衡乗算 (513-1025 ワード): Toom-4,2 特化で不均衡のピースサイズを揃える。

つまり「除算を速くする」には「除算の内側の乗算を速くする」のが本筋である。 524K bit の除算が大きく短縮された要因も、 5-smooth NTT (第9章) の導入と Toom-8 の最適化 (2026-04-19) にあった。

5.9 まとめ

• Knuth D と 3-by-2 推測: 小サイズでは学校式 + Möller-Granlund の事前計算逆数が最速。 推測誤差が $\{0, 1\}$ に収まるよう正規化が前提。
• Burnikel-Ziegler: 分割統治で $O(M(n) \log n)$ を達成。 sangi は 24 ワード (BZ_THRESHOLD) から BZ に切り替える。
• Newton 逆数反復: 除算を乗算に還元、$O(M(n))$ 達成。 二次収束なので倍精度初期値から数回で最終精度に到達。 精度追跡 (effective_bits) を誤ると商が壊れる。
• Exact Division: 余り 0 が保証される場合は Jebelean の $p$-進版で $O(n)$。 Toom-Cook 補間の小定数除算に多用される。
• 性能の主因は乗算: 除算の速度は大部分「内部で呼ばれる乗算」の速度で決まる。 sangi が 524K bit で達成した改善も、NTT と Toom-Cook の最適化に由来する。

次章では、除算の特殊形であるモジュラー算術 ($a \bmod m$) を扱う。 RSA や NTT のホットループで現れる Barrett 還元・Montgomery 乗算は、 本章の除算を「除算しない」方向に再構成する技術である。

参考文献

• Knuth, D.E. The Art of Computer Programming, Vol. 2: Seminumerical Algorithms, 3rd ed., Addison-Wesley, 1997. §4.3.1 (Algorithm D)。
• Burnikel, C., Ziegler, J. "Fast Recursive Division", MPI-I-98-1-022, Max-Planck-Institut für Informatik, 1998.
• Möller, N., Granlund, T. "Improved Division by Invariant Integers", IEEE Trans. Computers, 60(2), pp. 165-175, 2011.
• Brent, R.P., Zimmermann, P. Modern Computer Arithmetic, Cambridge University Press, 2010. 第1章。
• Jebelean, T. "An Algorithm for Exact Division", J. Symbolic Computation, 15(2), pp. 169-180, 1993.
• Hasselström, K. "Fast Division of Large Integers — A Comparison of Algorithms", Master's thesis, KTH, 2003.

関連章

• 第4章 高速乗算アルゴリズム — 除算が還元する対象としての乗算
• 第6章 モジュラー算術 — 「除算しない」モジュラー還元の技法
• 第7章 整数 GCD アルゴリズム — 除算を繰り返す Euclid 法とその改良
• 第8章 不均衡 Toom-Cook 乗算 — 不均衡乗算が除算性能に与える影響
• 第9章 NTT 実装比較 — 大サイズ乗算を支える NTT
• [多倍長浮動小数] 第10章 Binary Splitting と超幾何級数 — 最終結合での 1 回の大整数除算