高速剰余算 div2by1 実装してみた

div2by1というアルゴリズムがある -> https://gmplib.org/~tege/division-paper.pdf

これはBarrett reductionやMontgomery乗算と違い、

  • (k, k) -> 2k-bitの乗算器でk-bitの剰余算ができる(Barrett reductionは(2k, 2k) -> 4k-bitの乗算器を必要とする)
  • modが偶数でも動作する(Montgomery乗算はmodが奇数の必要がある)

という二つの性質を持つ。今回は次の6種類のプログラムを実装し、$(8 \times 10^{7})! \bmod 998244353$を計算して速度を比較した。

  • 32-bit Montgomery
  • 32-bit div2by1
  • mod < $2^{30}$という制約を仮定し32-bit div2by1をちょっと改造したもの
  • これら三種のアルゴリズムをAVX2で高速化したもの

ベンチマークのコード、およびdiv2by1改造版のpythonコードについては、div2by1 改造版 · GitHub

手元(Ryzen 5 5600X)、およびAtCoderのコードテストでの実行時間は次の通り。

Local AtCoder
naive 285ms 482ms
naive(const mod) 227ms 298ms
montgomery 202ms 257ms
montgomery AVX 28ms 51ms
div2by1 316ms 458ms
div2by1 AVX 42ms 107ms
my div2by1 292ms 383ms
my div2by1 AVX 52ms 83ms

次のようなことがわかる

  • そもそもnaive(=除算命令を使っているはず)が直観よりかなり速い。おそらくIce Lakeから除算が速くなったというやつ(Intel Ice Lakeのプロセッサは整数除算命令がアツい - chroot("/home/hibari"))。でもIce lakeの製造開始は2019年らしい、老人さん?笑
  • montgomery + AVX2が速すぎる。4倍速以上になるとは思ってなかった。
  • div2by1 あんまり速くない。montgomery速すぎ + AVXで高速化とかするならもうmod 998244353決め打ちで問題なさそう であることを考えると正直出番がなさそうな