高速なMOD演算 Barrett Reduction, Montgomery 乗算

% を剰余演算子とする。 $M$ を2以上の2^31未満の整数とする。
$a \times b \% M$ をたくさん求めたいときに、 $M$ について前処理をしておけば除算・剰余算を使わずに計算できる場合がある。
とくに、 $M$ が定数であるがコンパイル時に決定できない場合に高速化することができる（コンパイル時に決定しているならばそもそも今回の手法と同程度に高速）。
2つの方法 Barrett Reduction と Montgomery Reduction を実装した。

まとめ

コード https://ideone.com/sNBQTz
実装は4種: ナイーブな剰余演算子、Montgomery乗算、Barrett Reduction、Barrett Reductionの128bit intを使わない実装。
CodeForces, AtCoder, ideoneで、10種類のMOD数で、Modの乗算を10^8 回ずつ (合計10^9 回) 計算したときの実行時間 ms。

	CF	AtC	io
%	9921	10000?	10000?
MR	6318	4953	4520
BR	-	5798	4650
BR64	8330	6695	6000?

? が付いているのは実行時間制限のため、計算の回数を分割して実行した。CodeForcesは128bit intが利用できないのでBRを実行していない。

Montgomery 乗算

自分の実装では最速。普通の剰余算の2倍の速さ。ただし、 $M$ が奇数でないと高速化が全くできない。

Barrett Reduction

Montgomeryには負ける。 $M$ がどんな値でもよい。実装が短い。使い方も剰余演算子を書くべきところでReduce関数を当てるだけなのでわかりやすい。
計算の過程で $M^3$ 程度の大きさになるので、 $M$ が32bit でも128bit intの計算が必要になる

Barrett Reduction 64

Barrett Reductionの計算を64bitで行ったもの。CodeForcesは (32bit環境のせいか乗算も遅いので) 性能が悪い。

Barrett Reduction

Barrett reduction - Wikipedia

0 以上 $M^2$ 未満の $x$ について $x \% M$ を求めたい。
実数 $s = 1 / M$ を用いて
$x \% M = x - \lfloor x s \rfloor M$
と表すことができる。 $s$ の近似値 $m / 2^k$ を使うと、 $x s \fallingdotseq x m / 2^k$ と表され、整数の乗算とビットシフトで計算することができる。

$2^k > M^2$ となるような $k$ を選び、 $m := \lfloor 2^k / M \rfloor$ を最初に一度だけ計算しておく。 $m / 2^k \le 1 / M$ より、 $x - \lfloor xm/2^k \rfloor M$ の値の範囲は $[ 0, 2M )$ 。

using ULL = unsigned long long;

struct BarrettReduction {
    ULL MOD;
    __int128 m;

    BarrettReduction(ULL MOD_): MOD(MOD_) {
	m = (__int128(1)<<64) / MOD;
    }

    ULL reduce(ULL x) const {
	ULL q = (x*m)>>64;
	x -= q * MOD;
	return x < MOD? x: x-MOD;
    }
};

void test() {
    BarrettReduction BR(1000000007);
    ULL a = BR.reduce(12345678ULL * 87654321ULL);
    printf("%llu\n", a);
}

Montgomery 乗算

モンゴメリ乗算 - Wikipedia

$R = 2^k$ を用いて、計算の過程で $x$ は常に $x R$ (の mod M) として保存されている。 $x$ から $x R$ 求めたり、逆に $x R$ から $x$ を求める計算が剰余算なしで $M^2$ の大きさまで計算できる。

使い方は、計算する前に実際の値をgenerate関数でMontgomery表現に変換する、乗算したら reduce する、Montgomery表現から元の値に戻したいときは reduce する。

struct MontgomeryReduction {
    ULL MOD;
    // MOD * NEG_INV % (2^32) == (2^32) - 1;
    ULL NEG_INV;
    // R2 == (2^64) % MOD;
    ULL R2;
    MontgomeryReduction(ULL MOD_): MOD(MOD_) {
	NEG_INV = 0;
	ULL s = 1, t = 0;
	REP (i, 32) {
	    if (~t & 1) {
		t += MOD;
		NEG_INV += s;
	    }
	    t /= 2;
	    s *= 2;
	}

	R2 = (1ULL<<32) % MOD;
	R2 = R2 * R2 % MOD;
    }

    // return x * R % MOD;
    ULL generate(ULL x) const {
	assert(x < MOD);
	return reduce(x * R2);
    }

    // return x / R % MOD;
    ULL reduce(ULL x) const {
	assert(x < MOD * MOD);
	x = (x + ((unsigned)x * (unsigned)NEG_INV) * MOD) >> 32;
	return x < MOD? x: x-MOD;
    }
};

void test() {
    MontgomeryReduction MR(1000000007);
    ULL a = MR.generate(12345678);
    ULL b = MR.generate(87654321);
    ULL c = MR.reduce(a * b);
    printf("%llu\n", MR.reduce(c));
}