Garner (連立一次合同式 mod M)
(Mathematics/garner.hpp)

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

概要

連立一次合同式の解 $\pmod M$ を求めます。
似たような問題として、中国剰余定理があります。

関数

bool T pre_garner(const vector<T> & b, const vector<T> & m)

連立一次合同式

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

と同値な連立一次合同式であって、$m_i, m_j$ が互いに素、$0 \leq b_i < m_i$ となるように書き換えます。
この合同式に解が存在するならば $true$ を、存在しないのならば $false$ を返します。
連立一次合同式の解の存在条件や一意性については中国剰余定理を参照してください。

制約

T は int, long long
$0 < m_i$

計算量

$\mathcal{O}(K^2 \log{\max{m_i}} \log{\log{\max{m_i}}}))$

template<typename T, typename U> T garner(const vector<T> & b, const vector<T> & m, const T M)

$m_i, m_j (i \neq j)$ が互いに素な連立一次合同式

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

の解 $x \pmod lcm(m_1, m_2, \ldots, m_K)$ に対して、 $x \pmod M$ を返します。
内部の計算で使用する型を U に指定してください。

制約

T, U は int または long long
$\max(m_i, M)^2$ が U で表現可能
$0 < m_i, M$

計算量

$\mathcal{O}(K(K + \log{\max{m_i}}))$

使用例

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

int main() {
	// 17
	cout << "garner({3, 6}, {7, 11}, 10000) = " << tk::garner<int, int>({3, 6}, {7, 11}, 10000) << endl;
	// 2
	cout << "garner({3, 6}, {7, 11}, 5) = " << tk::garner<int, int>({3, 6}, {7, 11}, 5) << endl;
	
	vector<int> b{2, 3, 4}, m{10, 20, 30};
	if (!tk::pre_garner(b, m)) cout << "no" << endl; // no
	
	b = {2, 12, 22};
	if (tk::pre_garner(b, m)) {
		for (int i = 0; i < 3; ++i) cout << "b[i] = " << b[i] << ", m[i] = " << m[i] << endl;
		// 52
		cout << "garner(b, m, 1000000) = " << tk::garner<int, long long>(b, m, 1000000) << endl;
		// 2
		cout << "garner(b, m, 10) = " << tk::garner<int, int>(b, m, 10) << endl;
	}
}

参考

2020/05/05: https://qiita.com/drken/items/ae02240cd1f8edfc86fd

Depends on

Required by

任意 mod 数論変換 (任意 mod NTT) (Convolution/NumberTheoreticTransform_AnyMod.hpp)

Verified with

Code

#ifndef INCLUDE_GUARD_GARNER_HPP
#define INCLUDE_GUARD_GARNER_HPP

#include <vector>
#include <cassert>
#include <type_traits>

#include "Mathematics/euclid.hpp"
#include "Mathematics/mod_pow_inv.hpp"

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/garner.hpp
 */
namespace tk {
template<typename T>
bool pre_garner(std::vector<T> & b, std::vector<T> & m) noexcept {
	static_assert(std::is_integral<T>::value);
	static_assert(std::is_signed<T>::value);
	for (int i = 0; i < static_cast<int>(b.size()); ++i) {
		b[i] = b[i] % m[i] + (b[i] >= 0 ? 0 : m[i]);
		for (int j = 0; j < i; ++j) {
			T g = gcd(m[i], m[j]);
			if ((b[i] - b[j]) % g != 0) return false;
			m[i] /= g; m[j] /= g;
			T gi = gcd(g, m[i]), gj = g / gi;
			do {
				g = gcd(gi, gj);
				gi *= g; gj /= g;
			} while (g != 1);
			m[i] *= gi; m[j] *= gj;
			b[i] %= m[i]; b[j] %= m[j];
		}
	}
	return true;
}

template<typename T, typename U>
T garner(const std::vector<T> & b, const std::vector<T> & m, const T M) {
	static_assert(std::is_integral<T>::value);
	assert(b.size() == m.size());
	const int n = b.size();
	assert(n > 0);
	{
		T g = 0;
		for (auto v : m) {
			assert(v > 0);
			g = gcd(g, v);
		}
		assert(n == 1 || g == 1);
	}
	assert(M > 0);
	std::vector<T> sum(n + 1, 0), ip(n + 1, 1);
	for (int i = 0; i < n; ++i) {
		if (m[i] == 1) continue;
		U t = (b[i] % m[i] + static_cast<U>(2) * m[i] - sum[i]) % m[i] * mod_inv(ip[i], m[i]) % m[i];
		for (int j = i + 1; j < n; ++j) {
			sum[j] = (sum[j] + ip[j] * t) % m[j];
			ip[j] = static_cast<U>(ip[j]) * m[i] % m[j];
		}
		sum[n] = (sum[n] + ip[n] * t % M) % M;
		ip[n] = static_cast<U>(ip[n]) * m[i] % M; 
	}
	return sum.back();
}
} // namespace tk


