定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム

競プロキャンプ2026関東LT, 2026-06-06

Transcript

1. 定数除数整数除算の高速実装 近似商と補正による除算アルゴリズム Mizar/みざー 2026-06-06
2. 0. 自己紹介代わりの問題紹介 Accuracy of Integer Division Approximate Functions (2023年8月出題) [1] Max Weighted Floor of Linear (2025年12月出題) [2] Accuracy of Integer Division Approximate Functions 2 (2025年12月出題) [3] 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [1] Mizar, “No.2440 Accuracy of Integer Division Approximate Functions.” yukicoder, 2023, [Online]. Available: https://yukicoder.me/problems/no/2440. [URL] [2] Mizar, “No.3397 Max Weighted Floor of Linear.” yukicoder, 2025, [Online]. Available: https://yukicoder.me/problems/no/3397. [URL] [3] Mizar, “No.3398 Accuracy of Integer Division Approximate Function 2.” yukicoder, 2025, [Online]. Available: https://yukicoder.me/problems/no/3398. [URL] Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 2 / 56
3. Accuracy of Integer Division Approximate Functions (2023年8月出題) 個のクエリが与えられる。 番目のクエリ では、整数 が与えられる。 かつ となるような 整数 がいくつあるかを 数え上げ、 とおく。 この整数 を答えよ。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 3 / 56
4. Accuracy of Integer Division Approximate Functions の解法の要点 を求める。 ならば、答えは である。 ならば、 は の範囲 となりうる 上限 を用いて以下のように表せる。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 4 / 56
5. Max Weighted Floor of Linear (2025年12月出題) 個のテストケースがあります。 各テストケースについて、 正の整数 と 整数 が与えられます。 を求めてください。 つまり、整数 を 以上 未満で動かしたとき、式 が 取りうる値の中で最大のものを求めてください。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 5 / 56
6. Max Weighted Floor of Linear の解法の要点 求める値は である。 は非減少な階段関数であり、 となる各段では 目的関数は という一次関数になる。 そのため各段では端点だけを見ればよく、 なら右端、 なら左端が 候補になる。段の端点を floor 式で表すと、候補全体の最大化は同じ形の小問題に 帰着できる。正規化後の再帰では分母が と Euclid 型に減る ため、各テストケースを で処理できる。 別解として、 floor_sum の再帰をモノイド積に一般化した floor_prod を使う方 法もある。 の増加と floor 値の増加をそれぞれ重み のステップと見なし、 累積和と最大 prefix を持つモノイド上で計算すると、目的関数の最大値を同じ Euclid 型の反復で求められる。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 6 / 56
7. Accuracy of Integer Division Approximate Functions 2 (2025年12月出題) 各テストケースについて、 正の整数 と 非負整数 が与えられます。 次に示す関数 を定義します。 を満たす非負整数 が存在する場合は、 を満たす最小の非負整数 を答えてください。存在しない場合は、 を答えてください。 出力すべき値は、 64bit 整数で表記できる範囲を超えることがあります。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 7 / 56
8. Accuracy of Integer Division Approximate Functions 2 の解法の要点 とおく。 が の倍数に達したときだけ の値は増加しうるため、閾値を 初めて超える候補を調べるには、まず の形に注目する。 とした場合、 誤差の閾値条件は と同値変形できる。この不等式 を用いて、区間内に を 満たす候補が存在するかを Max Weighted Floor of Linear で判定し、その判定を用いて最小の を二分探索/指数探索する。また不等式 からは、 が を満たすときは を満たす候補が存在することもわかるため、 二分探索/指数探索の上限を に絞ることもできる。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 8 / 56
9. 1. 概要 整数除算では、被除数 と 除数 に対して、商 や剰余 を 求める。除算命令は、乗算やシフトに比べて高コストである。除数 が 定数なら、 の固定小数点近似 を事前計算し、除算を の 上位ビット取得と少数回の補正に置き換えられる。このスライドでは、 定数除数整数除算を 逆数近似の掛け算 + 軽い後処理 として実装する方法と、 その正確性条件を扱う。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 9 / 56
10. 除算によらない方法への変換 コンパイラは時に、整数除算を乗算すら用いない方法に変換する事もあります。 これは数学的には と変形しています。 https://rust.godbolt.org/z/oETab1ssr 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 10 / 56
11. 2. ハードウェア除算の手法とコスト 除算命令と乗算命令のコスト比較 なぜハードウェア除算は遅いのか？ SRT 除算器 Goldschmidt 除算 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 11 / 56
12. 除算命令と乗算命令のコスト比較 https://uops.info/ による Arrow Lake-P の 実測値を見ると、整数除算命令 DIV は、乗算命令 IMUL / MUL / MULX と比べてレイテンシ (Latency)・スループット (Throughput, TP) の 両面で重い。 DIV r64 はレイテンシが 16–19 cycles、 スループットが約 10.55 cycles である。一方、64- bit 乗算系命令は、 IMUL r64,r64 なら 3–4 cycles / 0.36 cycles、 MUL r64 や MULX r64,r64,r64 でも 3–5 cycles / 1.00 cycle 程度である。 したがって、除数が定数である場合に、除算を 乗算 + シフト + 少数回の補正 へ置き換える ことには、命令レベルでも明確な動機がある。 命令 内容 Latency TP μops DIV r32 EDX:EAX / r32 11–16 6.05 4 DIV r64 RDX:RAX / r64 16–19 10.55 4 IMUL r32,r32 op1 := op1 * op2 3 0.36 1 IMUL r64,r64 op1 := op1 * op2 3–4 0.36 1 MUL r32 EDX:EAX := EAX * r32 3–4 1.00 3 MUL r64 RDX:RAX := RAX * r64 3–5 1.00 2 MULX r32,r32,r32 op1:op2 := EDX * op3 (BMI2) 4 0.94 3 MULX r64,r64,r64 op1:op2 := RDX * op3 (BMI2) 3–5 1.00 2 SHRD op1 := 4 1.00 1 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 12 / 56
13. なぜハードウェア除算は遅いのか？ 64bit整数の高性能な除算器としては、SRT除算器を用いた実装が一般的である。SRT除算器は、いわば「除算の商を 2,4,6bitずつ段階的に求めていく」回路である。SRT除算器は、加算器やシフト回路を複数段重ねる構造をしているため、 遅延時間が ～ となる。 -bit 整数演算の遅延時間と回路規模の例を以下に示す。 演算 代表的実装 遅延時間 回路規模 加算・減算 桁上げ伝播加算器 Ripple Carry Adder (RCA) 加算・減算 桁上げ先読み加算器・接頭辞加算器 Carry Look-ahead Adder (CLA) Prefix Adder シフト バレルシフタ Barrel Shifter 乗算 Wallace/Dadda Tree + 高速加算器 除算 SRT除算器 + 高速加算器など 除算 Goldschmidt / Newton-Raphson 回の乗算 乗算器依存 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [1] [2] [1] クリティカルパスの遅延時間のオーダー。実際のサイクル数などは、回路規模の制限・実装技術などにより異なる。 [2] AMD-K7(Athron, 1999年) など。四倍精度浮動小数点数以上の精度向けに候補として挙げられることがある。 Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 13 / 56
14. SRT 除算器 SRTは考案者 Sweeney, Robertson, Tocher の頭文字であり、 2進の筆算型除算を冗長表現付きに拡張した ハードウェア実装の手法である。 Radix-R SRT では、各段で 商の1桁を冗長な桁集合から後続桁で補正可能な範囲で 近似的に選ぶ。Radix-4で最大冗長表現を用いた場合、 {-3,...,+3} の候補から商桁を求める。このため1段で商を bit 相当分進められる。また、少しの拡張でSRT平方 根器も作れる。 近年のCPUの整数除算器は、 SRT Radix-4 除算器を 2~3段重ねて実行し、1サイクルで商の桁を 4bit, 6bit 相 当進める設計もあり、 Radix-16, Radix-64 と称される こともある。(Radix-16 は 2006年の Intel Penryn など) 1993年、Intel Pentium プロセッサで SRT 除算器の実装ミスに より、浮動小数点数の除算が特定の被除数・除数の組に対して 誤った値を返すバグが発生した。これにより、Intelは大規模な リコールを行うことになった。 右図は Radix-4 の桁選択表の例で、被除数（縦軸）と除数 （横軸）の上位bitから次の商桁（各グリッドセル）を求める。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 14 / 56
15. 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 15 / 56
16. 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 16 / 56
17. Goldschmidt 除算 の分子と分母に共通の補正係数を掛けることを繰り返し、分母を に近づけていくのが Goldschmidt 除算の基本的な考え方である。例として、被除数 と除数 に共通の 冪スケーリングを行い、除数を または の範囲に正規化済みであるとする。 とおくと なので、 は に収束する。したがって、 と展開できる。実装上は、 とおき、各段で を計算する。このとき である。有限段で打ち切ると残る誤差は の部分に対応する。 は各段で二乗されるため急速に に近づき、十分な段数の後では が の近似値になる。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 17 / 56
18. 3. 数学的基礎と定数除数整数除算のアルゴリズム 床関数と天井関数 定数除数整数除算の高速化 商の近似関数・近似誤差の定義 中国剰余定理・Garner のアルゴリズムと Montgomery reduction Barrett Reduction GM method 上側近似整除法の正確性条件 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 18 / 56
19. 床関数と天井関数 床関数 は、 以下の最大の整数である。すなわち、 である。 不等式で書くと、 である。 天井関数 は、 以上の最小の整数である。すなわち、 である。 不等式で書くと、 である。 基本的な不等式として、 、および が成り立つ。 整数 に対して、 が成り立つ。 床関数と天井関数は、 により相互に変換できる。 正の整数 と 実数 に対して、 が成り立つ。 実際、 と書くと、 な ので、 である。 非負整数 と正の整数 に対して、 と書けば、 である。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 19 / 56
20. 定数除数整数除算の高速化 を毎回ハードウェア除算命令で求めるのは、乗算やシフトに比べて高コ ストである。 今回は、除数 が定数である場合に、 という固定小数点近似を使い、除算を 乗算 + シ フト + 少数回の補正 へ置き換えるアルゴリズムを主に扱う。 は整数型のビット幅、 、 を想定する。 場合によっては、対象範囲すべてで近似商が真の商と一致するように と を選ぶ。 典型的には、正しい商 を得てから剰余 を求める。 近似商 から を作り、 と を補正して正しい商 と剰余 を得る方法もある。 剰余 のみが欲しい場合、商を明示的に求めずに剰余を直接計算する方法もある。 被除数・除数の範囲や性質、欲しい値が商だけ、剰余だけ、商と剰余の両方なのかに応じて、様々なア ルゴリズムが存在する。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 20 / 56
21. 商の近似関数・近似誤差の定義 は非負整数、 は正の整数 として、商 を計算するための近似関数・近似誤差関数を 定義する。小さな整数型では主に とした議論を行う。 と の両方について、床関数を 取ると下側近似、天井関数を取ると上側近似とここでは呼ぶ。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 21 / 56 の下側近似 商の下側近似 の上側近似 商の上側近似 の時に生じる下側近似の誤差 商の真値 商の下側近似 の時に生じる上側近似の誤差 商の上側近似 商の真値
22. 中国剰余定理・Garner のアルゴリズムと Montgomery reduction 互いに素な正整数 に対し、 とする。このとき、整数 から を作ると、 かつ かつ が成り立つ。これは 2 つの法に対する中国剰余定理の解を、片側から順に復元する形になっている。 ここで とおく。 とすれば、 は かつ を満たす。したがって は で 割り切れ、 とおくと、 に対して となる。 これは Montgomery reduction [4] の REDC の中核である。 Montgomery reduction に関連して、 における乗法逆元を高速に計算する手法 [5] や、偶数 に対する Montgomery 乗算の手法 [6] も存在する。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [4] P. L. Montgomery, “Modular multiplication without trial division,” Mathematics of Computation, vol. 44, no. 170, pp. 519–521, 1985. [5] J. Hurchalla, “An Improved Integer Modular Multiplicative Inverse (modulo 2 ).” 2022, [Online]. Available: https://arxiv.org/abs/2204.04342. [URL] [6] Nachia, Mathenachia, “mod 偶数の乗算にモンゴメリ乗算を使いたい？.” [Online]. Available: https://www.mathenachia.blog/even-mod-montgomery-impl/. [URL] w Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 22 / 56
23. 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 23 / 56
24. Barrett Reduction Paul Barrett の CRYPTO’86 論文 [7] は、RSA を標準的な DSP 上で高速に実装するため、 多倍長除算を避ける modulo reduction を提案した。基数 の 桁の整数 で除算を 行うとき、 を事前計算しておき、 を除算ではなく、乗算・桁の切り捨て・減算で求める。 現在は Barrett reduction と呼ばれる、固定された除数 に対する商・剰余計算を、 事前計算した 逆数近似・乗算・桁の切り捨て・補正減算で行う古典的手法である。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [7] P. Barrett, “Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor,” in Advances in Cryptology — CRYPTO’ 86, 1987, pp. 311–323. Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 24 / 56
25. Barrett Reduction の誤差評価 前提・定義 は整数、 主張 証明 と書くと、 より まず より したがって 一方、 かつ より よって また より したがって よって ここで なので、 したがって 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 25 / 56
26. GM method 除数 の逆数 を として近似し、 がすべての で成り立つように と を選ぶ方法が [8] で提案されました。 論文の Theorem 4.2 では、次の十分条件を与えています。 と を求める方法は、 [9][10] にも記載されています。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [8] T. Granlund and P. L. Montgomery, “Division by invariant integers using multiplication,” in Proceedings of the ACM SIGPLAN 1994 Conference on Programming Language Design and Implementation, 1994, pp. 61–72. [9] Jr. Warren Henry S., Hacker’s Delight, 2nd ed. Addison-Wesley, 2013. [10] 著ヘンリー・S・ウォーレン、ジュニア and 滝沢徹, ハッカーのたのしみ : 本物のプログラマはいかにして問 題を解くか, オンデマンド版. エスアイビー・アクセス, 2014. [URL] Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 26 / 56
27. GM method の正確性条件 Theorem 4.2: unsigned division を乗算で正確に計算する条件 bit 非負整数 を定数 で割ることを考える。 主張 を非負整数、 とし、 が成り立つとする。このとき、任意の整数 について つまり、 を で近似しても、 bit 入力に対する 商は変わらない。 証明の要点 とおくと、仮定より 、ただし と書く。 このとき 右辺は非負であり、さらに より したがって よって 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 27 / 56
28. 定理: 上側近似整除法の正確性条件 非負整数の被除数 と、正の整数の除数 による整除法を考える。正の整数 に対し とおくと、以下の必要十分条件がそれぞれ成り立つ。 かつ この一般化した必要十分条件は、以降の 32bit/64bit 実装の正当性に帰着する。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 28 / 56
29. 定理: 上側近似整除法の正確性条件 非負整数の被除数 と、正の整数の除数 による整除法を考える。正の整数 に対し とおくと、以下の必要十分条件が それぞれ成り立つ。 基本変形 とおく。 より ここで とおくと さらに 商の正確性 また、(2) より よって 商と剰余の同時復元 商が正しい、すなわち のとき、 (1) より したがって、(2) から ここで より また なので、 なら であり、上の商条件より 商も正しい。よって 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 29 / 56
30. 4. 32bit/64bit 定数除数整数除算の実装例 32bit整数の除算器 実装例とベンチマーク 64/32bit整数の除算器 実装例とベンチマーク 64bit整数の除算器 Up , Even , Down , Overflow の分岐の正当性 実装例とベンチマーク ac_library に見る実装の住み分け 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 30 / 56
31. 32bit整数の除算器 #[derive(Debug, Clone, Copy)] enum DivisorU32Kind { Up, Shift } #[derive(Debug, Clone)] pub struct DivisorU32(DivisorU32Kind, u64, u32); impl DivisorU32 { pub const fn new(d: u32) -> Self { assert!(d > 0); if d.is_power_of_two() { Self(DivisorU32Kind::Shift, d.trailing_zeros() as u64, d) } else { Self(DivisorU32Kind::Up, u64::MAX / d as u64 + 1, d) } } /// return (x / d, x % d) pub const fn divrem(&self, x: u32) -> (u32, u32) { match self { &Self(DivisorU32Kind::Up, m, d) => { let t = m as u128 * x as u128; ((t >> 64) as u32, ((t as u64 as u128 * d as u128) >> 64) // [!annotate label="Up" note=" 右の系で N = L = 32, β = 2^ δ ≤ 2^L とした形に相当。 D = 1 の場合は M^+ = 2^64 であるため実装に注意。 } &Self(DivisorU32Kind::Shift, s, d) => (x >> (s as u32), x & } } } 系: による -bit 整除法 を整数、 を非負整数とし、 とおく。 を 仮定すると、任意の について が成り立つ。剰余の式は [11] の Theorem 1 の形に対応する。 証明 より以下が成り立つため、 前命題より結論を得る。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [11] D. Lemire, O. Kaser, and N. Kurz, “Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries.” 2019, [Online]. Available: https://arxiv.org/abs/1902.01961. [URL] Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 31 / 56
32. 32-bit同士の剰余算の実装(1) として、 を計算する 2つの実装例を比較する。 の場合は であるため、実装上の注意が必要である。 pub struct DivisorU32U32A(u64, u32); impl DivisorU32U32A { pub const fn new(d: u32) -> Self { assert!(d > 1); Self((u64::MAX / d as u64).wrapping_add(1), } /// returns (x % d) #[target_feature(enable="bmi2")] pub const fn rem(&self, x: u32) -> u32 { (self.0.wrapping_mul(x as u64) as u128 * self.1 as u128 >> 64) as u32 } } pub struct DivisorU32U32B(u64, u32); impl DivisorU32U32B { pub const fn new(d: u32) -> Self { assert!(d > 1); Self((u64::MAX / d as u64).wrapping_add(1), } /// returns (x % d) #[target_feature(enable="bmi2")] pub const fn rem(&self, x: u32) -> u32 { let q = ((self.0 as u128 * x as u128) >> 64) x.wrapping_sub(self.1.wrapping_mul(q as u32) } } 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 32 / 56
33. 32-bit同士の剰余算の実装(2) example::DivisorU32U32A::rem::hb2832088297c2a3c: mov edx, esi ; edx = x imul rdx, qword ptr [rdi] ; rdx = rdx * M^+ mov eax, dword ptr [rdi + 8] ; eax = D mulx rax, rax, rax ; rax = (rax * rdx ret ; return eax example::DivisorU32U32B::rem::hcd863e667530d983: mov eax, esi ; eax = x mov edx, esi ; edx = x mulx rcx, rcx, qword ptr [rdi] ; rcx = (M^+ * rd imul ecx, dword ptr [rdi + 8] ; ecx = ecx * D sub eax, ecx ; eax = eax - ecx ret ; return eax 除数 に関して前計算済みなら、実行時は転送・乗算・減算命令のみで剰余を 計算できるため、除算命令を使うより高速になることが期待される。ただし、この 手法では 64bit 乗算を必要とするため、 SIMD 化との相性はあまり良くない。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 33 / 56
34. 64/32-bitの剰余算の実装(1) として、 を計算する 2つの実装例を比較する。 pub struct DivisorU64U32A(u128, u32); impl DivisorU64U32A { const MASK96: u128 = u128::MAX >> 32; pub const fn new(d: u32) -> Self { assert!(d > 0); Self((Self::MASK96 / d as u128) + 1, d) } /// returns (x % d) #[target_feature(enable = "bmi2")] pub const fn rem(&self, x: u64) -> u32 { ((self.0.wrapping_mul(x as u128) & Self::MAS self.1 as u128 >> 96) as u32 } } pub struct DivisorU64U32B(u64, u32); impl DivisorU64U32B { pub const fn new(d: u32) -> Self { assert!(d > 1); Self((u64::MAX / d as u64).wrapping_add(1), } /// returns (x % d) #[target_feature(enable = "bmi2")] pub const fn rem(&self, x: u64) -> u32 { let q = ((self.0 as u128 * x as u128) >> 64) let (r, b) = x.overflowing_sub(q * self.1 as r.wrapping_add((b as u64) * self.1 as u64) a } } 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 34 / 56
35. 64/32-bitの剰余算の実装(2) example::DivisorU64U32A::rem::h15a37c417d12d6b4: mov rdx, rsi ; rdx = x mulx rcx, rax, qword ptr [rdi] ; [rcx:rax] = M^+[6 mov esi, dword ptr [rdi + 8] ; esi = M^+[95:64] imul esi, edx ; esi = esi * edx add ecx, esi ; ecx = ecx + esi mov esi, dword ptr [rdi + 16] ; esi = D mov rdx, rax ; rdx = rax mulx rax, rax, rsi ; rax = (rsi * rdx) imul rcx, rsi ; rcx = rcx * rsi add rax, rcx ; rax = rax + rcx shr rax, 32 ; rax = rax >> 32 ret ; return rax example::DivisorU64U32B::rem::h85014351b09bd16a: mov rdx, rsi ; rdx = x mulx rcx, rcx, qword ptr [rdi] ; rcx = (M^+ * rdx) mov esi, dword ptr [rdi + 8] ; esi = D imul rcx, rsi ; rcx = rcx * rsi xor eax, eax ; eax = 0 sub rdx, rcx ; rdx = rdx - rcx cmovb eax, esi ; if borrow then ea add eax, edx ; eax = eax + edx ret ; return eax 左の手法では、 128bit 乗算を必要とするため、必要な乗算命令数の増加が性能に 悪影響を与える。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 35 / 56
36. 64bit整数の除算器 /// returns floor(a * b / 2^64) const fn mulh(a: u64, b: u64) -> u64 { ((a as u128 * b as u128) >> 64) as u64 } #[derive(Debug, Clone, Copy)] enum DivisorU64Kind { Shift { d: u64, s: u32 }, Cmp { d: u64 }, Up { d: u64, m: u64, s: u32 }, Even { d: u64, m: u64, s0: u32, s1: u32 }, Down { d: u64, m: u64, s: u32 }, Overflow { d: u64, m: u64, s: u32 }, } #[derive(Debug, Clone)] pub struct DivisorU64(DivisorU64Kind); impl DivisorU64 { pub const fn new<const USE_DOWN: bool>(d: u64) -> Self { assert!(d > 0); let s = d.ilog2(); if d.is_power_of_two() { Self(DivisorU64Kind::Shift { d, s }) } else if s == u64::BITS - 1 { Self(DivisorU64Kind::Cmp { d }) } else { let b = u128::MAX >> (u64::BITS - s); let m = (b / d as u128 + 1) as u64; if m.wrapping_mul(d).saturating_mul(u64::MAX / d) < Self(DivisorU64Kind::Up { d, m, s }) } else if d % 2 == 0 { let s0 = d.trailing_zeros(); let s1 = s - s0; Self(DivisorU64Kind::Even { d, m, s0, s1 }) } else if USE_DOWN { Self(DivisorU64Kind::Down { d, m: m - 1, s }) } else { let b = u128::MAX >> (u64::BITS - 1 - s); let m = (b / d as u128 + 1) as u64; Self(DivisorU64Kind::Overflow { d, m, s }) } } } /// return (x / d, x % d) pub const fn divrem(&self, x: u64) -> (u64, u64) { match self { &Self(DivisorU64Kind::Up { d, m, s }) => { let q = mulh(m, x) >> s; (q, x - d * q) } &Self(DivisorU64Kind::Even { d, m, s0, s1 }) => { let q = mulh(m, x >> s0) >> s1; (q, x - d * q) } &Self(DivisorU64Kind::Down { d, m, s }) => { let q = mulh(m, x.saturating_add(1)) >> s; (q, x - d * q) } &Self(DivisorU64Kind::Overflow { d, m, s }) => { let y = mulh(m, x); let q = (((x - y) >> 1) + y) >> s; (q, x - d * q) } &Self(DivisorU64Kind::Cmp { d }) => { let q = (x >= d) as u64; (q, x - d * q) } &Self(DivisorU64Kind::Shift { d, s }) => { (x >> s, x & (d - 1)) } } } } 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 36 / 56
37. DivisorU64 の正当性: 全体像 分岐 magic number 商 の式 条件 Shift 不要 Cmp 不要 Up Even Down Overflow 条件が重複するものは、性能が高い方を優先して選ぶとよい。この後、 Up , Even , Down , Overflow 分岐の正当性を与える。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [1] [2] [2] [1] は が真のとき 1、偽のとき 0 を返す アイバーソン括弧 (Iverson bracket) である。 [2] . が非2冪で Up 条件を満たさなければ、 Down 条件を満たす。 Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 37 / 56
38. DivisorU64 の正当性: Up 分岐 完全ブロック側の最大を見ると さらに より、 が非2冪であれば なので、最後のブロックも同じ条件で 抑えられる。 と から を用いて、 これで、 Up 条件 が、上側近似整除法の正確性を保証する十分条件として自然に現れる。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 38 / 56
39. DivisorU64 の正当性: Even 分岐 が非2冪偶数なら、奇数部分 を取り出せる。さらに とおくと、 から なので である。したがって、 元の偶数除数での除算は、右シフト後の奇数除数 での除算に帰着される。 また である。ここで とおくと、 かつ であるから、前節の系を , , , に適用できる。よって 以上より 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 39 / 56
40. DivisorU64 の正当性: Down 分岐 Down は を用いる。 である。したがって、全入力で を示せばよく、特に なので、十分条件として があればよい。 非 2 冪では であり、 Up が失敗すると なので よって Down 条件 が成り立つ。さらに Shift でも Up でもない場合は である。 である。 したがって、 のときだけ を に飽和しても、 得られる商は のままである。 [12] 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [12] ridiculous fish, “Labor of Division, Episode III: Faster Unsigned Division by Constants.” Oct. 2011, [Online]. Available: https://ridiculousfish.com/blog/posts/labor-of-division-episode-iii.html. [URL] Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 40 / 56
41. DivisorU64 の正当性: Overflow 分岐 が非2冪なら より また , とすると したがって では となり、上側近似の正確性条件より とすると から が従い、 である ことから、商 は 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 41 / 56
42. ac_library に見る実装の住み分け ac_library [13] の modint 実装は、法がいつ決まるかで方針が分かれる。 static_modint : 法がコンパイル時定数 コンパイラの定数除算最適化を使いやすく、前節の 64bit 整数除算器のうち mulh + shift 型の実装に近い dynamic_modint : 法が実行時に決まる Barrett Reduction が自然で、 とする上側近似の実装例に近い 補正回数が最悪 1 回必要 Barrett 側で補正が 1 回で足りる本質は、対象範囲で が成り立つことである。 たとえば , , なら、 は真の商より高々 1 大きい。 したがって を作ったあと、必要なら を 1 回だけ足せばよい。 このあと見る被除数128bitの話は、こうした の Barrett 型よりさらに一般の の fixed-width 実装へ進むものである。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム [13] AtCoder, “AtCoder Library.” [Online]. Available: https://github.com/atcoder/ac-library/. [URL] Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 42 / 56
43. 5. 被除数128bit・除数64bitの除算器 被除数128bit の除算はコンパイラ最適化されない? 被除数128bit・除数64bitの除算器 上位商・下位商 下側近似の一般形 実装方針の分岐 -bit divmod の代表的方針 系列の一般処理 mwf_lr , max_error_lower , max_error_upper の実装 max_error_lower が見ている量 max_error_lower の具体例 おおまかな設計指針 Rustでの実装例 性能比較 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 43 / 56
44. 被除数128bit の除算はコンパイラ最適化されない? gcc 16.1, clang22.1.0 での 128 ビット除算は、 コンパイラが内部で使用する 128 ビット除算を 計算するためのビルトイン関数 __udivmodti4 , __udivti3 などを呼び出し、除数が 定数であっても（2冪などの例外を除いた）大抵の 除数に対して最適化は抑止されています。 ただ、 clang trunk ビルドで試した所、 乗算命令等への置き換えをする最適化が 入るようになりました（右コード） 。 次期 llvm/clang に反映されるかもしれませんし、 gcc も追随するかもしれません。 #include <bits/stdc++.h> std::pair<std::uint64_t, std::uint64_t> divrem_1e19_u128(__uint128_t x) { static constexpr std::uint64_t D = 10'000'000'000'000'000 return {x / D, x % D}; } divrem_1e19_u128(unsigned __int128): movabs rax, 8507059173023461586 movabs r9, -7667109045778114189 mov rdx, rdi mulx r8, rcx, rax mulx r10, r10, r9 mov rdx, rsi mulx rax, r11, rax mulx rsi, rdx, r9 add r10, rcx movabs rcx, -8446744073709551616 adc r8, 0 add rdx, r10 adc rsi, r8 adc rax, 0 add rsi, r11 adc rax, 0 shld rax, rsi, 2 imul rcx, rax sub rdi, rcx mov rdx, rdi ret 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 44 / 56
45. 被除数128bit・除数64bitの除算器 (1) 上位商・下位商 とし、除数 を とする。 被除数を をそれぞれ -bit 非負整数とする。 商を をそれぞれ -bit 非負整数とする。 剰余を とする。 1. を で割ると商の上位 と仮の剰余 が得られる。 2. 仮の剰余 と被除数の下位 を結合して新しい被除数 を作る。 3. 新しい被除数を で割ると商の下位 と剰余 が得られる。 4. 以上の操作で、商 と剰余 が得られる。 最初の 1. の操作は既に紹介した 64bit 整数の除算器で実装できる。 以降では、 3. の操作に着目し、 被除数 が の範囲 かつ が 2 の冪ではない場合に ついて説明する。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 45 / 56
46. 被除数128bit・除数64bitの除算器 (2) 下側近似の一般形 , とし、 の範囲で を使って商を近似する。 ここで α > 1 を導入する理由は、被除数が 64bit レジスタに収まらない数を扱いたいからである。 α = 1 は被除数が 64bit 幅に収まるときに自然だが、128bit 被除数では、まず を取り出して 64bit 乗算器で扱いやすい大きさまで落としたくなる。 実装の基本手順は共通である。 1. をシフトで求める 2. を乗算とシフトで求める 3. を作る 4. の間、 から を引き、商を 1 増やして補正する したがって、下側近似を用いた設計上の論点は 「 を何 bit 幅で保持するか」 「 を何 bit 幅で保持するか」 「補正回数を 何回まで見込むか」 の 3 点に集約される。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 46 / 56
47. 128bit/64bit の除算: 実装方針の分岐 の -bit divmod は、除数 の形で大きく 3 通りに分かれる。 なら近似は不要で、商は x >> t 、剰余は x & (D - 1) で直接求まる が 2 の冪でない場合、まず候補になるのは または このスライドでは後者 を中心に見る 理由は、 を 以下に取るので となり、 magic number が、 が 2 の累乗数でないとき を満たすからである。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 47 / 56
48. -bit divmod の代表的方針 とした時の、 に対する の最大値(=補正回数の上界)を、代表的な方針ごとにまとめると次のようになる。 方針 の幅 の幅 補正回数 備考 は 2 の冪 不要 不要 シフト演算で実装可能 は 2 の冪でない, 約80% の で は 2 の冪でない, 約40% の で は 2 の冪でない, 約80% の で は 2 の冪でない, の の約70%がこの条件に該当 注) による補正回数の割合については、全数探索ではなく、 をランダムにサンプリングした評価による。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 48 / 56
49. mwf_lr , max_error_lower , max_error_upper の実装 個別の に対し、 の範囲における 誤差の最大値を厳密に求めたい場合、右のような Python 実装で計算できる。 mwf_lr(l, r, m, a, b, c, d) 下側最大誤差 max_error_lower(d, a, b, n) 上側最大誤差 max_error_upper(d, a, b, n) def mwf_lr(l: int, r: int, m: int, a: int, b: int, c: int, d: in assert l < r and m > 0 n = r - l d += c * l s = a * l t = s + b * (d // m) while True: s += b * (d // m) a += b * (c // m) c %= m d %= m ymax = (c * (n - 1) + d) // m t = max(t, s, s + a * (n - 1) + b * ymax) if ymax == 0 or a * b >= 0: return t if a < 0: s += a + b n, m, a, b, c, d = ymax, c, b, a, m, m - d - 1 def max_error_lower(d: int, a: int, b: int, n: int) -> int: assert d > 0 and a > 0 and b > 0 and n > 0 m, r = (a * b) // d, (n + d - 1) // d return (mwf_lr(0, r, a, b, -m, d, 0) + b - 1) // b def max_error_upper(d: int, a: int, b: int, n: int) -> int: assert d > 0 and a > 0 and b > 0 and n > 0 if n == 1: return 0 m, r = (a * b + d - 1) // d, (n - 1 + a - 1) // a return max(0, m + mwf_lr(0, r, d, m, -b, a, 1)) // b 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 49 / 56
50. max_error_lower が見ている量 とおくと、下側誤差の最大は に帰着できる。 とおくと、 mwf_lr(0, r, a, b, -m, d, 0) はちょうど を返すので、最後に で天井除算すれば下側最大誤差になる。 被除数 ・除数 の用途では を入れればよい。 一方、具体的な定数除数については n をその用途の被除数上限に合わせてそのまま入れればよい。たとえばまとめで挙げる の特化ケースも、対応する に対する max_error_lower の評価から、それぞれの補正回数上界が読める。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 50 / 56
51. max_error_lower の具体例 max_error_lower は解析用の関数なので、 や が 64bit に収まらない場合も、そのまま整数演算として評価できる。 例えば次の通りである。 max_error_lower(10**19, 2**63, 2**64, 2**128) == 2 max_error_lower(10**32, 2**64, 2**64, 2**128) == 1 max_error_lower(10**16, 2**52, 2**64, 10**32) == 1 したがって、対応する 商の真値 と 近似商 と その差 この差の最大値、つまり補正回数上界は、それぞれ では では では と読める。 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 51 / 56
52. おおまかな設計指針 が 2 の冪ではない、 系列を使うときの流れは次の 2 段階でよい。 1. 一般論としては が 2 の累乗数でないので だが、 は bit になりうる。補正回数はまず とみなす 2. 実際の除数 が決まったら 計算量が問題にならない範囲で max_error_lower(D, α, β, 2^W D) を評価し、その値をその まま必要補正回数の上限として採用する def conceptual_divmod(x: int, d: int) -> (int, int): assert 1 <= d < (1 << 64) and 0 <= x < (d << 64) s = (d.bit_length() - 1) # = floor(log2(d)) if d & (d - 1) == 0: # d is a power of two q, r = x >> s, x & (d - 1) else: # general case m = (1 << (64 + s)) // d # = floor(2^(64+s)/d) < 2^64 q = (m * (x >> s)) >> 64 # Note: (x >> s) < 2^65 r = x - q * d while r >= d: // correction loop, max 3 iterations q += 1 r -= d return q, r 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 52 / 56
53. 被除数128bit・除数64bitの除算器 (3) Rustでの実装例 入力 は の範囲であると仮定し、 (x / 1e19, x % 1e19) を計算する。 /// returns (x / 1e19, x % 1e19) for x < 2^64 * 1e19 pub fn divrem_1e19_rep2(x: u128) -> (u64, u64) { debug_assert!((x >> 64) < 10u128.pow(19)); const D: u64 = 10u64.pow(19); const M: u64 = (u128::MAX / 2 / D as u128) as u64; let mut q = (((x >> 127) > 0) as u64) * M + ((((x >> 63) as u64 as u128) * (M as u128) >> 64 let r = x - (q as u128 * (D as u128)); let f = (r >> 64) as u64 > 0 || (r as u64) >= D; let r = (r as u64).wrapping_sub((f as u64) * D); q += f as u64; let t = r.checked_sub(D); let r = t.unwrap_or(r); q += t.is_some() as u64; (q, r) } example::divrem_1e19_rep2::h8c930d3d88b36963: mov rax, rsi shr rax, 63 mov rdx, rsi shld rdx, rdi, 1 movabs rcx, -1432625727662628443 mulx rdx, rdx, rcx imul rax, rcx add rdx, rax movabs rcx, 8446744073709551616 mulx r8, rax, rcx sub r8, rdx add rax, rdi adc r8, rsi setne sil movabs rdi, -8446744073709551617 cmp rax, rdi seta r8b or r8b, sil movzx esi, r8b test sil, sil lea r8, [rax + rcx] cmove r8, rax add rcx, r8 xor eax, eax cmp r8, rdi seta al cmovbe rcx, r8 add rdx, rsi add rax, rdx mov rdx, rcx ret 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 53 / 56
54. 被除数128bit・除数64bitの除算器 (4) 性能比較 example::impl_divrem_1e19_asm::hc9f92614b40cf push rax mov rdx, rsi mov rax, rdi movabs rcx, -8446744073709551616 div rcx pop rcx ret DIV 命令を使った実装 294,881 ns/iter 測定は AMD Ryzen 5 7600X 6-Core Processor にて。 の範囲で被除数を ランダムにサンプリングした 131072 個の入力の処理時間を評価。 divrem_1e19_u128(unsigned __int128): movabs rax, 8507059173023461586 movabs r9, -7667109045778114189 mov rdx, rdi mulx r8, rcx, rax mulx r10, r10, r9 mov rdx, rsi mulx rax, r11, rax mulx rsi, rdx, r9 add r10, rcx movabs rcx, -8446744073709551616 adc r8, 0 add rdx, r10 adc rsi, r8 adc rax, 0 add rsi, r11 adc rax, 0 shld rax, rsi, 2 imul rcx, rax sub rdi, rcx mov rdx, rdi ret clang trunk による最適化実装 226,113 ns/iter example::divrem_1e19_rep2::h8c930d3d88b36963: mov rax, rsi shr rax, 63 mov rdx, rsi shld rdx, rdi, 1 movabs rcx, -1432625727662628443 mulx rdx, rdx, rcx imul rax, rcx add rdx, rax movabs rcx, 8446744073709551616 mulx r8, rax, rcx sub r8, rdx add rax, rdi adc r8, rsi setne sil movabs rdi, -8446744073709551617 cmp rax, rdi seta r8b or r8b, sil movzx esi, r8b test sil, sil lea r8, [rax + rcx] cmove r8, rax add rcx, r8 xor eax, eax cmp r8, rdi seta al cmovbe rcx, r8 add rdx, rsi add rax, rdx mov rdx, rcx ret 近似商+最大2回補正の実装 211,356 ns/iter 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 54 / 56
55. まとめ: 商の求め方 被除数 除数 条件 商の近似 誤差 備考 除数が2冪 商が高々1 32bit 64/32bit:Barrett 64bit:Up 64bit:Even 64bit:Down 64bit:Overflow 128/64bit:一般 128/64bit:条件付 128bit:特化 128bit:特化 128bit:特化 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 55 / 56
56. References 定数除数整数除算の高速実装 / 近似商と補正による除算アルゴリズム Mizar, “No.2440 Accuracy of Integer Division Approximate Functions.” yukicoder, 2023, [Online]. Available: https://yukicoder.me/problems/no/2440. [URL] [1] Mizar, “No.3397 Max Weighted Floor of Linear.” yukicoder, 2025, [Online]. Available: https://yukicoder.me/problems/no/3397. [URL] [2] Mizar, “No.3398 Accuracy of Integer Division Approximate Function 2.” yukicoder, 2025, [Online]. Available: https://yukicoder.me/problems/no/3398. [URL] [3] P. L. Montgomery, “Modular multiplication without trial division,” Mathematics of Computation, vol. 44, no. 170, pp. 519–521, 1985. [4] J. Hurchalla, “An Improved Integer Modular Multiplicative Inverse (modulo 2 ).” 2022, [Online]. Available: https://arxiv.org/abs/2204.04342. [URL] [5] w Nachia, Mathenachia, “mod 偶数の乗算にモンゴメリ乗算を使いたい？.” [Online]. Available: https://www.mathenachia.blog/even-mod-montgomery-impl/. [URL] [6] P. Barrett, “Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor,” in Advances in Cryptology — CRYPTO’ 86, 1987, pp. 311–323. [7] T. Granlund and P. L. Montgomery, “Division by invariant integers using multiplication,” in Proceedings of the ACM SIGPLAN 1994 Conference on Programming Language Design and Implementation, 1994, pp. 61–72. [8] Jr. Warren Henry S., Hacker’s Delight, 2nd ed. Addison-Wesley, 2013. [9] 著ヘンリー・S・ウォーレン、ジュニア and 滝沢徹, ハッカーのたのしみ : 本物のプログラマはいかにして問題を解くか, オンデマンド版. エスアイビー・アクセ ス, 2014. [URL] [10] D. Lemire, O. Kaser, and N. Kurz, “Faster Remainder by Direct Computation: Applications to Compilers and Software Libraries.” 2019, [Online]. Available: https://arxiv.org/abs/1902.01961. [URL] [11] ridiculous fish, “Labor of Division, Episode III: Faster Unsigned Division by Constants.” Oct. 2011, [Online]. Available: https://ridiculousfish.com/blog/posts/labor-of-division-episode-iii.html. [URL] [12] AtCoder, “AtCoder Library.” [Online]. Available: https://github.com/atcoder/ac-library/. [URL] [13] Mizar/みざー, 競プロキャンプ2026関東LT, 2026-06-06 56 / 56