Constant integer division faster than compiler-generated code

Kernel/VM探検隊@東京 No18
https://kernelvm.connpass.com/event/355100/
定数整数除算最適化再考
https://github.com/herumi/constdiv/

1. 定数整数除算最適化再考 光成滋生 @herumi 2025/08/09 Kernel/VM探検隊@東京 No18 1 / 17
2. 自己紹介&概要 @herumi サイボウズ・ラボで暗号や高速実装のR&D 2004年Google Open Source Peer Bonus受賞 Ethereum Foundation grant獲得x3 2016~2025 Microsoft MVPなど 定数整数除算最適化って何? uint32_t な整数定数 で割る関数 uint32_t udiv(uint32_t x) { return x/d; } を高速に求める話 Cの a/b はPythonの a//b で floor(a/b)) を示す 以下コード中を除いて本文で は厳密な分数、 を 整数への切り捨てを表す 2 / 17
3. さっそくですが簡単なクイズです クイズ1. 次のx64 Windows用アセンブリコードは何をする関数? func1: shr ecx, 3 mov eax, ecx ret 以下アセンブリコードをasmと略します 3 / 17
4. 答え 8で割る関数 3ビット論理右シフトと同じ x64 Windowsでは関数の第一引数はecxで戻り値はeax uint32_t func1(uint32_t x) { return x / 8; // return x >> 3; } ちなみにARM64なら第一引数と戻り値はw0 (32bit)なので func1: lsr w0, w0, #3 ret 4 / 17
5. ちょっと難しい クイズ2. 次のx64 Windows用asmは何をする関数? func2: mov eax, 0xaaaaaaab mul ecx shr edx, 1 mov eax, edx ret 5 / 17
6. ヒント Cで書くなら uint32_t func2(uint32_t x) { return (x * uint64_t(0xaaaaaaab)) >> 33; } 6 / 17
7. 答え 3で割る関数 uint32_t func2(uint32_t x) { return x / 3; } 1/3を近似する 定理 任意の について が成り立つ 7 / 17
8. もうちょい難しい クイズ3. 次のx64 Windows用asmは何をする関数? func3: mov eax, ecx imul rax, rax, 0x24924925 shr rax, 32 sub ecx, eax shr ecx, 1 add eax, ecx shr eax, 2 ret 8 / 17
9. ヒント Cで書くなら uint32_t func3(uint32_t x) { uint64_t v = x * uint64_t(0x24924925); v >>= 32; uint32_t t = x - v; t >>= 1; t += v; t >>= 2; return t; } 9 / 17
10. 答え 7で割る関数 uint32_t func3(uint32_t x) { return x / 7; } 1/7を近似する 疑問 のときと同じなのになぜコードが複雑になってる? この定数（シフトの数）はどうやって決まる? 10 / 17
11. ざっくりとした歴史 定数整数除算の最適化 定数整数で割る関数を乗算とシフトで近似するアイデアは1970年代からあった [GM94] Granlund, Montgomery, "Division by invariant integers using multiplication", SIGPLAN 1994 でコンパイラ用に実用化 Warren, "Hacker's Delight", 2002 で精密化 Lemire, et.al., "Integer division by constants: optimal bounds", 2021 で同値な定式化 定理 , , を で割った余りが となる を超えない最大の数 をとり、 , としたとき ならば全ての について . 定理を満たす を探せばOK 細かい話はblog記事定数除算最適化1参照 11 / 17
12. c L c H 1-bit x 32-bit c L x x × ＋ v v+x 7で割る関数の詳細 なぜ3で割る関数と違うのか? 定数 が32ビットを超えているから uint32_t な と掛けると65ビット を上位1ビット と下位32ビット に分ける このあと35ビット右シフトなので の下位32ビットは捨ててよい の上位32ビットを とする . は33ビットなので避ける . なので 細かい話はblog記事定数除算最適化2参照 12 / 17
13. 前振り終わり ここから本題 13 / 17
14. 現状のコンパイラを超えろ 私が使ってるCPUは64ビットだ [GM94]のテクニックは古くないかい?（30年前から変わってない） gcc 13.3.0 (x64), Apple clang 17.0.0 (ARM64), cl.exe 19.44.35213 (x64) で確認 x64は64bit×64bit=128bit乗算命令があるのでそれを使えばよいのでは 試したコード typedef __attribute__((mode(TI))) unsigned int uint128_t; uint32_t div7(uint3_t x) { uint64_t c = 0x124924925; uint32_t a = 35; uint128_t v = (x * uint128_t(c)) >> a; return uint32_t(v); } 14 / 17
15. ベンチマーク ループアンロールを防ぎたい 最適化は有効で測定する 関数ポインタで比較すると呼び出しコストの影響が大きい コンパイラの影響を除去するためXbyak(x64)/Xbyak_aarch64(ARM64)でJIT生成した 詳細はconstdiv参照 結果 CPU org (msec) mul64 (msec) rate=org/mul64 Xeon w9-3495X 129.237 112.496 1.15 M1 178.125 156.734 1.14 M4 92.944 75.298 1.23 1割ほど速くなった（M4速いな） 15 / 17
16. 改善 問題点 x64 128ビットのシフトshrdは遅い（Skylakeでimulと同じスループット1・レイテンシ3） ARM64 一度に128ビットを得る乗算は無い. 上位64ビット/下位64ビットを得る乗算を2回発行 改善 128ビット乗算はせずに のところだけ64ビットレジスタを使う uint32_t div7opti(uint32_t x) { uint64_t v = x * uint64_t(0x24924925); v >>= 32; v += x; return uint32_t(v >> 3); } 16 / 17
17. ベンチマーク2 更によくなった CPU org (msec) mul64 (msec) opti rate=org/opti Xeon w9-3495X 129.237 112.496 101.203 1.28 M1 178.125 156.734 142.202 1.25 M4 92.944 75.298 72.856 1.28 生成されたJITコード div7opti_x64: div7opti_aarch64: mov eax, 0x24924925 mov w9, #0x4925 imul rax, rdi movk w9, #0x2492, lsl #16 shr rax, 0x20 umull x9, w0, w9 add rax, rdi add x0, x0, x9, lsr #32 shr rax, 0x3 lsr x0, x0, #3 ret ret 17 / 17