定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム

定数除数整数除算の高速化逆数近似による除算アルゴリズム Mizar/ みざー CS 集会Vket2026Summer 2026-07-14 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム
Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 1

今日のゴール主な対象は、CPU 命令やレジスタ幅まで意識して定数除数除算を高速化したい低レイヤ実装者。証明も用意したが、発表としては結論のみ用いる。定数除数による整数除算を、除算命令ではなく逆数近似の乗算を中心に実装する。この発表では、次の
3 本を軸に見る。被除数 32bit ・除数 32bit: 上側近似整除法で商と剰余を同時に復元する被除数 64bit ・除数 64bit: DivisorU64 の分岐で定数除数を処理する被除数 128bit ・除数 64bit: Barrett 型下側近似で最大誤差を評価する定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 2

除算命令はなぜ避けたいか整数除算命令 DIV は、乗算命令と比べて重い。 https://uops.info/ 例: x86_64 Arrow Lake-P の
64bit 整数演算では、おおまかに命令役割 Lat TP 傾向 DIV 64bit 除算 16-19 10.55 高レイテンシ・低スループット MUL , IMUL , MULX 64bit 乗算 3-5 1.00 除算より軽い SHRD シフト 4 1.00 SHRD は乗算と同程度 SHR シフト 1 0.33 TP0.33 ≒ 3 並列 ADD , SUB 加減算 1 0.17 TP0.17 ≒ 6 並列除数が固定なら、の固定小数点近似を前計算し、実行時の DIV を避けたい。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 3

定数除数除算の基本アイデアを整数型のビット幅とする。除数が定数なら、あらかじめとなる magic number を計算しておく。すると、商はのように近似できる。この発表で扱う問いは次の 2
種類である。どのような条件で近似商が真の商と一致するか？近似商と真の商との最大差はいくつか？定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 4

軸1 の導入: 32/32bit 上側近似最初の軸では、被除数も除数も 32bit の場合を扱う。上側近似を前計算しておくと、実行時には 64bit
乗算（128bit 積の上位64bit ・下位64bit の取り出し）を 2 回で商と剰余の両方を求められる。正の整数と非負整数がのときとおくと、定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 5

上側近似整除法の正確性条件非負整数と、正の整数による整除法で商と余りを得ることを考える。正の整数に対しとおくと、以下の必要十分条件がそれぞれ成り立つ。かつ定数除数整数除算の高速化
/ 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 6

上側近似整除法の正確性条件の証明非負整数の被除数と、正の整数の除数による整除法を考える。正の整数に対しとおくと、以下の必要十分条件がそれぞれ成り立つ。基本変形
とおく。よりここでとおくとさらに商の正確性また、(2) よりよって商と剰余の同時復元商が正しい、すなわちのとき、(1) よりしたがって、(2) からここでよりまたなので、ならであり、上の商条件より商も正しい。よって定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 7

Lemire-Bartlett-Kaser の Theorem 1/2 との対応と読むと、前ページの証明は [^lbk] の multiply-divide 形式を点ご
とに見たものになっている。 Theorem 1: 商だけ前ページの商条件は範囲で一様に成り立たせる最悪ケースはの最大ので、したがってこれは Theorem 1 の上界である。 Theorem 2: 商と剰余商が正しい上で剰余も復元する点ごとの条件はである。範囲で一様に成り立たせるにはで十分かつ必要であり、と同値である。これは Theorem 2 の上界である。前ページの条件は各に対する局所条件であり、 Theorem 1/2 はそれを全体へ一様化した最適境界と見られる。 [^lbk]: Daniel Lemire, Colin Bartlett, Owen Kaser, “Integer Division by Constants: Optimal Bounds,” arXiv:2012.12369. https://arxiv.org/abs/2012.12369 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 8

