Montgomery 64

components/montgomery.hpp · GitHub で見る

共通ページ

#include "lib/components/montgomery.hpp"

Montgomery64 m(1000000007);
auto a = m.in(12), b = m.in(34);
auto value = m.reduce(m.mul(a, b));

操作	内容・返り値	計算量
Montgomery64(mod)	64 bit の奇数法を指定して構築	\(O(1)\)
in(x)	通常の整数を Montgomery 表現に変換	\(O(1)\)
reduce(x)	Montgomery 縮約。reduce(in(a)) で通常の整数に戻す	\(O(1)\)
add(a,b) / mul(a,b)	Montgomery 表現で加算／乗算	\(O(1)\)
pow(a,k)	Montgomery 表現で非負整数乗	\(O(\log k)\)

法は 1 より大きい 64 bit 符号なし奇数です。in(x) で Montgomery 表現に変換し、reduce(x) で戻します。add(a,b)・mul(a,b)・pow(a,k) は Montgomery 表現で計算します。入力は法未満とします。

構築・変換・加算・乗算は \(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 MONTGOMERY_HPP
#define MONTGOMERY_HPP

#include <cassert>
#include <cstdint>

struct Montgomery64 {
	using U = std::uint64_t;
	using W = __uint128_t;
	U mod, inverse = 1, r2, one;
	explicit Montgomery64(U n) : mod(n), r2(-W(n) % n) {
		for (int repeat_index = 0; repeat_index < 6; ++repeat_index) inverse *= 2 - n * inverse;
		assert(n > 1 && (n & 1));
		one = reduce(r2);
	}
	U reduce(W x) const {
		U a = x >> 64, b = (W(U(x) * inverse) * mod) >> 64;
		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;
		for (; e; e >>= 1, a = mul(a, a))
			if (e & 1) r = mul(r, a);
		return r;
	}
};

#endif

Back to top page