一、一元线性同余方程组
形如
{x≡a1(modn1)x≡a2(modn2)x≡a3(modn3)⋮x≡ak(modnk)
\left\{\begin{aligned}
x & \equiv a_1\pmod {n_1} \\
x & \equiv a_2\pmod {n_2} \\
x & \equiv a_3\pmod {n_3} \\
&\vdots \\
x &\equiv a_k\pmod {n_k} \\
\end{aligned}
\right.
⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧xxxx≡a1(modn1)≡a2(modn2)≡a3(modn3)⋮≡ak(modnk)
的方程被称作一元线性同余方程组。中国剩余定理可以用来解决一元线性同余方程组中 n1,n2,...,nkn_1, n_2,..., n_kn1,n2,...,nk 两两互质的情况。
二、中国剩余定理
中国剩余定理的前提条件是 n1,n2,n3,...,nkn_1,n_2,n_3,...,n_kn1,n2,n3,...,nk 两两互质。
一元线性同余方程的解不唯一,任意解 ±y×lcm(n1,n2,...,nk),y≠0\pm y\times lcm(n_1,n_2,...,n_k),y\neq 0±y×lcm(n1,n2,...,nk),y=0 后依然是方程的解。在中国剩余定理的前提条件下,lcm(n1,n2,...,nk)=Πi=1knilcm(n_1,n_2,...,n_k)=\Pi_{i=1}^{k}n_ilcm(n1,n2,...,nk)=Πi=1kni ,因此找到满足条件的最小解就可以得到解的通式。
中国剩余定理的算法流程如下:
- 计算所有 nin_ini 的乘积,记为 nnn
- 对于第 iii 个方程:
a. 计算 mi=nnim_i=\frac{n}{n_i}mi=nin
b. 计算 mim_imi 在模 nin_ini 意义下的逆元 mi−1m_i^{-1}mi−1
c. 计算 ci=mimi−1c_i=m_im_i^{-1}ci=mimi−1 (不要对 nin_ini 取模) - 方程的唯一解为 x=Σi=1kaicimod nx=\Sigma_{i=1}^{k}a_ic_i\mod nx=Σi=1kaicimodn
三、正确性证明
中国剩余定理是构造性的算法,因此证明最终结果代入原方程均满足即可。
证明:
对于 i=1,2,3,...,ki=1,2,3,...,ki=1,2,3,...,k ,有:
① 当 i≠ji\neq ji=j 时,
ajcj=ajmjmj−1=ajnnia_jc_j=a_jm_jm_j^{-1}=a_j\frac{n}{n_i}ajcj=ajmjmj−1=ajnin
此时 nni\frac{n}{n_i}nin 一定是 njn_jnj 的倍数,因此 ajcj≡0(modni)a_jc_j\equiv 0\pmod {n_i}ajcj≡0(modni) 。
② 当 i=ji=ji=j 时,
ajcj=aici=aimimi−1≡ai(modni)
a_jc_j = a_ic_i=a_im_im_i^{-1}\equiv a_i\pmod{n_i}
ajcj=aici=aimimi−1≡ai(modni)
因此
x=Σj=1kajcj=Σj=1,j≠ikajcj+aici≡ai(modni)
x=\Sigma_{j=1}^{k}a_jc_j=\Sigma_{j=1,j\neq i}^{k}a_jc_j+a_ic_i\equiv a_i\pmod{n_i}
x=Σj=1kajcj=Σj=1,j=ikajcj+aici≡ai(modni)
四、扩展中国剩余定理
在中国剩余定理中,我们要求模数两两互质,但往往这种条件过于苛刻,因此考虑模数并不两两互质的情况。仅有一个同余防城港时,退化成线性同余方程,可以使用扩展欧几里得求解,因此从两个同余方程开始。即:
{x≡a1(modn1)x≡a2(modn2)
\left\{\begin{aligned}
x\equiv a_1\pmod{n_1} \\
x\equiv a_2\pmod{n_2}
\end{aligned}\right .
{x≡a1(modn1)x≡a2(modn2)
使用类似线性同余方程的处理方式,得到:
{x=k1n1+a1x=k2n2+a2(1)
\left\{
\begin{aligned}
x = k_1n_1+a_1 \tag{1} \\
x = k_2n_2 + a_2
\end{aligned}
\right.
{x=k1n1+a1x=k2n2+a2(1)
两式做差,得:
n1k1−n2k2=a1−a2
n_1k_1 - n_2k_2 = a_1-a_2
n1k1−n2k2=a1−a2
其中,k1,k2k_1,k_2k1,k2 是未知量,上式是一个二元一次方程,使用扩展欧几里得算法求解得到一组特解 k1∗,k2∗k_1^*,k_2^*k1∗,k2∗ ,写出 k1,k2k_1,k_2k1,k2 的通解形式为:
{k1=k1∗+n2gcd(n1,n2)tk2=k2∗+n1gcd(n1,n2)t
\left\{
\begin{aligned}
k_1=k_1^*+\frac{n_2}{\gcd(n_1,n_2)}t \\
k_2=k_2^*+\frac{n_1}{\gcd(n_1,n_2)}t
\end{aligned}
\right.
⎩⎪⎪⎨⎪⎪⎧k1=k1∗+gcd(n1,n2)n2tk2=k2∗+gcd(n1,n2)n1t
将上式代入公式 (1)(1)(1) 中的任意一个式子可得(代入两个式子得到的是一样的结果,这里代入第一个式子):
x=(k1∗+n2gcd(n1,n2)t)n1+a1=k1∗n1+n1n2gcd(n1,n2)t+a1
x=(k_1^*+\frac{n_2}{\gcd(n_1,n_2)}t)n_1+a_1 =k_1^*n_1+\frac{n_1n_2}{\gcd(n_1,n_2)}t+a_1
x=(k1∗+gcd(n1,n2)n2t)n1+a1=k1∗n1+gcd(n1,n2)n1n2t+a1
这里记 k1∗n1+a1k_1^*n_1+a_1k1∗n1+a1 为 x0x_0x0 ,得:
x=x0+t∗lcm(n1,n2)
x=x_0+t*lcm(n_1,n_2)
x=x0+t∗lcm(n1,n2)
即
x≡x0(modlcm(n1,n2))
x\equiv x_0 \pmod {lcm(n_1,n_2)}
x≡x0(modlcm(n1,n2))
根据以上步骤,可以将两个同余方程合并成一个,由此推广到多个同余方程的情况,即将多个同余方程两两合并直到剩余一个方程。
五、例题和代码
(1)洛谷P1495 【模板】中国剩余定理(CRT)/ 曹冲养猪
中国剩余定理模板题。
#include<bits/stdc++.h>
using namespace std;
const int maxn = 15;
int N;
int n[maxn], a[maxn];
long long c[maxn];
void exgcd(long long a, long long b, long long &x, long long &y){
if (!b){
x = 1, y = 0;
return ;
}
exgcd(b, a % b, y, x);
y -= (a / b) * x;
}
long long get_inv(long long a, long long p){
long long x, y;
exgcd(a, p, x, y);
if (x < 0) x += p;
return x;
}
int main(){
scanf("%d", &N);
long long muln = 1;
for (int i = 1; i <= N; i ++) scanf("%d%d", &n[i], &a[i]), muln *= n[i];
for (int i = 1; i <= N; i ++){
long long m = muln / n[i];
long long m_inv = get_inv(m, n[i]);
c[i] = m * m_inv;
}
long long ans = 0;
for (int i = 1; i <= N; i ++) ans = (ans + a[i] * c[i]) % muln;
printf("%lld\n", ans);
return 0;
}
(2)洛谷P4777 【模板】扩展中国剩余定理(EXCRT)
扩展中国剩余定理模板题。
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e5 + 5;
int N;
long long a[maxn], n[maxn];
__int128 gcd(__int128 a, __int128 b){
return (b == 0) ? a : gcd(b, a % b);
}
void exgcd(__int128 a, __int128 b, __int128 &x, __int128 &y){
if (!b){
x = 1, y = 0;
return ;
}
exgcd(b, a % b, y, x);
y -= (a / b) * x;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);
cin >> N;
for (int i = 1; i <= N; i ++) cin >> n[i] >> a[i];
long long n1 = n[1], a1 = a[1];
for (int i = 2; i <= N; i ++){
// 方程为 k1n1 - kin[i] = a[i] - a1
__int128 d = gcd(n1, n[i]);
if ((a[i] - a1) % d != 0) return -1;
__int128 k1_star, ki_star;
exgcd(n1, n[i], k1_star, ki_star);
k1_star = ((k1_star * (a[i] - a1) / d) % (n[i] / d) + (n[i] / d)) % (n[i] / d);
__int128 x_0 = k1_star * n1 + a1;
a1 = x_0;
n1 = n1 * (n[i] / d);
}
cout << a1 << endl;
return 0;
}
上面的代码是按照之前的推导一步步写的,变量名和推导时一样,所以显得不太美观,如果想找一个便于使用的板子,建议使用下一份的代码。
需要注意的是在执行完 exgcd,将其解扩大若干倍时,可能出现爆 long long 的情况,使用 __int128 可以解决这个问题,当然使用龟速乘也是可以的,代码如下。
#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e5 + 5;
int N;
long long a[maxn], n[maxn];
long long sMul(long long a, long long b, long long MOD){
long long ans = 0;
while (b){
if (b & 1) ans = (ans + a) % MOD;
a = (a + a) % MOD;
b >>= 1;
}
return (ans + MOD) % MOD;
}
long long gcd(long long a, long long b){
return (b == 0) ? a : gcd(b, a % b);
}
void exgcd(long long a, long long b, long long &x, long long &y){
if (!b){
x = 1, y = 0;
return ;
}
exgcd(b, a % b, y, x);
y -= (a / b) * x;
}
long long exCRT(){
long long n1 = n[1], a1 = a[1];
for (int i = 2; i <= N; i ++){
long long n2 = n[i], a2 = a[i];
long long d = gcd(n1, n2);
if ((a2 - a1) % d != 0) return -1;
long long x, y;
exgcd(n1, n2, x, y);
x = (sMul(x, ((a2 - a1) / d % (n2 / d) + (n2 / d)) % (n2 / d), (n2 / d)) + (n2 / d)) % (n2 / d);
a1 = x * n1 + a1;
n1 = n1 * (n2 / d); // 括号必须加,细节防爆
}
return a1;
}
int main(){
cin >> N;
for (int i = 1; i <= N; i ++) cin >> n[i] >> a[i];
cout << exCRT() << endl;
return 0;
}
上面代码使用龟速乘,比使用 __int128 更通用一点,并且对 exCRT 的计算部分进行了封装,比较适合当模板。