#endif // INCLUDE_GUARD_GARNER_HPP

#line 1 "Mathematics/garner.hpp"



#include <vector>
#include <cassert>
#include <type_traits>

#line 1 "Mathematics/euclid.hpp"



#line 5 "Mathematics/euclid.hpp"
#include <utility>
#include <tuple>
#line 8 "Mathematics/euclid.hpp"
#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 1 "Mathematics/mod_pow_inv.hpp"



#line 6 "Mathematics/mod_pow_inv.hpp"

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/mod_pow_inv.hpp
 */
namespace tk {
template<typename T>
constexpr T mod_pow(T x, T n, T m) noexcept {
	static_assert(std::is_integral<T>::value);
	assert(m > 0);
	assert(n >= 0);
	x = x % m + (x >= 0 ? 0 : m);
	T res = 1 % m;
	while (n > 0) {
		if (n & 1) res = res * x % m;
		x = x * x % m;
		n >>= 1;
	}
	return res;
}

template<typename T>
constexpr T mod_inv(T x, T m) noexcept {
	static_assert(std::is_integral<T>::value);
	static_assert(std::is_signed<T>::value);
	assert(m > 0);
	x = x % m + (x >= 0 ? 0 : m);
	T x1 = 1, y = m, y1 = 0;
	while (y > 0) {
		const T q = x / y;
		T tmp = x - q * y; x = y; y = tmp;
		tmp = x1 - q * y1; x1 = y1; y1 = tmp;
	}
	assert(x == 1);
	if (x1 == m) x1 = 0;
	if (x1 < 0) x1 += m;
	return x1;
}
} // namespace tk



#line 10 "Mathematics/garner.hpp"

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/garner.hpp
 */
namespace tk {
template<typename T>
bool pre_garner(std::vector<T> & b, std::vector<T> & m) noexcept {
	static_assert(std::is_integral<T>::value);
	static_assert(std::is_signed<T>::value);
	for (int i = 0; i < static_cast<int>(b.size()); ++i) {
		b[i] = b[i] % m[i] + (b[i] >= 0 ? 0 : m[i]);
		for (int j = 0; j < i; ++j) {
			T g = gcd(m[i], m[j]);
			if ((b[i] - b[j]) % g != 0) return false;
			m[i] /= g; m[j] /= g;
			T gi = gcd(g, m[i]), gj = g / gi;
			do {
				g = gcd(gi, gj);
				gi *= g; gj /= g;
			} while (g != 1);
			m[i] *= gi; m[j] *= gj;
			b[i] %= m[i]; b[j] %= m[j];
		}
	}
	return true;
}

template<typename T, typename U>
T garner(const std::vector<T> & b, const std::vector<T> & m, const T M) {
	static_assert(std::is_integral<T>::value);
	assert(b.size() == m.size());
	const int n = b.size();
	assert(n > 0);
	{
		T g = 0;
		for (auto v : m) {
			assert(v > 0);
			g = gcd(g, v);
		}
		assert(n == 1 || g == 1);
	}
	assert(M > 0);
	std::vector<T> sum(n + 1, 0), ip(n + 1, 1);
	for (int i = 0; i < n; ++i) {
		if (m[i] == 1) continue;
		U t = (b[i] % m[i] + static_cast<U>(2) * m[i] - sum[i]) % m[i] * mod_inv(ip[i], m[i]) % m[i];
		for (int j = i + 1; j < n; ++j) {
			sum[j] = (sum[j] + ip[j] * t) % m[j];
			ip[j] = static_cast<U>(ip[j]) * m[i] % m[j];
		}
		sum[n] = (sum[n] + ip[n] * t % M) % M;
		ip[n] = static_cast<U>(ip[n]) * m[i] % M; 
	}
	return sum.back();
}
} // namespace tk

Garner (連立一次合同式 mod M) (Mathematics/garner.hpp)

概要

関数