中国剰余定理 (Chinese remainder theorem)
(include/emthrm/math/chinese_remainder_theorem.hpp)

• View this file on GitHub
• Last update: 2023-02-23 21:59:12+09:00
• Include: #include "emthrm/math/chinese_remainder_theorem.hpp"

連立合同式 $x \equiv b_i \pmod{m_i}$ ($i = 1, 2,\ldots, n$) が与えられる。

任意の $i, j \in \lbrace 1, 2, \ldots, n \rbrace$ に対して $m_i \perp m_j$ が成り立つならば、連立合同式を満たす $x \bmod{\prod_{i = 1}^n m_i}$ が一意に存在する。

これを任意の $i, j \in \lbrace 1, 2, \ldots, n \rbrace$ に対して $b_i \equiv b_j \pmod{\gcd(m_i, m_j)}$ が成り立つときに拡張すると、連立合同式を満たす $x \bmod{\mathrm{lcm}(m_1, m_2,…, m_n)}$ が一意に存在する。

時間計算量

$O(N \log{\mathrm{lcm}(m_1, m_2, \ldots, m_N)})$

仕様

名前	戻り値	要件
template <typename T> std::pair<T, T> chinese_remainder_theorem(std::vector<T> b, std::vector<T> m);	$x \equiv b_i \pmod{m_i}$ を満たす $x$ と $\mathrm{lcm}(m_1, m_2,…, m_n)$ の組。ただし存在しないときは $(0, 0)$ を返す。	$0 \leq x < \mathrm{lcm}(m_1, m_2,…, m_n)$

参考文献

• https://rsk0315.hatenablog.com/entry/2021/01/18/065720#crt
• https://github.com/atcoder/ac-library/blob/master/atcoder/math.hpp#L33-L80
• https://qiita.com/drken/items/ae02240cd1f8edfc86fd
• https://mathtrain.jp/remainder

TODO

• Garner のアルゴリズム
  • https://qiita.com/drken/items/ae02240cd1f8edfc86fd
  • https://github.com/drken1215/algorithm/blob/master/MathNumberTheory/garner.cpp
  • https://www.creativ.xyz/ect-gcd-crt-garner-927/
  • http://kirika-comp.hatenablog.com/entry/2017/12/18/143923
  • http://sugarknri.hatenablog.com/entry/2018/06/28/144407
  • https://lumakernel.github.io/ecasdqina/math/number-theory/Garner
  • https://mugen1337.github.io/procon/Math/Garner.cpp
  • https://sotanishy.github.io/cp-library-cpp/math/garner.cpp
  • https://cp-algorithms.com/algebra/chinese-remainder-theorem.html
  • http://flex.phys.tohoku.ac.jp/~maru/implementations/garner.php

Submissons

https://yukicoder.me/submissions/630411

Depends on

• 逆元 (inverse element) (include/emthrm/math/mod_inv.hpp)

Verified with

• 数学/中国剰余定理 (test/math/chinese_remainder_theorem.test.cpp)

Code

#ifndef EMTHRM_MATH_CHINESE_REMAINDER_THEOREM_HPP_
#define EMTHRM_MATH_CHINESE_REMAINDER_THEOREM_HPP_

#include <numeric>
#include <utility>
#include <vector>

#include "emthrm/math/mod_inv.hpp"

namespace emthrm {

template <typename T>
std::pair<T, T> chinese_remainder_theorem(std::vector<T> b, std::vector<T> m) {
  const int n = b.size();
  T x = 0, md = 1;
  for (int i = 0; i < n; ++i) {
    if ((b[i] %= m[i]) < 0) b[i] += m[i];
    if (md < m[i]) {
      std::swap(x, b[i]);
      std::swap(md, m[i]);
    }
    if (md % m[i] == 0) {
      if (x % m[i] != b[i]) return {0, 0};
      continue;
    }
    const T g = std::gcd(md, m[i]);
    if ((b[i] - x) % g != 0) return {0, 0};
    const T u_i = m[i] / g;
    x += (b[i] - x) / g % u_i * mod_inv(md / g, u_i) % u_i * md;
    md *= u_i;
    if (x < 0) x += md;
  }
  return {x, md};
}

}  // namespace emthrm

#endif  // EMTHRM_MATH_CHINESE_REMAINDER_THEOREM_HPP_

#line 1 "include/emthrm/math/chinese_remainder_theorem.hpp"



#include <numeric>
#include <utility>
#include <vector>

#line 1 "include/emthrm/math/mod_inv.hpp"



#line 6 "include/emthrm/math/mod_inv.hpp"

namespace emthrm {

long long mod_inv(long long a, const int m) {
  if ((a %= m) < 0) a += m;
  if (std::gcd(a, m) != 1) return -1;
  long long x = 1;
  for (long long b = m, u = 0; b > 0;) {
    const long long q = a / b;
    std::swap(a -= q * b, b);
    std::swap(x -= q * u, u);
  }
  x %= m;
  return x < 0 ? x + m : x;
}

}  // namespace emthrm


#line 9 "include/emthrm/math/chinese_remainder_theorem.hpp"

namespace emthrm {

template <typename T>
std::pair<T, T> chinese_remainder_theorem(std::vector<T> b, std::vector<T> m) {
  const int n = b.size();
  T x = 0, md = 1;
  for (int i = 0; i < n; ++i) {
    if ((b[i] %= m[i]) < 0) b[i] += m[i];
    if (md < m[i]) {
      std::swap(x, b[i]);
      std::swap(md, m[i]);
    }
    if (md % m[i] == 0) {
      if (x % m[i] != b[i]) return {0, 0};
      continue;
    }
    const T g = std::gcd(md, m[i]);
    if ((b[i] - x) % g != 0) return {0, 0};
    const T u_i = m[i] / g;
    x += (b[i] - x) / g % u_i * mod_inv(md / g, u_i) % u_i * md;
    md *= u_i;
    if (x < 0) x += md;
  }
  return {x, md};
}

}  // namespace emthrm

Back to top page