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

6.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·β + d0)⌋ − β  (Möller-Granlund の 3/2 逆数、invert_pi1)
    // 反復: 各位置 j で 3-by-2 商推測 + addmul_1 による減算 + 補正
}

6.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 回)
}

6.3 Burnikel-Ziegler 再帰除算

Knuth D は $O(n^2)$ のワード乗算を要する。 サイズが数百 ワードに達すると、乗算を内部で高速化しても全体が $O(n^2)$ に縛られる。 Burnikel & Ziegler (1998) は、分割統治により除算のコストを乗算で表せるようにした（乗算の速さに応じて $O(M(n))$ または $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$ ワード $(A_3, A_2, A_1)$ を $B$ で割る（$3n/2$-by-$n$）。中身は、上位 $n$ ワード $(A_3, A_2)$ を $B_1$ で割る半サイズの再帰（$n$-by-$n/2$）で仮の商 $\hat{Q}_1$ を得て、$A_1$ を下ろして $\hat{Q}_1 B_0$ を引き、負なら $B$ を足し戻す補正（最大 2 回）である。
2. 下半分の商 $Q_0$: 中間余り $R_1$ と $A_0$ を並べた $3n/2$ ワードを、同じ手順で $B$ で割る。
3. 結合: $Q = Q_1 \beta^{n/2} + Q_0$。

再帰の各段で、サイズが半分の除算（$n$-by-$n/2$）2 回と、$\hat{Q} B_0$ を作る $n/2 \times n/2$ 乗算 2 回が発生する。乗算コストを $M(n)$ とすると:

この漸化式の解は $M(n)$ の増大の速さで変わる。$M(n) = \Theta(n^{\alpha})$（$\alpha > 1$。Toom-3 なら $\alpha \approx 1.465$）のときは、各段のコストが等比級数的に減るので Master 定理の第 3 の場合にあたり $D(n) = O(M(n))$、すなわち除算は乗算と同じ漸近計算量になる。 一方 $M(n) = \Theta(n \log n)$（FFT）のように準線形のときは第 2 の場合にあたり $D(n) = O(M(n) \log n) = O(n \log^2 n)$ で、乗算より $\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 に入ると、最初の再帰分割で 12〜16 ワードになり、基底閾値 20 を下回って学校式に切替わる」自然な遷移になる。なお BZ_THRESHOLD は当初 64 だったが、2026-05-08 に sbpi1 (学校式) ループの高速化を反映して 24 に下げられた。

6.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$ (たとえば倍精度浮動小数。53 ビットの仮数を、丸めの余裕を見て 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;                          // 精度フィールドの初期化（反復は correct_bits = 48 から数える）

// 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, step_bits);  // step_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 補正である。

6.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$-進的に使うことである（逆元が存在するのは $B$ が奇数、すなわち $\gcd(B, \beta) = 1$ のとき。偶数除数の扱いは下の「divexact_by_odd の威力」を参照）:

これを商のワード数 $n$ 回繰り返せば商が得られる。各ステップは $q_0$ を決める 1 ワード乗算と、 $B$（$m$ ワード）の addmul_1 で、商の推測も補正も要らない。Möller-Granlund の逆数推測も不要。結果として:

用途

• Toom-Cook 補間: 評価値を Vandermonde 逆行列で戻すとき、固定の小整数で割る (第9章 参照)。 Toom-8 では /3, /5, /7, /11, /13 といった小素数（とその積）での除算が補間に 60 回余り入り、sangi はそれを 8 回の divexact_by_odd に折り畳む。
• 二項係数: $\binom{n}{k}$ を逐次更新 $\binom{n}{k} = \binom{n}{k-1} \cdot (n-k+1) / k$ で 求めると、各ステップは exact division。
• Binary Splitting: 有理級数の結合で共通因子の除去 (第10章 参照)。

6.6 特殊定数の除算

2 の冪での除算

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

10 の冪での除算

10 進出力や有効桁管理では $A / 10^k$ が頻繁に必要になる。 $10^k \mid A$ が分かっているとき（末尾の 10 進ゼロを落とす等）は、$10 = 2 \cdot 5$ なので $A / 10^k = (A / 2^k) / 5^k$。右シフト後に 固定除数 $5^k$ で exact division すればよい。 一般の $\lfloor A / 10^k \rfloor$（基数変換や桁の切り捨て）は割り切れないので exact division は使えない（黙って誤った商を返す）。 商と余りを通常の除算で求めるか、$5^k$ が 1 ワードに収まるなら下の小定数での除算を使う。

小定数での除算

$d \in \{3, 5, 7, \ldots\}$ のような小定数での除算は、 Granlund–Montgomery (1994) の division by invariant integers 法で実装する（$a$ は 64 ビット、$\ell = \lceil \log_2 d \rceil$）:

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

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

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

なぜ Int では $\mu$-div を使い Float では Newton を使うのか — Int の商は整数のため 余り $R$ を同時に計算しなければならず、Barrett 還元ベースの $\mu$-div が向いている。 一方 Float は商だけを必要とし、厳密な余りを返す必要がなく必要な精度まで近似できればよいため、Newton が自然である。

6.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 (第10章) の導入と Toom-8 の最適化 (2026-04-19) にあった。

6.9 まとめ

• Knuth D と 3-by-2 推測: 小サイズでは学校式 + Möller-Granlund の事前計算逆数が最速。 正規化が前提で、3-by-2 の候補は多くとも 2 回の補正で確定する。
• Burnikel-Ziegler: 分割統治。$M(n) = \Theta(n^{\alpha})$（$\alpha>1$）なら $O(M(n))$、FFT のような準線形乗算なら $O(M(n) \log n)$。 sangi は 24 ワード (BZ_THRESHOLD) から BZ に切り替える。
• Newton 逆数反復: 除算を乗算に還元、$O(M(n))$ 達成。 二次収束なので正しい桁数が毎回ほぼ倍増し、倍精度の初期値から $n$ ビットまで $\lceil \log_2(n/48) \rceil$ 回程度で到達する（実用サイズでは数回）。 精度追跡 (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.
• Granlund, T., Montgomery, P.L. "Division by Invariant Integers using Multiplication", Proc. ACM SIGPLAN PLDI '94, pp. 61-72, 1994.
• 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.

関連章

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