Chinese Remainder Theorem (中國剩餘定理)
(library/math/crt.hpp)

• View this file on GitHub
• Last update: 2023-04-21 21:20:30+08:00
• Include: #include "library/math/crt.hpp"

中國剩餘定理

找 $x$ 使得 $x \equiv r[i] \pmod{m[i]}$。如果無解的話回傳 $(0, 0)$。

時間複雜度：$O(n \log{\mathrm{lcm}(m[i])})$

References

AtCoder Library document

Depends on

• library/math/inv-gcd.hpp
• library/math/safe-mod.hpp

Code

#pragma once
#include <vector>

#include <tuple>

#include <cassert>

#include <algorithm>

#include "safe-mod.hpp"

#include "inv-gcd.hpp"


namespace felix {

// (rem, mod)

template<class T>
std::pair<T, T> crt(const std::vector<T>& r, const std::vector<T>& m) {
	assert(r.size() == m.size());
	int n = (int) r.size();
	// Contracts: 0 <= r0 < m0

	T r0 = 0, m0 = 1;
	for(int i = 0; i < n; i++) {
		assert(1 <= m[i]);
		T r1 = internal::safe_mod(r[i], m[i]);
		T m1 = m[i];
		if(m0 < m1) {
			std::swap(r0, r1);
			std::swap(m0, m1);
		}
		if(m0 % m1 == 0) {
			if(r0 % m1 != r1) {
				return {0, 0};
			}
			continue;
		}
		T g, im;
		std::tie(g, im) = internal::inv_gcd(m0, m1);
		T u1 = (m1 / g);
		if((r1 - r0) % g) {
			return {0, 0};
		}
		T x = (r1 - r0) / g % u1 * im % u1;
		r0 += x * m0;
		m0 *= u1;
		if(r0 < 0) {
			r0 += m0;
		}
	}
	return {r0, m0};
}

} // namespace felix

#line 2 "library/math/crt.hpp"
#include <vector>

#include <tuple>

#include <cassert>

#include <algorithm>

#line 2 "library/math/safe-mod.hpp"

namespace felix {

namespace internal {

template<class T>
constexpr T safe_mod(T x, T m) {
	x %= m;
	if(x < 0) {
		x += m;
	}
	return x;
}

} // namespace internal


} // namespace felix
#line 3 "library/math/inv-gcd.hpp"

namespace felix {

namespace internal {

template<class T>
constexpr std::pair<T, T> inv_gcd(T a, T b) {
	a = safe_mod(a, b);
	if(a == 0) {
		return {b, 0};
	}
	T s = b, t = a;
	T m0 = 0, m1 = 1;
	while(t) {
		T u = s / t;
		s -= t * u;
		m0 -= m1 * u;
		auto tmp = s;
		s = t;
		t = tmp;
		tmp = m0;
		m0 = m1;
		m1 = tmp;
	}
	if(m0 < 0) {
		m0 += b / s;
	}
	return {s, m0};
}

} // namespace internal


} // namespace felix

#line 8 "library/math/crt.hpp"

namespace felix {

// (rem, mod)

template<class T>
std::pair<T, T> crt(const std::vector<T>& r, const std::vector<T>& m) {
	assert(r.size() == m.size());
	int n = (int) r.size();
	// Contracts: 0 <= r0 < m0

	T r0 = 0, m0 = 1;
	for(int i = 0; i < n; i++) {
		assert(1 <= m[i]);
		T r1 = internal::safe_mod(r[i], m[i]);
		T m1 = m[i];
		if(m0 < m1) {
			std::swap(r0, r1);
			std::swap(m0, m1);
		}
		if(m0 % m1 == 0) {
			if(r0 % m1 != r1) {
				return {0, 0};
			}
			continue;
		}
		T g, im;
		std::tie(g, im) = internal::inv_gcd(m0, m1);
		T u1 = (m1 / g);
		if((r1 - r0) % g) {
			return {0, 0};
		}
		T x = (r1 - r0) / g % u1 * im % u1;
		r0 += x * m0;
		m0 *= u1;
		if(r0 < 0) {
			r0 += m0;
		}
	}
	return {r0, m0};
}

} // namespace felix

Back to top page