定数整数除算・剰余算最適化再考

Constant Integer Division and Modulo Optimization Revisited
SCIS2026 https://www.iwsec.org/scis/2026/ 3B1-1発表資料

Transcript

1. 定数整数除算・剰余算最適化再考 サイボウズ・ラボ 光成滋生 2026/1/ SCIS2026 1 / 17
2. 本発表の概要 32ビット定数除算のコンパイラ生成コードよりもよい手法の提案 uint32_t divd(uint32_t x) { return x / 7; } などのCコードを x86-64, Apple M4用GCC, Clang, MSVCが生成するコードより36~39%高速な手法を提案 定数剰余算における, より小さい乗数の提案 MLKEMなどで使われる に対して既存のBarrett Reductionによる剰余算よりも 小さい乗数を使う手法の提案 select命令の追加が必要 2 / 17
3. 準備 記法 : 被除数の最大値, : 定数除数, : 被除数 商: , 余り: のとき : で割って余りが となる 以下の最大の整数 : 条件 cond が真なら x , 偽なら y を返す選択関数 u32 , u64 など: 符号無し32/64ビット定数 の場合わけ が2のべき乗 のとき として のとき それ以外 → 本発表の本題 3 / 17
4. GM法の定式化 定理1 について , とするとき ならば 任意の に対して が成り立つ 後述の定理2のためにオリジナルと同値な異なる定式化 の場合 となる が存在する の場合 , 共に u32 変数で乗算結果 は u64 に収まるので でよい の場合 GM法では中間変数の値が32ビットに収まる方法を提案（次のスライド） 4 / 17
5. x/7 の状況 Apple/x64用コンパイラの生成コード div7: ; Apple Clangの生成コード ; x64用gcc-14の生成コード mov w8, #18725 ; =0x4925 ; mov eax, edi movk w8, #9362, lsl #16 umull x8, w0, w8 ; x8 = x * cL ; imul rax, rax, 613566757 lsr x8, x8, #32 ; y = x8 = (x * cL) >> 32 ; shr rax, 32 sub w9, w0, w8 ; w9 = x - y ; sub edi, eax add w8, w8, w9, lsr #1 ; t = ((x - y) >> 1) + y ; shr edi / add eax, edi lsr w0, w8, #2 ; t >> (a - 33) ; shr eax, 2 ret ; ret Cのコードに書き直したもの（ 、 ） u32 udivd(u32 x) { u32 cL = c & 0xffffffff; u32 y = (u64(x) * cL) >> 32; u32 t = ((x - y) >> 1) + y; return t >> (a - 33); } 5 / 17
6. 乗算後の中間変数を32ビットに収める がu64 に収まらない と3段階のシフトに分解 ( ) とすると 32ビット×32ビット→64ビット(mul32)を利用して を計算 とすると の計算 は u32 に収まらない を利用して c L 1 1-bit x 32-bit x c L x x 2 32 + (x c L ) × ＋ c 6 / 17
7. 64ビットCPUの場合の改良案 中間変数を32ビットに収める必要はない 1. 64ビット×64ビット乗算(mul64)と128ビット整数シフトを利用する (mul64+sft) を直接計算する x64は imul の1命令で128ビットを得られる M4 (AArch64) は umulh / mul の2命令で上位/下位64ビットを得る u32 udivd1(u32 x) { return (x * u128(c)) >> a; } 2. 前述のmul32を利用して を算出後, 64ビット整数で計算する u32 udivd2(u32 x) { u64 cL = c & 0xffffffff; u64 y = (x * cL) >> 32; return (x + y) >> (a - 32); // 32ビットのオーバーフロー処理が無いのでシンプルになる } 7 / 17
8. ベンチマーク方法 評価方法 対象が小さいため関数呼び出しでなくinline展開する 複数回処理させて差分を評価する lp=1, 2, 3は処理単位内で複数回処理した回数（詳細は予稿集参照） lpが増えたときの差分 lat. は概ね処理時間のレイテンシに相当 u32 loop(FuncType f) { u32 sum = 0; // このループは展開しない for (u32 x = 0; x < N; x++) { u32 t = x; // f(t) = t//d のコードをlp回埋め込む for (int i = 0; i < lp; i++) { sum += f(t); t += sum; } ... } return sum; } 8 / 17
9. ベンチマーク結果 Clangとudiv1, udiv2との比較 Xeon w9-3495X - 処理時間（単位: msec） 方式 lp=1 lp=2 lp=3 lat. clang 60.62 239.26 452.89 196.14 mul64+sft 42.61 219.96 408.95 183.17 mul32 52.03 195.53 363.37 155.67 Apple M4 Pro 方式 lp=1 lp=2 lp=3 lat. clang 40.94 267.47 511.61 235.34 mul64+sft 33.29 201.42 382.63 174.67 mul32 33.33 223.04 423.37 195.02 lat.=((lp=3)-(lp=2))/2 x86-64は128ビット整数シフト (shrd) が遅いため mul32のレイテンシが小さい M4は128ビット整数シフト (extr) が速いため mul64+sftのレイテンシが小さい 9 / 17 命令 レイテンシ スループット add/sub/or 1 0.25 shr 1 0.5 imul 3 1 shrd 3 1
10. 提案方式 考慮すべきこと 128ビット乗算（64ビット×64ビット→128ビット）を使いたい 128ビットシフトは避けたい 提案方式 乗数を大きくしてシフトを除去する ・ ・ ・ 128ビット整数は64ビット整数2個で表現されるため 上位64ビット整数を取り出すコストは0 x64/M4で値設定を除いてアセンブリ言語の乗算1命令になる 10 / 17 div7_optimized: mov edx, edi mov rax, 0x24924924a0000000 mulx rax, rax, rax ret
11. 提案方式とベンチマーク結果 前述の比較に追加 Xeon w9-3495X - 処理時間（単位: msec） 方式 lp=1 lp=2 lp=3 lat. clang 60.62 239.26 452.89 196.14 mul64+sft 42.61 219.96 408.95 183.17 mul32 52.03 195.53 363.37 155.67 mul64 43.61 158.61 283.94 120.16 Apple M4 Pro 方式 lp=1 lp=2 lp=3 lat. clang 40.94 267.47 511.61 235.34 mul64+sft 33.29 201.42 382.63 174.67 mul32 33.33 223.04 423.37 195.02 mul64 30.21 156.40 289.56 129.68 x86-64で39%, Apple M4で36%の高速化を達成 レイテンシは60~80%の高速化 11 / 17
12. 定数整数剰余 Barrett Reduction , とする よりも少し小さい を求めて が 以上ならば を減らす u32 barrett(u32 x) { u32 q = (x * c) >> a; u32 r = x - q * d; while (r >= d) r -= d; // 定数のとり方によってはループは最大2回実行される return r; } MLKEM768の実装mlkem-nativeでは , に対して s16 ref3329(s16 x) { const int d = 3329; int t = (x * 20159 + (1<<25)) >> 26; return x - t * d; } 12 / 17
13. 提案方式 定理2 について , , とするとき 「 かつ 」 ならば 全ての について が成り立つ． 「」 が定理1の「 」から緩和されて が小さくなり, 除算の誤差が高々1になる よりも高々1大きい を求めて が負なら を加える 実装 u32 な に対して は必ず32ビットに収まる（後述） u32 uremd_mywork(u32 x) { u32 q = (x * c) >> a; u32 r = x - q * d; if (r < 0) r += d; return r; } 13 / 17
14. はGM法やBarrett Reductionより小さくなる のいくつかの例と分布 (定理1) (定理1) (定理2) (定理2) 7 0x124924925 35 0x24924925 32 824480341 0x14d653545 62 3 31 1239864366 0x1bb6663f9 63 4 32 のとき に対する の 89%が 特に で84% を占める のときは と最適化できる も乗算をビットシフト+加減算1回に置き換え可能 c 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16以上 合計 割合(%) 0.0 0.0 0.0 48.7 5.8 19.8 0.0 3.5 0.0 6.6 0.0 3.1 0.0 1.3 0.0 0.4 11.0 100.0 14 / 17
15. MLKEMの素数 ( ) の場合 Barrett Reductionによる方法（再掲）とMontgomery Reduction const int q = 3329; s16 ref3329(s16 x) { // ((q-1)/2)^2 * c > 2^32 int t = (x * 20159 + (1<<25)) >> 26; return x - t * q; } int mont(int x, int y) { int t = x * y; int m = ((t & 0xffff) * 3327) & 0xffff; t = (t + m * q) >> 16; if (t >= q) t -= q; return t; } 提案方式（定理2） を定理2に適用すると , となる（予稿集の +(1<<12) は不要だった） (x * (5<<2))>>16 とすれば16ビット乗算の上位16ビットを返すSIMD(PMULHW)が利用可能 なので組み込み系で乗算コストが軽減される可能性 s16 my3329(s16 x) { // ref3329の出力と完全一致 int t = (x * 5) >> 14; // (x*20)>>16でもOK t = x - t * q; if (t > q/2) t -= q; return t; // [-(q//2), q//2]に収めるため } 15 / 17
16. 定理2の証明 証明 について , に対して , とすると よって ならば または が成り立つ は非負定数で が動いたとき（連動して も動く）の を考える 最大値を与えるのは か より より よって条件 「」 が成り立てば となり が成り立つ M 0 x M d d-1 r (M,r 0 ) ... (M d ,d-1 ) 16 / 17
17. まとめ 定数整数除算最適化 u32 / u32除算を従来のGM法に対して64ビット整数乗算1回でできる手法を提案 既存のコンパイラGCC, Clang, MSVCに対してx86-64で39%, Apple M4で36%の高速化を達成 定数整数剰余算の改善 Barrett Reductionの亜種を提案 既存手法に比べて乗数が小さくなる手法を提案 64ビット乗算を避けられるためSIMD環境に適用しやすい MLKEMで使われる素数に対して剰余算が軽減される可能性 今後の課題 既存コンパイラへの本手法の提案 SIMD環境などでの評価 64~320ビット整数や符号つき整数に対する詳細な評価 17 / 17