Montgomery 32

components/montgomery32.hpp · GitHub で見る

共通ページ

#include "lib/components/montgomery32.hpp"

Montgomery32 m(1000000007);
auto x = m.in(12), y = m.in(34);
auto product = m.reduce(m.mul(x, y));

操作	内容・返り値	計算量
Montgomery32(mod)	\(1<\mathrm{mod}<2^{32}\) の奇数法で構築	\(O(1)\)
mod / inverse / r2 / one	法／法の \(2^{32}\) に対する逆元／\(2^{64}\bmod\mathrm{mod}\)／Montgomery 表現の 1	\(O(1)\)
in(x)	法未満の値を Montgomery 表現に変換	\(O(1)\)
reduce(x)	\(0\le x<2^{32}\mathrm{mod}\) を Montgomery 縮約。reduce(in(a)) で通常の値に戻す	\(O(1)\)
mul(a,b) / add(a,b)	法未満の Montgomery 表現で乗算／加算	\(O(1)\)
pow(a,k)	Montgomery 表現で非負整数乗。指数は uint32_t	\(O(\log k)\)

型 U は uint32_t、W は uint64_t です。表現は常に法未満です。演算途中の加算が 32 bit を超える法にも対応します。より大きい法には Montgomery64 を使います。

構築・変換・加算・乗算は \(O(1)\)、累乗は \(O(\log k)\) 時間、空間は \(O(1)\) です。

検証コード

• test/discrete_log.test.cpp
• test/factorize.test.cpp
• test/is_prime.test.cpp
• test/is_prime_light.test.cpp
• test/mod_pow_log.test.cpp
• test/mod_sqrt.test.cpp
• test/multiplicative_convolution.test.cpp
• test/multiplicative_convolution_light.test.cpp

コード

#ifndef MONTGOMERY32_HPP
#define MONTGOMERY32_HPP

#include <cassert>
#include <cstdint>

struct Montgomery32 {
	using U = std::uint32_t;
	using W = std::uint64_t;
	U mod, inverse = 1, r2, one;

	explicit Montgomery32(U n) : mod(n), r2(-W(n) % n) {
		for (int i = 0; i < 5; ++i) inverse *= 2 - n * inverse;
		assert(n > 1 && (n & 1));
		one = reduce(r2);
	}
	U reduce(W x) const {
		U a = x >> 32, b = (W(U(x) * inverse) * mod) >> 32;
		return a < b ? a - b + mod : a - b;
	}
	U in(U x) const { return reduce(W(x) * r2); }
	U mul(U a, U b) const { return reduce(W(a) * b); }
	U add(U a, U b) const {
		U c = a + b;
		return c < a || c >= mod ? c - mod : c;
	}
	U pow(U a, U e) const {
		U r = one;
		if (mod < (U(1) << 30)) {
			U inverse = -this->inverse;
			auto product = [&](U a, U b) {
				W x = W(a) * b;
				return U((x + W(U(x) * inverse) * mod) >> 32);
			};
			for (; e; e >>= 1, a = product(a, a)) if (e & 1) r = product(r, a);
			return r < mod ? r : r - mod;
		}
		for (; e; e >>= 1, a = mul(a, a)) if (e & 1) r = mul(r, a);
		return r;
	}
};

#endif

Back to top page