定数時間 mod 演算の内部部品

components/mod_fast_core.hpp · GitHub で見る

共通ページ

ModInv と ModPow が使う内部部品です。通常は各ライブラリのヘッダーを include します。名前空間は mod_fast_detail です。

操作	内容・返り値	計算量
Reducer(m)	\(1\le m<2^{31}\) の法を指定し、逆数を前計算	\(O(1)\)
Reducer::reduce(x)	任意の uint64_t の剰余を uint32_t で返す	\(O(1)\)
Reducer::mul(a,b)	uint32_t の積の剰余	\(O(1)\)
Reducer::pow(a,e)	非負の uint64_t 指数による繰り返し二乗法	\(O(\log e)\)
Reducer::mod / Reducer::reciprocal	法／縮約用の逆数。参照用	\(O(1)\)
RationalTable(p)	\(2\le p<2^{31}\) の素数に対する有理数近似の表を構築	\(O(p^{2/3})\)
RationalTable::direct()	小さい法に対して全要素の表を使うか	\(O(1)\)
RationalTable::limit()	分子の絶対値の上限	\(O(1)\)
RationalTable::fraction(x)	\(x\equiv a/b\pmod p\) を満たす {a,b}。\(1\le x<p\)、\(a\ne0\)、\(b>0\)	\(O(1)\)

RationalTable は Farey 数列を走査し、各バケットの中央で分子の絶対値が小さい分数を保存します。分母上限とバケット幅を法に応じて調整し、クエリには右シフト・配列参照・積と差だけを使います。

Reducer の構築・縮約・積は \(O(1)\) 時間・空間、累乗は \(O(\log e)\) 時間です。RationalTable の構築は \(O(p^{2/3})\) 時間・空間、各クエリは \(O(1)\) 時間です。

検証コード

• test/mod_pow_log.test.cpp

コード

#ifndef MOD_FAST_CORE_HPP
#define MOD_FAST_CORE_HPP

#include <algorithm>
#include <cassert>
#include <bit>
#include <cstdint>
#include <utility>
#include <vector>

namespace mod_fast_detail {
struct Reducer {
	std::uint32_t mod;
	std::uint64_t reciprocal;
	explicit Reducer(std::uint32_t m) : mod(m), reciprocal(m > 1 ? std::uint64_t((__uint128_t(1) << 64) / m) : 0) {
		assert(1 <= m && m < (std::uint32_t(1) << 31));
	}
	std::uint32_t reduce(std::uint64_t x) const {
		if (mod == 1) return 0;
		std::uint64_t q = (__uint128_t(x) * reciprocal) >> 64;
		std::uint32_t r = x - q * mod;
		return r < mod ? r : r - mod;
	}
	std::uint32_t mul(std::uint32_t a, std::uint32_t b) const { return reduce(std::uint64_t(a) * b); }
	std::uint32_t pow(std::uint32_t a, std::uint64_t e) const {
		std::uint32_t r = mod == 1 ? 0 : 1;
		for (; e; e >>= 1, a = mul(a, a)) if (e & 1) r = mul(r, a);
		return r;
	}
};

class RationalTable {
	struct Fraction { std::uint16_t a, b; };
	std::vector<Fraction> fractions_;
	std::uint32_t p_, shift_ = 0, limit_;

  public:
	explicit RationalTable(std::uint32_t p) : p_(p), limit_(p - 1) {
		assert(2 <= p && p < (std::uint32_t(1) << 31));
		if (p <= 4096) return;
		std::uint32_t d = 1;
		while (std::uint64_t(d) * d * d < p) d *= 2;
		std::uint32_t width = std::bit_ceil((p + d * d - 1) / (d * d));
		shift_ = std::countr_zero(width);
		limit_ = (p + d - 1) / d + d * width / 2;
		assert(limit_ < p);
		fractions_.resize((p + width - 1) / width);
		std::uint32_t a = 0, b = 1, c = 1, e = d;
		for (std::uint32_t i = 0; i < fractions_.size(); ++i) {
			std::uint64_t x = std::min((2ULL * i + 1) * width, 2ULL * p - 1);
			while (2ULL * p * c < x * e) {
				std::uint32_t k = (d + b) / e, next_a = k * c - a, next_b = k * e - b;
				a = c, b = e, c = next_a, e = next_b;
			}
			if (x * b - 2ULL * p * a <= 2ULL * p * c - x * e)
				fractions_[i] = {std::uint16_t(a), std::uint16_t(b)};
			else fractions_[i] = {std::uint16_t(c), std::uint16_t(e)};
		}
	}
	bool direct() const { return fractions_.empty(); }
	std::uint32_t limit() const { return limit_; }
	std::pair<int, std::uint32_t> fraction(std::uint32_t x) const {
		assert(0 < x && x < p_);
		if (direct()) return {int(x), 1};
		auto [a, b] = fractions_[x >> shift_];
		std::int64_t t = std::int64_t(x) * b - std::int64_t(a) * p_;
		assert(t != 0 && -std::int64_t(limit_) <= t && t <= limit_);
		return {int(t), b};
	}
};
}

#endif

Back to top page