洛谷题单指南-进阶数论-P1495 【模板】中国剩余定理（CRT）/ 曹冲养猪

原题链接：https://www.luogu.com.cn/problem/P1495

题意解读：求方程组x ≡ b_i (mod a_i), i∈[1,n]的最小正整数解，所有的a_i互质。

解题思路：

1、中国剩余定理

设方程组为（a1，a2，a3互质）：

• x ≡ b1 (mod a1)
• x ≡ b2 (mod a2)
• x ≡ b3 (mod a3)

要找到x满足上面三种情况，设a = a1 * a2 * a3

x需要拆解成三个数：x = (x1 + x2 + x3) % a，

必须满足：

• x1 % a1 = b1, x1 % a2 = 0, x1 % a3 = 0
• x2 % a2 = b2, x2 % a1 = 0, x2 % a3 = 0
• x3 % a3 = b3, x3 % a1 = 0, x3 % a2 = 0

可以看出，

x1是a2、a3的倍数，不妨设x1 = a2 * a3 * k1，a2 * a3 * k1 ≡ b1 (mod a1)，

转化为不定方程a2 * a3 * k1 + a1 * p1 = b1，

由于gcd(a2 * a3, a1) = 1，求a2 * a3 * k1 + a1 * p1 = 1的一个k1的解再乘以b1倍即可

又有a2 * a3 * k1 ≡ 1 (mod a1)，k1的解就是a2 * a3模a1的逆元，再乘以b1

因此：

k1 = (a2 * a3)^-1 * b1

x1 = a2 * a3 * (a2 * a3)^-1 * b1

同理分析，可得到：

x2 = a1 * a3 * (a1 * a3)^-1 * b2

x3 = a1 * a2 * (a1 * a2)^-1 * b3

x = a2 * a3 * (a2 * a3)^-1 * b1 + a1 * a3 * (a1 * a3)^-1 * b2 + a1 * a2 * (a1 * a2)^-1 * b3

= (a / a1) * (a / a1)^-1 * b1 + (a / a2) * (a / a2)^-1 * b2 + (a / a3) * (a / a3)^-1 * b3

x最后要模a

推而广之：

设a = a1 * a2 * ... * an

ti = a / ai，ti^-1是ti在模ai意义下的逆元

x = (∑ ti * ti^-1 * bi) % a

100分代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

const int N = 15;
LL a[N], b[N];
LL n, A = 1, ans;

LL exgcd(LL a, LL b, LL &x, LL &y)
{
    if(b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y = y - a / b * x;
    return d;
}

LL mul(LL a, LL b, LL mod)
{
    LL res = 0;
    while(b)
    {
        if(b & 1) res = (res + a) % mod;
        b >>= 1;
        a = (a + a) % mod;
    }
    return res;
}

int main()
{
    cin >> n;
    for(int i = 1; i <= n; i++)
    {
        cin >> a[i] >> b[i];
        A *= a[i];
    }
    for(int i = 1; i <= n; i++)
    {
        LL m = A / a[i];
        LL x, y;
        exgcd(m, a[i], x, y); //求m模a[i]的逆元
        x = (x % a[i] + a[i]) % a[i];
        //ans = (ans + m * x * b[i]) % A; //可能溢出
        ans = (ans + mul(m, x * b[i], A)) % A; //快速乘避免溢出
    }
    cout << ans;
    return 0;
}

2、扩展：同时适用于模数不互质的情况

设方程组：

• x ≡ b1 (mod a1)
• x ≡ b2 (mod a2)
• ...
• x ≡ bn (mod an)

先取前两个方程，展开

x = a1 * s + b1

x = a2 * t + b2

a1 * s - a2 * t = b2 - b1

调整负数为整数，因为t是任意值，所以系数正负没关系，正数便于计算

a1 * s + a2 * t = b2 - b1

用扩展欧几里得算法求a1 * s + a2 * t = gcd(a1, a2)的一个特解s0

如果(b2 - b1) % d != 0，则方程组无解！

设d = gcd(a1, a2), tmp1 = (b2 - b1) / d，tmp2 = a2 / d

s的通解为：s = s0 * tmp1 + tmp2 * k，k为任意整数 

注意：在计算s0 * tmp1时如果数据范围过大可考虑模tmp2的快速乘！前提是tmp1是非负数，即tmp1 = (tmp1 % tmp2 + tmp2) % tmp2

将s的通解代入x = a1 * s + b1得到x = a1 * (s0 * tmp1 + tmp2 * k) + b1，展开得到

x = a1 * tmp2 * k + (a1 * s0 * tmp1 + b1) 

将方程x = a1 * tmp2 * k + (a1 * s0 * tmp1 + b1) 看作x = a1 * s + b1，可以将a1 = a1 * tmp2, b1 = a1 * s0 * tmp1 + b1

将第三个方程x = a3 * t + b3看作x = a2 * t + b2

继续上述合并方程的过程，最后处理完所有的方程，最终b1就是方程合并之后的一个解，对(b1 % a1 + a1) % a1即得到最小正整数解。

100分代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;

const int N = 15;
LL a[N], b[N];
LL n, A = 1, ans;

LL exgcd(LL a, LL b, LL &x, LL &y)
{
    if(b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y = y - a / b * x;
    return d;
}

LL mul(LL a, LL b, LL mod)
{
    LL res = 0;
    while(b)
    {
        if(b & 1) res = (res + a) % mod;
        b >>= 1;
        a = (a + a) % mod;
    }
    return res;
}

int main()
{
    cin >> n;
    for(int i = 1; i <= n; i++)
    {
        cin >> a[i] >> b[i];
    }
    LL a1 = a[1], b1 = b[1];
    for(int i = 2; i <= n; i++)
    {
        int a2 = a[i], b2 = b[i];
        LL s, t;
        LL d = exgcd(a1, a2, s, t);
        if((b2 - b1) % d)
        {
            //无解的情况此题不用考虑
        }
        LL tmp1 = (b2 - b1) / d, tmp2 = a2 / d;
        tmp1 = (tmp1 % tmp2 + tmp2) % tmp2;
        s = mul(s, tmp1, tmp2);
        b1 = a1 * s + b1;
        a1 = a1 * tmp2;
    }
    cout << (b1 % a1 + a1) % a1;
    return 0;
}