中国剰余定理
(math/crt.hpp)

• View this file on GitHub
• Last update: 2024-06-26 03:18:46+09:00
• Include: #include "math/crt.hpp"

概要

中国剰余定理

x ≡ b1 mod m1, x ≡ b2 mod m2をみたすxがx ≡ r mod mと書けるとき，

(r, m)を返す，解がなければ(0, -1)を返す

使い方

x ≡ bi (mod mi) を満たすlong long型のvector b, mを入力．

long long型のr, mを返す．

制約

x ≡ b1 mod m1, x ≡ b2 mod m2をみたすx ≡ r mod mが存在するには，

b1 ≡ b2 mod gcd(m1, m2)が必要

計算

d = gcd(m1, m2)とおくと，m1p + m2q = dを満たす(p, q)が構成できる

b1 ≡ b2 mod dより，b2 - b1はdで割り切れる．

s = (b2 - b1) / dとおくと，sm1p + sm2q = b2 - b1となる．

式変形し，b1 + sm1p = b2 - sm2qとし，左辺をxとおくと，

x ≡ b1 mod m1, x ≡ b2 mod m2を満たす．

Depends on

• 拡張ユークリッドの互除法 (math/ext_gcd.hpp)

Verified with

• verify/yukicoder/yukicoder_186.test.cpp

Code

#include "ext_gcd.hpp"

inline long long safe_mod(long long a, long long m) {
    return (a % m + m) % m;
}

pair<long long, long long> crt(const vector<long long>& b, const vector<long long>& m) {
    assert(b.size() == m.size());
    long long r = 0, M = 1;
    for(int i = 0; i < (int)b.size(); i++) {
        auto [d, p, q] = ext_gcd(M, m[i]);
        if((b[i] - r) % d != 0) return {0, -1};
        long long tmp = (b[i] - r) / d * p % (m[i] / d);
        r += M * tmp;
        M *= m[i] / d;
    }
    r = safe_mod(r, M);
    return {r, M};
}

#line 1 "math/ext_gcd.hpp"
tuple<long long, long long, long long> ext_gcd(long long a, long long b) {
    if(b == 0) return {a, 1, 0};
    auto [d, y, x] = ext_gcd(b, a % b);
    y -= a / b * x;
    return {d, x, y};
}
#line 2 "math/crt.hpp"

inline long long safe_mod(long long a, long long m) {
    return (a % m + m) % m;
}

pair<long long, long long> crt(const vector<long long>& b, const vector<long long>& m) {
    assert(b.size() == m.size());
    long long r = 0, M = 1;
    for(int i = 0; i < (int)b.size(); i++) {
        auto [d, p, q] = ext_gcd(M, m[i]);
        if((b[i] - r) % d != 0) return {0, -1};
        long long tmp = (b[i] - r) / d * p % (m[i] / d);
        r += M * tmp;
        M *= m[i] / d;
    }
    r = safe_mod(r, M);
    return {r, M};
}

Back to top page