中国剰余定理
(Mathematics/chinese_remainder.hpp)

• View this file on GitHub
• Last update: 2021-06-05 13:56:14+09:00
• Include: #include "Mathematics/chinese_remainder.hpp"
• Link: https://tkmst201.github.io/Library/Mathematics/chinese_remainder.hpp

概要

中国剰余定理です。連立一次合同式を解きます。
解 $x \pmod M$ を求めたい場合は Garner を使用してください。

関数

pair<T, T> chinese_remainder(T b1, T m1, T b2, T m2)

連立一次合同式

\[x \equiv b_1 \pmod{m_1}\] \[x \equiv b_2 \pmod{m_2}\]

の解の一つを $0 \leq x < lcm(m_1, m_2)$ の範囲で求め、$(x, lcm(m_1, m_2))$ を返します。
解が存在しない場合は $(0, 0)$ を返します。
また、この合同式の解は、$b_1 \equiv b_2 \pmod{gcd(m_1, m_2)}$ であるときに限り存在し、周期 $lcm(m_1, m_2)$ の中で一意です。

制約

• T は int, long long
• $lcm(m_1, m_2)$ が T で表現可能
• $0 < m_1, m_2$

計算量

• $\mathcal{O}(\log{\min(m1, m2)})$

---

使用例

#include <bits/stdc++.h>
#include "Mathematics/chinese_remainder.hpp"
using namespace std;

int main() {
	auto [a, l] = tk::chinese_remainder(2, 3, 5, 12);
	cout << "chinese_remainder(2, 3, 5, 12) = (" << a << ", " << l << ")" << endl; // (5, 12)
	
	tie(a, l) = tk::chinese_remainder(3, 7, 6, 11);
	cout << "chinese_remainder(3, 7, 6, 11) = (" << a << ", " << l << ")" << endl; // (17, 77)
	
	tie(a, l) = tk::chinese_remainder(0, 10, 0, 4);
	cout << "chinese_remainder(0, 10, 0, 4) = (" << a << ", " << l << ")" << endl; // (0, 20)
	
	tie(a, l) = tk::chinese_remainder(2, 10, 1, 8);
	cout << "chinese_remainder(2, 10, 1, 8) = (" << a << ", " << l << ")" << endl; // (0, 0)
}

参考

2021/03/18: https://qiita.com/drken/items/ae02240cd1f8edfc86fd
2020/09/21: AC Library

Depends on

• ユークリッドの互除法 (Mathematics/euclid.hpp)

Verified with

• Test/chinese_remainder.test.cpp

Code

#ifndef INCLUDE_GUARD_CHINESE_REMAINDER_HPP
#define INCLUDE_GUARD_CHINESE_REMAINDER_HPP

#include "Mathematics/euclid.hpp"

#include <cassert>
#include <utility>
#include <type_traits>

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/chinese_remainder.hpp
 */
namespace tk {
template<typename T>
constexpr std::pair<T, T> chinese_remainder(T b1, T m1, T b2, T m2) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(m1 > 0);
	assert(m2 > 0);
	if (m1 < m2) { std::swap(b1, b2); std::swap(m1, m2); }
	b1 = b1 % m1 + (b1 >= 0 ? 0 : m1);
	b2 = b2 % m2 + (b2 >= 0 ? 0 : m2);
	auto [g, x, _] = ext_gcd(m1, m2);
	if ((b2 - b1) % g != 0) return {0, 0};
	const T pm2 = m2 / g;
	if (x < 0) x += pm2;
	const T t = ((b2 - b1) / g % pm2 + pm2) % pm2 * x % pm2;
	return {b1 + t * m1, m1 * pm2};
}
} // namespace tk


#endif // INCLUDE_GUARD_CHINESE_REMAINDER_HPP

#line 1 "Mathematics/chinese_remainder.hpp"



#line 1 "Mathematics/euclid.hpp"



#include <cassert>
#include <utility>
#include <tuple>
#include <type_traits>
#include <cmath>

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/euclid.hpp
 */
namespace tk {
template<typename T>
constexpr T gcd(T a, T b) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(a >= 0);
	assert(b >= 0);
	while (b != 0) {
		const T t = a % b;
		a = b; b = t;
	}
	return a;
}

template<typename T>
constexpr T lcm(T a, T b) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(a >= 0);
	assert(b >= 0);
	if (a == 0 || b == 0) return 0;
	return a / gcd(a, b) * b;
}

template<typename T>
constexpr std::tuple<T, T, T> ext_gcd(T a, T b) noexcept {
	static_assert(std::is_integral<T>::value);
	static_assert(std::is_signed<T>::value);
	assert(a != 0);
	assert(b != 0);
	T a1 = (a > 0) * 2 - 1, a2 = 0, b1 = 0, b2 = (b > 0) * 2 - 1;
	a = std::abs(a);
	b = std::abs(b);
	while (b > 0) {
		const T q = a / b;
		T tmp = a - q * b; a = b; b = tmp;
		tmp = a1 - q * b1; a1 = b1; b1 = tmp;
		tmp = a2 - q * b2; a2 = b2; b2 = tmp;
	}
	return {a, a1, a2};
}
} // namespace tk



#line 5 "Mathematics/chinese_remainder.hpp"

#line 9 "Mathematics/chinese_remainder.hpp"

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/chinese_remainder.hpp
 */
namespace tk {
template<typename T>
constexpr std::pair<T, T> chinese_remainder(T b1, T m1, T b2, T m2) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(m1 > 0);
	assert(m2 > 0);
	if (m1 < m2) { std::swap(b1, b2); std::swap(m1, m2); }
	b1 = b1 % m1 + (b1 >= 0 ? 0 : m1);
	b2 = b2 % m2 + (b2 >= 0 ? 0 : m2);
	auto [g, x, _] = ext_gcd(m1, m2);
	if ((b2 - b1) % g != 0) return {0, 0};
	const T pm2 = m2 / g;
	if (x < 0) x += pm2;
	const T t = ((b2 - b1) / g % pm2 + pm2) % pm2 * x % pm2;
	return {b1 + t * m1, m1 * pm2};
}
} // namespace tk

Back to top page