除算・剰余算の高速化
除算・剰余算の高速化のテクニック。
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 := \left\lfloor{\log_2{(M - 1)}}\right\rfloor$ とする。また、$X := \left\lceil{\frac{2^{b+s}}{M}}\right\rceil$ とおく。このとき、$X < 2^{b}$ であり、入力 $0 \le N < 2^{b-1}$ に対して、$$\left\lfloor{\frac{NX}{2^{b+s}}}\right\rfloor = \left\lfloor{\frac{N}{M}}\right\rfloor$$ が成立する。
証明
まず、$$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{aligned} \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{aligned} $$ により示せる($r < M \le 2^{s+1}$, $N < 2^{b-1}$ を使う)。
実装
例えば、以下のように実装できる。$(10^8 + 6)! \bmod (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 \left\lfloor{\frac{2^{b}}{M}}\right\rfloor $$ が成り立つ。ここで、$\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 Wiki などを参照できる。
効率の面で書くと、$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! \bmod (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 の乗算・除算・剰余算を効率よく行うための方法。逆数乗算に近い。
実装
略
長所と短所
長所は、除算が高速に行える、偶数でも問題ないなどがある。短所は、Montgomery 乗算より若干遅い点。
多倍長の数に対して 64-bit 除算を行う場合などに適する。
long double での近似
32-bit 環境で 62-bit の乗算 + 剰余算 ($0 \le a, b < M < 2^{62}$ で $ab \bmod 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 でも容易に実装できるのが良い点。
短所は、結果を信じるのが難しい (?)