一般化 CRT

math/crt.hpp · GitHub で見る

共通ページ

連立合同式 \(x\equiv r_i\pmod {m_i}\) を解きます。法が互いに素でなくても使え、剰余は負でも構いません。

#include "lib/math/crt.hpp"
#include <cassert>

const auto a = crt(2, 6, 5, 9);
assert(a && a->residue == 14 && a->modulus == 18);

const auto b = crt({2, 5, 14}, {6, 9, 18});
assert(b == a);
assert(!crt(1, 2, 0, 4));

操作・型	引数・返り値
crt(r0, m0, r1, m1)	std::int64_t の 2 組の剰余・法を渡す。std::optional<CRTSolution> を返す。constexpr 対応
crt(residues, moduli)	const std::vector<std::int64_t>& の剰余列と法の列を渡す。同じ添字の組を合同式とし、std::optional<CRTSolution> を返す
solution.residue	最小非負解。型は std::int64_t
solution.modulus	全解の周期で、全ての法の最小公倍数。型は std::int64_t

各法は \(1\le m_i\le\mathrm{INT64\_MAX}\)、剰余は符号付き 64 bit 整数の全範囲で指定できます。解がなければ std::nullopt を返します。解があれば、全ての整数解は

\[x=\mathrm{residue}+k\,\mathrm{modulus}\quad(k\in\mathbb{Z}),\qquad 0\le\mathrm{residue}<\mathrm{modulus}\]

です。空の連立式は {0, 1}、すなわち全ての整数を表します。同じ合同式を複数指定しても構いません。

配列長の不一致と正でない法は std::invalid_argument を送出します。合同式を入力順に併合する途中で周期が INT64_MAX を超えた場合は std::overflow_error を送出します。この場合、それ以降の合同式は処理しません。途中の積には 128 bit 整数を使います。

一次合同式 を使って 2 つずつ併合します。合同式の数を \(n\)、最大の法を \(M\) として、時間は \(O(n\log M)\)、追加空間は \(O(1)\) です。2 式版は \(O(\log\min(m_0,m_1))\) 時間です。法 1 の処理は \(O(1)\) 時間です。

依存ライブラリ

• 一次合同式

コード

#ifndef CRT_HPP
#define CRT_HPP

#include "linear_congruence.hpp"
#include <cstdint>
#include <limits>
#include <optional>
#include <stdexcept>
#include <utility>
#include <vector>

struct CRTSolution {
	std::int64_t residue, modulus;
	bool operator==(const CRTSolution &) const = default;
};

constexpr std::optional<CRTSolution> crt(std::int64_t r0, std::int64_t m0, std::int64_t r1, std::int64_t m1) {
	if (m0 <= 0 || m1 <= 0) throw std::invalid_argument("crt modulus");
	r0 %= m0;
	r1 %= m1;
	if (r0 < 0) r0 += m0;
	if (r1 < 0) r1 += m1;
	if (m0 < m1) {
		std::swap(r0, r1);
		std::swap(m0, m1);
	}
	if (m0 % m1 == 0) {
		if (r0 % m1 != r1) return std::nullopt;
		return CRTSolution{r0, m0};
	}
	auto solution = linear_congruence(m0, r1 - r0, m1);
	if (!solution) return std::nullopt;
	__int128_t modulus = __int128_t(m0) * solution->period;
	if (modulus > std::numeric_limits<std::int64_t>::max()) throw std::overflow_error("crt modulus");
	return CRTSolution{std::int64_t(r0 + __int128_t(m0) * solution->residue), std::int64_t(modulus)};
}

inline std::optional<CRTSolution> crt(const std::vector<std::int64_t> &residues, const std::vector<std::int64_t> &moduli) {
	if (residues.size() != moduli.size()) throw std::invalid_argument("crt size");
	for (auto m : moduli) if (m <= 0) throw std::invalid_argument("crt modulus");
	CRTSolution result{0, 1};
	for (std::size_t i = 0; i < residues.size(); ++i) {
		auto next = crt(result.residue, result.modulus, residues[i], moduli[i]);
		if (!next) return std::nullopt;
		result = *next;
	}
	return result;
}

#endif

Back to top page