中国剩余定理

中国剩余定理

2023-05-27

CRT / exCRT

upd on 04/27/25：发现该文没写完

给定下列关于 \(x\) 的一元同余方程组：

\[ \begin {cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \\ \quad \quad \vdots \\ x \equiv a_k \pmod {m_k} \end {cases} \] 其中对于 \(\forall \, 1\le i, j\le k\, (i\ne j)\)，满足 \(m_i \, \bot\, m_j\)。下面是求解该方程组的具体方法：

令 \(M = \prod \limits _{i = 1} ^ k m_i, \, M_i=\dfrac M {m_i}\)，则 \(x = \sum \limits _{i = 1} ^k a_i\times M_i\times M_i ^{-1} \pmod M\)。

其中，\(M_i^{-1}\) 是 \(M_i\) 在模 \(m_i\) 意义下的逆元（所以 \(M_i\times M_i^{-1}\) 的值并不是视觉上的 \(1\)）。

正确性证明：对于 \(\forall \, i\in \{1, 2, \cdots, k\}\)，有 \(M_i\mid M\)；对于 \(\forall \, j\in \{1, 2, \cdots, k\}\) 且 \(j\ne i\)，有 \(m_i\mid \dfrac M{m_j}\) 即 \(m_i\mid M_j\)，即 \(M_j\equiv 0\pmod {m_i}\)。

又 \(m_i\mid M\) 那么就有 \(x\equiv a_i\cdot M_i\cdot {M_i}^-1\equiv a_i\pmod {n_i}\)。

时间复杂度 \(\mathcal O(n\log n)\)。注意求解逆元时要用到 exgcd。

inline int CRT(int k, int *a, int *m) {
    int M = 1, x = 0;
    static int M1[maxn], Mi[maxn]; 
    for (int i = 1; i <= k; ++i)
        M *= m[i];
    for (int i = 1; i <= k; ++i) {
        M1[i] = M / m[i];
        Mi[i] = getinv(M1[i], m[i]);
        (x += a[i] * M1[i] * Mi[i]) %= M;
    }
    return x;
}

A. 曹冲养猪

http://222.180.160.110:1024/contest/3642/problem/1

板。

#define int long long
namespace XSC062 {
using namespace fastIO;
const int maxn = 15;
int n;
int a[maxn], b[maxn];
int exgcd(int a, int &x, int b, int &y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int u = exgcd(b, x, a % b, y);
    int t = x;
    x = y, y = t - (a / b) * y;
    return u;
}
inline int getinv(int a, int p) {
    int k1, k2;
    int g = exgcd(a, k1, p, k2);
    if (g != 1)
        return -1;
    return (k1 % p + p) % p;
}
inline int CRT(int k, int *a, int *m) {
    int M = 1, x = 0;
    static int M1[maxn], Mi[maxn]; 
    for (int i = 1; i <= k; ++i)
        M *= m[i];
    for (int i = 1; i <= k; ++i) {
        M1[i] = M / m[i];
        Mi[i] = getinv(M1[i], m[i]);
        (x += a[i] * M1[i] * Mi[i]) %= M;
    }
    return x;
}
int main() {
    read(n);
    for (int i = 1; i <= n; ++i)
        read(a[i]), read(b[i]);
    print(CRT(n, b, a));
    return 0;
}
} // namespace XSC062
#undef int

B. 阿九大战朱最学

http://222.180.160.110:1024/contest/3642/problem/2

它甚至连样例都不愿意改一下，就是上一道题啊。

C. 韩信点兵

http://222.180.160.110:1024/contest/3642/problem/3

还是板。

#define int long long
namespace XSC062 {
using namespace fastIO;
const int maxn = 15;
int n, res;
int a[maxn], b[maxn];
int exgcd(int a, int &x, int b, int &y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int u = exgcd(b, x, a % b, y);
    int t = x;
    x = y, y = t - (a / b) * y;
    return u;
}
inline int getinv(int a, int p) {
    int k1, k2;
    int g = exgcd(a, k1, p, k2);
    if (g != 1)
        return -1;
    return (k1 % p + p) % p;
}
inline int CRT(int k, int *a, int *m) {
    int M = 1, x = 0;
    static int M1[maxn], Mi[maxn]; 
    for (int i = 1; i <= k; ++i)
        M *= m[i];
    for (int i = 1; i <= k; ++i) {
        M1[i] = M / m[i];
        Mi[i] = getinv(M1[i], m[i]);
        (x += a[i] * M1[i] * Mi[i]) %= M;
    }
    return x;
}
int main() {
    n = 3;
    a[1] = 3, a[2] = 5, a[3] = 7;
    read(b[1]), read(b[2]), read(b[3]);
    res = CRT(n, b, a);
    if (res < 10 || res > 100)
        puts("no answer");
    else print(res);
    return 0;
}
} // namespace XSC062
#undef int

扩展中国剩余定理

设有如下同余方程组：

\[ \begin {cases} x \equiv a_1 \pmod {m_1} \\ x \equiv a_2 \pmod {m_2} \end {cases} \]

不保证 \(m_1 \, \bot \, m_2\)。

设 \(x = m_1 \times p_1 + a_1 = m_2 \times p_2 + a_2\)，其中 \(p_i \in \mathbb Z\)，则有 \(p_1 \times m_1 - p_2 \times m_2 = a_2 - a_1\)。

\(p_1\) 的值可以用 exgcd 求解，则原方程的解满足 \(x \equiv m_1\times p_1 + a_1 \pmod {M}\)，其中 \(M=\operatorname{lcm}(m_1, m_2)\)。

这样我们就得到了一个新的同余方程。对于 \(k>2\) 的情况，我们不断合并两个同余方程即可得到最终同余方程。

此时根据题目要求（一般求最小值即 \(m_k\times p_k + a_k\)）求解答案即可。

时间复杂度 \(\mathcal O(n\log n)\)，其中 \(\log\) 来自 exgcd。

inline int calc(int m1, int a1,
                    int m2, int a2) {
    int p1, p2;
    if (a2 - a1 < 0)
        swap(a1, a2), swap(m1, m2);
    int g = exgcd(m1, p1, m2, p2);
    if ((a2 - a1) % g)
        return -1;
    p1 *= (a2 - a1) / g, m2 /= g;
    p1 = (p1 % m2 + m2) % m2;
    return p1 * m1 + a1;
}
inline int exCRT(int k, int *a, int *m) {
    for (int i = 2; i <= k; ++i) {
        a[i] = calc(m[i - 1],
                a[i - 1], m[i], a[i]);
        if (a[i] == -1)
            return -1;
        m[i] = lcm(m[i - 1], m[i]);
    }
    return a[k] % m[k];
}

D. Strange Way to Express Integers

http://222.180.160.110:1024/contest/3642/problem/4

板。但是要开 __int128。不知道为什么智力只用开 long long 就能跑过，我和揭哥就不行。

#define int __int128
namespace XSC062 {
using namespace fastIO;
const int maxn = 1e5 + 5;
int n;
int a[maxn], m[maxn];
int gcd(int x, int y) {
    return y ? gcd(y, x % y) : x;
}
inline int lcm(int x, int y) {
    return x / gcd(x, y) * y;
}
inline void swap(int &x, int &y) {
    x ^= y ^= x ^= y;
    return;
}
int exgcd(int a, int &x, int b, int &y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int u = exgcd(b, x, a % b, y);
    int t = x;
    x = y, y = t - (a / b) * y;
    return u;
}
inline int calc(int m1, int a1,
                    int m2, int a2) {
    int p1, p2;
    if (a2 - a1 < 0)
        swap(a1, a2), swap(m1, m2);
    int g = exgcd(m1, p1, m2, p2);
    if ((a2 - a1) % g)
        return -1;
    p1 *= (a2 - a1) / g, m2 /= g;
    p1 = (p1 % m2 + m2) % m2;
    return p1 * m1 + a1;
}
inline int exCRT(int k, int *a, int *m) {
    for (int i = 2; i <= k; ++i) {
        a[i] = calc(m[i - 1],
                a[i - 1], m[i], a[i]);
        if (a[i] == -1)
            return -1;
        m[i] = lcm(m[i - 1], m[i]);
    }
    return a[k] % m[k];
}
int main() {
    while (read(n)) {
        for (int i = 1; i <= n; ++i)
            read(m[i]), read(a[i]);
        print(exCRT(n, a, m), '\n');
    }
    return 0;
}
} // namespace XSC062
#undef int

言论

学习小组

二分图