除算・剰余算の高速化

除算・剰余算の高速化のテクニック。

$ \newcommand{\floor}[1]{\left\lfloor{#1}\right\rfloor}
\newcommand{\ceil}[1]{\left\lceil{#1}\right\rceil} $CPU での除算や剰余算は、加算や乗算に比べるとかなり遅い。ただ、多くの場合は高速化することができる。

ここでは、よくある高速化の方法

• 逆数乗算
• Exact Division
• Montgomery 乗算
• Div2by1

についてまとめる。

逆数乗算（Barrett reduction に近い）

除算 を 乗算 + シフト演算 に変換して計算する方法。$k$-bit CPU で $(k/2-1)$-bit の乗算・除算・剰余算を行うときに使える。

（64-bit 環境で）被除数が 64-bit で除数が定数（例えば、$10^9 + 7$）の場合、gcc はこの種の最適化を行う。Montgomery 乗算の恩恵をあまり得られない場合の理由の 1 つとして、この最適化がある。

内容

$b$ をレジスタの bit 幅、$2 \le M < 2^b$ を除数とし、$s := \floor{\log_2{(M - 1)}}$ とする。また、$X := \ceil{\frac{2^{b+s}}{M}}$ とおく。このとき、$X < 2^{b}$ であり、入力 $0 \le N < 2^{b-1}$ に対して、
$$
\floor{\frac{NX}{2^{b+s}}} = \floor{\frac{N}{M}}
$$
が成立する。

証明

まず、
$$
2^s \le M - 1 < 2^{s+1},\,2^{s} < M \le 2^{s+1}
$$
が成立する。また、ある $0 \le r < M$ が存在して、$X = \frac{2^{b+s} +r}{M}$ とかける。

このとき、
$$
X \le \frac{2^b(M-1) + r}{M} < \frac{2^b M}{M} = 2^b
$$
である。後者については、
$$
0 \le \frac{NX}{2^{b+s}} - \frac{N}{M} < \frac{1}{M}
$$
が示せればよいが、最初の不等号は $X$ の定義から成り立ち、2 つ目の不等号は
$$
\begin{align*}
\left(\frac{N (2^{b+s} +r)}{M \cdot 2^{b+s}} - \frac{N}{M}\right) - \frac{1}{M} &= \frac{N \cdot r}{M \cdot 2^{b+s}} - \frac{1}{M} \\
&= \frac{N \cdot r - 2^{b+s}}{M \cdot 2^{b+s}} \\
&< 0
\end{align*}
$$
により示せる（$r < M \le 2^{s+1}$, $N < 2^{b-1}$ を使う）。

実装

例えば、以下のように実装できる。$(10^8 + 6)! \pmod{10^8 + 7}$ の計算で、大体 3.1 倍程度早くなる（https://ideone.com/dpkoJW, 1.332 秒 / 0.435 秒）。$M = 1$ も動作するようにすると、少し遅くなる（2.5 倍程度）。

#include <cstdio>
#include <iostream>

using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;

/*
// 1 でも動作する
struct FastDiv {
  FastDiv() {}
  FastDiv(u64 n) : m(n) {
    s = (n == 1) ? 0 : 64 + std::__lg(n - 1);
    x = ((__uint128_t(1) << s) + n - 1) / n;
  }
  friend u64 operator / (u64 n, FastDiv d) { return __uint128_t(n) * d.x >> d.s; }
  friend u64 operator % (u64 n, FastDiv d) { return n - n / d * d.m; }
  u64 m, x; int s;
};
*/

// 1 は動作しない
struct FastDiv {
  FastDiv() {}
  FastDiv(u64 n) : m(n) {
    s = std::__lg(n - 1);
    x = ((__uint128_t(1) << (s + 64)) + n - 1) / n;
  }
  friend u64 operator / (u64 n, FastDiv d) { 
    return (__uint128_t(n) * d.x >> d.s) >> 64; 
  }
  friend u64 operator % (u64 n, FastDiv d) { return n - n / d * d.m; }
  u64 m, x; int s;
};

inline u64 mod64_32_small(u64 a, u32 b) {
  u32 q, r;
  __asm__ (
    "divl\t%4"
    : "=a"(q), "=d"(r)
    : "0"(u32(a)), "1"(u32(a >> 32)), "rm"(b)
  );
  return r;
}

int fact_slow(int n, int mod) {
  int ret = 1;
  for (int i = 1; i <= n; ++i) ret = i64(ret) * i % mod; // idivq になる
  return ret;
}

int fact_slow_asm(int n, int mod) {
  int ret = 1;
  for (int i = 1; i <= n; ++i) ret = mod64_32_small(u64(ret) * i, mod); // divl になる
  return ret;
}

int fact_fast(int n, int mod) {
  auto fd = FastDiv(mod);
  int ret = 1;
  for (int i = 1; i <= n; ++i) ret = i64(ret) * i % fd;
  return ret;
}

int main() {
  int mod; scanf("%d", &mod);
  clock_t beg = clock();
  int a1 = fact_slow(mod - 1, mod);
  clock_t t1 = clock();
  int a2 = fact_slow_asm(mod - 1, mod);
  clock_t t2 = clock();
  int a3 = fact_fast(mod - 1, mod);
  clock_t t3 = clock();
  printf("%d %.3f\n", a1, double(t1 - beg) / CLOCKS_PER_SEC);
  printf("%d %.3f\n", a2, double(t2 - t1) / CLOCKS_PER_SEC);
  printf("%d %.3f\n", a3, double(t3 - t2) / CLOCKS_PER_SEC);
  return 0;
}

長所と短所

長所としては、

• modulo ${M}$ の箇所の書き換えがほとんど不要
• Montgomery 乗算のように変換後の値で計算しなくていい

などがある。

短所は、（実用上では）除数 $M$ が $(k/2-1)$-bit 程度に制限されてしまう点である。Montgomery 乗算は、ここが $(k-1)$-bit で済む。

その他の資料

https://gmplib.org/~tege/divcnst-pldi94.pdf にも色々載っている。

Exact Division

剰余算の結果が 0 かどうかのみを判別したい場合、または、$M$ での剰余算が 0 となるとわかっている場合に $M$ で除算した場合に使える方法。例えば、素因数分解の前処理（試し割り）に使える。

この最適化は、gcc (7.2.0) では行われないようである。

内容

$b$ をレジスタの bit 長、$M$ を奇数とし、$M^{'}$ を $M \cdot M^{'} \equiv 1 \pmod{2^b}$ となる数で、$1 \le M' < 2^b$ を満たすものとする。このとき、$0 \le N < 2^b$　に対して、
$$
N \equiv 0 \pmod{M} \Leftrightarrow \mathrm{rem}(N \cdot M^{'}, 2^b) \le \floor{\frac{2^{b}}{M}}
$$ が成り立つ。ここで、$\mathrm{rem}(a, b)$ は $a$ を $b$ で割った時の余りとする。

実装

例えば、以下のように $[10^{10}+1, 10^{10} + 2 \cdot 10^5]$ の範囲で試し割りで素数判定を行う場合、ideone 上（https://ideone.com/SWG5lj）で大体 8〜9 倍程度早くなる（4.152 秒 / 0.415 秒）。キャッシュ効率にも依存するので、メモリの走査範囲が大きくなると少し恩恵が減る。

#include <cstdio>
#include <iostream>
 
using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;
 
template <typename T>
struct ExactDiv {
  ExactDiv() {}
  ExactDiv(T n) : t(T(-1) / n), i(mul_inv(n)) {}
  T mul_inv(T n) {
    T x = n;
    for (int i = 0; i < 5; ++i) x *= 2 - n * x;
    return x;
  }
  friend T operator / (T n, ExactDiv d) { return n * d.i; };
  bool divide(T n) { return n * this->i <= this->t; }
  T t, i;
};
 
ExactDiv<u64> ediv[200000];
 
bool is_odd_prime_slow(i64 n) {
  for (i64 i = 3; i * i <= n; i += 2) {
    if (n % i == 0) return false;
  }
  return true;
}
 
bool is_odd_prime_fast(i64 n) {
  for (i64 i = 3; i * i <= n; i += 2) {
    if (ediv[i].divide(n)) return false;
  }
  return true;
}
 
int main() {
  for (int i = 1; i <= int(2e5); i += 2) ediv[i] = ExactDiv<u64>(i);
 
  i64 beg = 1e10 + 1;
  i64 end = 1e10 + 2e5;
 
  clock_t t0 = clock();
  int ans1 = 0;
  for (i64 i = beg; i < end; i += 2) ans1 += is_odd_prime_slow(i);
  clock_t t1 = clock();
  int ans2 = 0;
  for (i64 i = beg; i < end; i += 2) ans2 += is_odd_prime_fast(i);
  clock_t t2 = clock();
 
  printf("%d %.3f sec\n", ans1, double(t1 - t0) / CLOCKS_PER_SEC);
  printf("%d %.3f sec\n", ans2, double(t2 - t1) / CLOCKS_PER_SEC);
  return 0;
}

長所と短所

長所は、除算・剰余算 が 乗算（+比較）で済むようになる点。短所は、適用可能な場面が少ない点。

（実用性を重視した）素因数分解・素数判定の前処理、一部の篩系のアルゴリズムの高速化に適する。

Montgomery 乗算

$k$-bit CPU で $(k-1)$-bit の 奇数 の乗算 + 剰余算を効率よく行うための方法。

内容

Montgomery modular multiplication - Wikipedia や モンゴメリ乗算 - yukicoder などを参照できる。

効率の面で書くと、$y = (x + \mathrm{rem}(xM^\prime, 2^b)\cdot M) / 2^b$ を使って計算するのではなく、$y = (x - \mathrm{rem}(x M^{-1}, 2^b) \cdot M) / 2^b$ を使って計算する方が、carry の計算が不要で、かつ $y$ が負かどうかを判定するだけでよいため、若干高速である。

実装

例えば、$10^8! \pmod{10^{18} + 3}$ の計算だと ideone 上（https://ideone.com/lqWWEj）で 7 倍ほど早くなる（2.723 秒 / 0.392 秒）。

#include <cstdio>
#include <cassert>

#include <iostream>

using i64 = long long;
using u64 = unsigned long long;
using u128 = __uint128_t;

struct Mod64 {
  Mod64() : n_(0) {}
  Mod64(u64 n) : n_(init(n)) {}
  static u64 modulus() { return mod; }
  static u64 init(u64 w) { return reduce(u128(w) * r2); }
  static void set_mod(u64 m) {
    mod = m; assert(mod & 1);
    inv = m; for (int i = 0; i < 5; ++i) inv *= 2 - inv * m;
    r2 = -u128(m) % m;
  }
  static u64 reduce(u128 x) {
    u64 y = u64(x >> 64) - u64((u128(u64(x) * inv) * mod) >> 64);
    return i64(y) < 0 ? y + mod : y;
  }
  Mod64& operator += (Mod64 rhs) { n_ += rhs.n_ - mod; if (i64(n_) < 0) n_ += mod; return *this; }
  Mod64 operator + (Mod64 rhs) const { return Mod64(*this) += rhs; }
  Mod64& operator *= (Mod64 rhs) { n_ = reduce(u128(n_) * rhs.n_); return *this; }
  Mod64 operator * (Mod64 rhs) const { return Mod64(*this) *= rhs; }
  u64 get() const { return reduce(n_); }
  static u64 mod, inv, r2;
  u64 n_;
};
u64 Mod64::mod, Mod64::inv, Mod64::r2;

inline u64 mod128_64_small(u128 a, u64 b) {
  u64 q, r;
  __asm__ (
    "divq\t%4"
    : "=a"(q), "=d"(r)
    : "0"(u64(a)), "1"(u64(a >> 64)), "rm"(b)
  );
  return r;
}

u64 fact_mod_fast(int N, u64 mod) {
  Mod64::set_mod(mod);
  Mod64 ret = Mod64(1), one = ret, t = one;
  for (int i = 1; i <= N; ++i) {
    ret *= t;
    t += one;
  }
  return ret.get();
}

u64 fact_mod_slow(int N, u64 mod) {
  u64 ret = 1;
  for (int i = 1; i <= N; ++i) ret = mod128_64_small(u128(ret) * i, mod);
  return ret;
}

int main() {
  int N; u64 M; scanf("%d %llu", &N, &M);
  clock_t t0 = clock();
  u64 ans1 = fact_mod_slow(N, M);
  clock_t t1 = clock();
  u64 ans2 = fact_mod_fast(N, M);
  clock_t t2 = clock();
  printf("%llu %.3f sec\n", ans1, double(t1 - t0) / CLOCKS_PER_SEC);
  printf("%llu %.3f sec\n", ans2, double(t2 - t1) / CLOCKS_PER_SEC);
  return 0;
}

補足

• divq は商が $2^{64}$ 以上になると floating point exception を吐く。
• __int128_t, __uint128_t 系の剰余算 (%) は、環境によっては遅い実装になっていて、理想的な実装より 10 倍程度遅いことがある（手元の環境だとそう）ので注意。

長所と短所

長所は、$(k-1)$-bit までの乗算 + 剰余算を高速に行える点がある。短所は、除算には向いていない点、偶数が混じると中国剰余定理を使う必要がある点などがある。

PowMod, Miller-Rabin, Pollard Rho などのアルゴリズムと相性がいい。

Div2By1

$k$-bit CPU で $k$-bit の乗算・除算・剰余算を効率よく行うための方法。逆数乗算に近い。

内容

https://gmplib.org/~tege/division-paper.pdf を参照。

実装

略

長所と短所

長所は、除算が高速に行える、偶数でも問題ないなどがある。短所は、Montgomery 乗算より若干遅い点。

多倍長の数に対して 64-bit 除算を行う場合などに適する。

long double での近似

32-bit 環境で 62-bit の乗算 + 剰余算 ($0 \le a, b < M < 2^{62}$ で $ab \pmod{M}$) を高速に計算したい場合の方法。

long double を使って近似的に商 $q = ab/M$ を計算し、剰余を計算する。

実装

Matters Computational の p765（付近）に載っている方法。
64-bit 環境で使っても、ideone 上では __uint128_t を介するより高速 https://ideone.com/gznOtF （2.755 秒 / 1.398 秒）。

#include <cstdio>
#include <ctime>

using i64 = long long;
using u64 = unsigned long long;
using f80 = long double;

inline u64 mul_mod(u64 a, u64 b, u64 mod) {
  u64 y = u64(f80(a) * b / mod + 0.5);
  u64 r = a * b - y * mod;
  return i64(r) < 0 ? r + mod : r;
}

int main() {
  u64 N, M; scanf("%llu %llu\n", &N, &M);
  clock_t t0 = clock();
  u64 ans1 = 1;
  for (u64 i = 1; i <= N; ++i) ans1 = mul_mod(ans1, i, M);
  clock_t t1 = clock();
  u64 ans2 = 1;
  for (u64 i = 1; i <= N; ++i) ans2 = __uint128_t(ans2) * i % M;
  clock_t t2 = clock();
  printf("%llu %.3f sec\n", ans1, double(t1 - t0) / CLOCKS_PER_SEC);
  printf("%llu %.3f sec\n", ans2, double(t2 - t1) / CLOCKS_PER_SEC);
  return 0;
}

長所と短所

32-bit でも容易に実装できるのが良い点。

短所は、結果を信じるのが難しい (?)