原题链接:https://www.luogu.com.cn/problem/P1495
题意解读:求方程组x ≡ bi (mod ai), i∈[1,n]的最小正整数解,所有的ai互质。
解题思路:
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 * s + b3看作x = a2 * s + b2
继续上述合并方程求s的过程,最后处理完所有的方程,对(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;
}