32-bit 同士の商と剰余の実装として、商と剰余を計算する実装例を示す。 #[derive(Clone, Copy, Debug)] pub struct
DivisorU32U32A { m: u64, d: u64 } impl DivisorU32U32A { pub const fn new(d: u32) -> Self { assert!(d > 1); Self{ m: (u64::MAX / d as u64).wrapping_add(1), d: d as u64 } } /// returns (x / d, x % d) #[target_feature(enable = "bmi2")] pub const fn divrem(&self, x: u32) -> (u32, u32) { let mx = self.m as u128 * x as u128; ((mx >> 64) as u32, (mx as u64 as u128 * self.d as u128 >> 64) as u32) } } divrem: mov edx, esi ; edx = x mulx rax, rdx, qword ptr [rdi] ; rax:rdx = rdx * M⁺ ; rax:rdx = x * M⁺ ; rax = floor(M⁺x / 2^64) = q ; rdx = M⁺x mod 2^64 mulx rdx, rdx, qword ptr [rdi + 8] ; rdx:_ = rdx * D ; rdx = floor(((M⁺x mod 2^64) * D) / 2^64) = r ret ; return (eax, edx) の場合はであるため、実装上の注意が必要である。特に 32bit 入力では scalar 実装では短い命令列になりやすい一方で、32bit 整数を多数並べて x86_64 SIMD を使う用途では、32bit lane を保ったままこの形を高速に処理しづらい。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 9

軸1 のまとめ: 32/32bit 上側近似 32bit 同士の定数除数除算では、とした上側近似を使うと、商と剰余を同じ積から取り出せる。この軸で使った核は、上側近似整除法の条件である。これは、商だけでなく剰余も同時に復元できることを保証する。
実装面では、とを保持すれば、除算命令なしで短い命令列に落ちる。ただし、この手法では SIMD で 32bit lane を保ったまま処理するのは難しい。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 10

軸2 の導入: 64/64bit DivisorU64 次の軸では、被除数も除数も 64bit の一般的な定数除数除算を扱う。 32bit での手法とは異なり、除数の性質に応じて分岐を選ぶ DivisorU64
の Rust 実装例を紹介する。分岐の選択は静的コンパイルや JIT コンパイルにより最適化することができる。 Shift : 除数が 2 冪なら右シフト Cmp : 除数が大きい場合は比較だけで商が決まる Up : 上側近似で直接計算 Even : 偶数除数を奇数部分へ帰着 Down : round-down algorithm Overflow : 1 bit 余分な精度を使うフォールバック定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 11

DivisorU64 実装例 /// 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) < m { 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/ みざー, CS 集会Vket2026Summer, 2026-07-14 12

DivisorU64 の正当性: 全体像分岐 magic number 商の式条件考え方
Shift 不要右シフトによる高速除算 Cmp 不要 [^1] 大きな除数の場合の簡易比較 Up [^2] 上側近似による除算 Even 偶数除数を右シフトして奇数除数へ帰着 Down [^2] 下側近似による除算 Overflow 1 bit 余分な精度を使う汎用フォールバック条件が重複するものは、性能が高い方を優先して選ぶとよい。は、がの約数でない事を示す。 [^1]: はが真のとき 1 、偽のとき 0 を返すアイバーソン括弧 (Iverson bracket) である。 [^2]: . が非2 冪で Up 条件を満たさなければ、 Down 条件を満たす。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 13

DivisorU64 の正当性: Up 分岐完全ブロック側の最大を見るとさらにより、が非2 冪であればなので、最後のブロックも同じ条件で抑えられ
る。とからを用いて、これで、 Up 条件が、上側近似整除法の正確性を保証する十分条件として自然に現れる。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 14

DivisorU64 の正当性: Even 分岐が非2 冪偶数なら、奇数部分を取り出せる。さらにとおくと、から
なのでである。したがって、元の偶数除数での除算は、右シフト後の奇数除数での除算に帰着される。またである。ここでとおくと、かつであるから、前節の系を , , , に適用できる。よって以上より定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 15

DivisorU64 の正当性: Down 分岐 Down はを用いる。である。したがって、全入力でを示せばよく、特になので、十分条件として
があればよい。非 2 冪ではであり、 Up が失敗するとなのでよって Down 条件が成り立つ。さらに Shift でも Up でもない場合はである。である。したがって、のときだけをに飽和しても、得られる商はのままである。 [^fish2011] [^fish2011]: ridiculous fish, “Labor of Division, Episode III: Faster Unsigned Division by Constants,” 2011. https://ridiculousfish.com/blog/posts/labor-of-division-episode-iii.html 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 16

DivisorU64 の正当性: Overflow 分岐が非2 冪ならよりまた , とすると
したがってではとなり、上側近似の正確性条件よりとするとからが従い、であることから、商は定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 17

関連研究: 直接剰余計算と round-down Lemire-Kaser-Kurz Lemire-Kaser-Kurz [^lemire2019] の Theorem 1 は、上側近
似逆数で商と剰余を同時に正しく復元できる条件を与えている。固定除数に対し、を前計算し、Algorithm 1 では積の下位側から剰余を直接求める。本発表の上側近似整除法における「商と剰余の同時復元」の点ごとの条件が、 Theorem 1 の条件に対応する。 [^lemire2019]: Daniel Lemire, Owen Kaser, Nathan Kurz, “Faster Remainder by Direct Computation,” arXiv:1902.01961. https://arxiv.org/abs/1902.01961 DivisorU64 との関係 Up : 上側近似で商を直接得る Even : 偶数除数を右シフトして奇数部分へ帰着 Overflow : bit の上側近似を bit 乗算で実装 Down : 上側近似ではなく round-down algorithm Overflow は前ページのように、 bit の上側近似を bit 乗算と加算・シフトへ変形して実装する。 Down 分岐は Lemire-Kaser-Kurz ではなく、とを使う ridiculous_fish の round-down algorithm [^fish2011] に対応する。 [^fish2011]: ridiculous fish, “Labor of Division, Episode III: Faster Unsigned Division by Constants,” 2011. https://ridiculousfish.com/blog/posts/labor-of-division- episode-iii.html 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 18

軸2 のまとめ: 64/64bit DivisorU64 DivisorU64 は、単一の magic number だけで押し切るのではなく、除数の形に応じて最も扱いやすい式へ分岐する。
分岐役割 Shift , Cmp 特殊な除数を乗算なしで処理する Up 上側近似がそのまま正しい場合に使う Even 偶数除数を右シフトして奇数除数へ帰着する Down 上側近似が選べない場合に下側近似で商を得る Overflow 1 bit 余分な精度を bit 乗算へ変形するこの軸の主眼は、条件分岐によって正確性と実装可能性を両立することである。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 19

軸3 の導入: 128/64bit Barrett 型下側近似最後の軸では、128bit 被除数を 64bit 定数除数で割る場合を扱う。ここでは、近似商が必ずしも真の商と一致するとは限らない。その代わり、下側近似で小さめ
の商を作り、剰余を見て補正する。見るべき問いは、近似商と真の商の差、すなわち商の最大誤差がどこまで抑えられるかである。非負整数について、では、つまり範囲最大も高々個別の除数では最大誤差を厳密に求められる定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 20

128/64bit: 概念実装なら商は未満であり、全体へ広げたい場合も筆算の要領で上位商・下位商の計算としてこの形を2 度適用すればよい。 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 が2冪の場合 q, r = x >> s, x & (d - 1) else: # d が非2冪の場合 m = 2 ** (64 + s) // d # = floor(2^(64+s)/d) < 2^64 q = (m * (x >> s)) >> 64 # (x >> s) は 65-bit になり得る r = x - q * d while r >= d: # 商の最大誤差3なので高々3回 q += 1 r -= d return q, r 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 21

Barrett 型下側近似の商の最大誤差は正の整数、は非負整数とし、とする。各行は、左列を満たす全ての被除数に対する商誤差の上界を示す。前頁の概念実装は、商の最大誤差が 3
以下になる、すなわちの場合に対応する。被除数の範囲・追加条件最大誤差上界仮の剰余一般条件: かつかつかつかつかつ定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 22

なら前提・定義は正の整数、は非負整数とし、を仮定する。このときである。主張
証明とおく。するとであり、よりである。まずとよりよってである。一方、とおく。よりこれにを掛けるとよりまたなのでは整数であるからまた、より以上よりしたがってである。最後に、かつなので、である。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 23

で最大誤差が高々2 になる追加条件前提・定義は正の整数、は非負整数とし、とする。さらにを仮定する。主張証明
まずとよりしたがってである。次にとおく。よりである。またなのでよりは単調なのでと書く。よりここでなので仮定よりであり、またよりである。したがってゆえには整数なので以上より。最後にから。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 24

最大誤差が高々1 になる追加条件前提・定義は正の整数、は非負整数とし、とする。このときである。さらに
を仮定する。主張証明とおくと、である。まずとよりしたがってである。ならであり、追加条件はとなる。よってである。以下ではとする。ならであり、よりである。残るの場合、よりであり、追加条件はと同値である。またかつなのでしたがってである。左辺は整数なので先ほどのと合わせて、以上の各場合で。最後にであり、 , より定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 25

かつなら商の最大誤差は0 前提・定義は正の整数、は非負整数とし、とする。さらにを仮定する。主張
証明と書ける。ただしは正整数である。よりなので、を消去してであり、と書ける。ここでも正整数である。したがってとおくと、である。なのでであり、なので、を加えても次のの倍数までは届かない。よってであり、定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 26

mwf , max_error_lower , max_error_upper の実装個別のに対し、の範囲における誤差の最大値を厳密に求めたい場
合、右のような Python 実装で計算できる。 mwf(l, r, m, a, b, c, d) 下側最大誤差 max_error_lower(d, a, b, n) 上側最大誤差 max_error_upper(d, a, b, n) def mwf(l: int, r: int, m: int, a: int, b: int, c: int, d: int) -> int: 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(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(0, r, d, m, -b, a, 1)) // b 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 27

max_error_lower と違反点 max_error_lower(D, α, β, n) は、範囲でとなる違反点があるかを判定できる。
とおく。最初に閾値を超える候補はの形に絞れる。実際、したがって max_error_lower を判定関数として使えば、最小の違反点を探索できる [^divapprox2] 。違反点の存在はで分かれる。なら、十分大きいで必ずとなるなら、として、違反点が存在する条件はつまりがで割り切れ、かつなら、指定したを超える違反点は存在しない。 [^divapprox2]: yukicoder No.3398 “Accuracy of Integer Division Approximate Functions 2” で問われる最小違反点探索に対応する。https://yukicoder.me/problems/no/3398 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 28

max_error_lower の具体例個別の定数除数では、最大誤差を前計算して補正上限に使える。これは、対応する真の商と近似商について、各被除数の範囲で最大誤差がになることを意味する。
定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 29

Rust 実装例: 入力範囲をに制限し、を計算する。この範囲では最大誤差が 2 以下なので、正しい商と剰余を得られる。
/// 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) as u64); 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/ みざー, CS 集会Vket2026Summer, 2026-07-14 30

性能比較 , 入力 example::impl_divrem_1e19_asm::hc9f92614b40cf972: push rax mov rdx, rsi mov
rax, rdi movabs rcx, -8446744073709551616 div rcx pop rcx ret DIV 命令を使った実装 294,881 ns/iter 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/ みざー, CS 集会Vket2026Summer, 2026-07-14 31

軸3 のまとめ: 128/64bit Barrett 型下側近似 128/64bit の定数除数除算では、下側近似で近似商を作ってから、残った剰余候補を補正する形が実装しやすい。下側近似なので、近似商は真の商を超えない。したがって、必要な補正上限はこの差の最大値、すなわち最大誤差で決まる。
では、という粗い上限を証明できる入力範囲やを絞ると、最大誤差を 2 や 1 へ落とせる個別定数では max_error_lower により、最大誤差を前計算できる実装例では、最大誤差 2 の近似商実装が DIV 命令より短い依存列で動く。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 32

まとめ定数除数整数除算は、逆数近似の乗算を中心に実装できる。本発表では、次の 3 本を軸に見た。 32/32bit 除算器: 上側近似整除法の正確性条件により、商と剰余を同時復元する
64/64bit 除算器: DivisorU64 の Up , Even , Down , Overflow 分岐で定数除数を処理する 128/64bit 除算器: 最大誤差を評価すると、 Barrett 型近似で近似商を求めた後、残った剰余候補を補正する形が実装しやすい定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 33

Accuracy of Integer Division Approximate Functions (2023 年8 月出題) https://yukicoder.me/problems/no/2440
個のクエリが与えられる。番目のクエリでは、整数が与えられる。かつとなるような整数がいくつあるかを数え上げ、とおく。この整数を答えよ。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 34

Accuracy of Integer Division Approximate Functions の解法の要点を求める。ならば、答えはである。
ならば、はの範囲となりうる上限を用いて以下のように表せる。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 35

Max Weighted Floor of Linear (2025 年12 月出題) https://yukicoder.me/problems/no/3397 個のテストケースがあります。
各テストケースについて、正の整数と整数が与えられます。を求めてください。つまり、整数を以上未満で動かしたとき、式が取りうる値の中で最大のものを求めてください。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 36

Max Weighted Floor of Linear の解法の要点求める値はである。は非減少な階段関数であり、となる各段では目的
関数はという一次関数になる。そのため各段では端点だけを見ればよく、なら右端、なら左端が候補になる。段の端点を floor 式で表すと、候補全体の最大化は同じ形の小問題に帰着できる。正規化後の再帰では分母がと Euclid 型に減るため、各テストケースをで処理できる。別解として、 floor_sum の再帰をモノイド積に一般化した floor_prod を使う方法もある。の増加と floor 値の増加をそれぞれ重みのステップと見なし、累積和と最大 prefix を持つモノイド上で計算すると、目的関数の最大値を同じ Euclid 型の反復で求められる。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 37

Accuracy of Integer Division Approximate Functions 2 (2025 年12 月出題)
https://yukicoder.me/problems/no/3398 各テストケースについて、正の整数と非負整数が与えられます。次に示す関数を定義します。を満たす非負整数が存在する場合は、を満たす最小の非負整数を答えてください。存在しない場合は、を答えてください。出力すべき値は、 64bit 整数で表記できる範囲を超えることがあります。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 38

Accuracy of Integer Division Approximate Functions 2 の解法の要点とおく。が
の倍数に達したときだけの値は増加しうるため、閾値を初めて超える候補を調べるには、まずの形に注目する。とした場合、誤差の閾値条件はと同値変形できる。この不等式を用いて、区間内にを満たす候補が存在するかを Max Weighted Floor of Linear で判定し、その判定を用いて最小のを二分探索/ 指数探索する。また不等式からは、がを満たすときはを満たす候補が存在することもわかるため、二分探索/ 指数探索の上限をに絞ることもできる。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 39

商の求め方: 整序を " " 、非整序を " " の記号で示す被除数除数
条件商の近似誤差備考除数が2 冪商が高々1 32bit 64/32bit:Barrett 64bit:Up 64bit:Even 64bit:Down 64bit:Overflow Barrett: 一般3 Barrett: 一般2 Barrett: 条件付 Barrett: 条件付 128bit: 特化 128bit: 特化 128bit: 特化定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 40

付録: 発表外資料 SRT 除算器の概要 Goldschmidt 除算の概要中国剰余定理・Garner のアルゴリズムと Montgomery reduction

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/ みざー, CS 集会Vket2026Summer, 2026-07-14 42

Goldschmidt 除算の分子と分母に共通の補正係数を掛けることを繰り返し、分母をに近づけていくのが Goldschmidt 除算の基本的な考え方である。例として、被除数と除数に共通の冪スケーリングを行い、除数を
またはの範囲に正規化済みであるとする。とおくとなので、はに収束する。したがって、と展開できる。実装上は、とおき、各段でを計算する。このときである。有限段で打ち切ると残る誤差はの部分に対応する。は各段で二乗されるため急速にに近づき、十分な段数の後ではがの近似値になる。定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 45

中国剰余定理・Garner のアルゴリズムと Montgomery reduction 互いに素な正整数に対し、とする。このとき、整数からを作ると、かつ
かつが成り立つ。これは 2 つの法に対する中国剰余定理の解を、片側から順に復元する形になっている。ここでとおく。とすれば、はかつを満たす。したがってはで割り切れ、とおくと、に対してとなる。これは Montgomery reduction [^1] の REDC の中核である。 Montgomery reduction に関連して、における乗法逆元を高速に計算する手法 [^2] や、偶数に対する Montgomery 乗算の手法 [^3] も存在する。 [^1]: P. L. Montgomery, “Modular multiplication without trial division,” Mathematics of Computation, vol. 44, no. 170, pp. 519– 521, 1985. [^2]: J. Hurchalla, “An Improved Integer Modular Multiplicative Inverse (modulo 2 ).” 2022, https://arxiv.org/abs/2204.04342. [^3]: Nachia, Mathenachia, “mod 偶数の乗算にモンゴメリ乗算を使いたい？.” https://www.mathenachia.blog/even-mod-montgomery-impl/. 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 46

Barrett Reduction Paul Barrett の CRYPTO’86 論文 [^1] は、RSA を標準的な
DSP 上で高速に実装するため、多倍長除算を避ける modulo reduction を提案した。基数の桁の整数で除算を行うとき、を事前計算しておき、を除算ではなく、乗算・桁の切り捨て・減算で求める。現在は Barrett reduction と呼ばれる、固定された除数に対する商・剰余計算を、事前計算した逆数近似・乗算・桁の切り捨て・補正減算で行う古典的手法である。 [^1]: 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/ みざー, CS 集会Vket2026Summer, 2026-07-14 48

Barrett Reduction の誤差評価前提・定義は整数、主張証明と書くと、よりまず
よりしたがって一方、かつよりよってまたよりしたがってよってここでなので、したがって定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 49

Granlund-Montgomery method 除数の逆数をとして近似し、がすべてので成り立つようにとを選ぶ方法が
[^1] で提案されました。論文の Theorem 4.2 では、次の十分条件を与えています。 [^1]: 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. 定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 50

Granlund-Montgomery method の正確性条件 Theorem 4.2: unsigned division を乗算で正確に計算する条件 bit
非負整数を定数で割ることを考える。主張を非負整数、とし、が成り立つとする。このとき、任意の整数についてつまり、をで近似しても、 bit 入力に対する商は変わらない。証明の要点とおくと、仮定より、ただしと書く。このとき右辺は非負であり、さらによりしたがってよって定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム Mizar/ みざー, CS 集会Vket2026Summer, 2026-07-14 51

定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム

定数除数整数除算の高速化 / 逆数近似による除算アルゴリズム

More Decks by Mizar

Other Decks in Technology

Featured

Transcript