一、一元线性同余方程组
形如
{
x
≡
a
1
(
m
o
d
n
1
)
x
≡
a
2
(
m
o
d
n
2
)
x
≡
a
3
(
m
o
d
n
3
)
⋮
x
≡
a
k
(
m
o
d
n
k
)
\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)
的方程被称作一元线性同余方程组。中国剩余定理可以用来解决一元线性同余方程组中
n
1
,
n
2
,
.
.
.
,
n
k
n_1, n_2,..., n_k
n1,n2,...,nk 两两互质的情况。
二、中国剩余定理
中国剩余定理的前提条件是
n
1
,
n
2
,
n
3
,
.
.
.
,
n
k
n_1,n_2,n_3,...,n_k
n1,n2,n3,...,nk 两两互质。
一元线性同余方程的解不唯一,任意解
±
y
×
l
c
m
(
n
1
,
n
2
,
.
.
.
,
n
k
)
,
y
≠
0
\pm y\times lcm(n_1,n_2,...,n_k),y\neq 0
±y×lcm(n1,n2,...,nk),y=0 后依然是方程的解。在中国剩余定理的前提条件下,
l
c
m
(
n
1
,
n
2
,
.
.
.
,
n
k
)
=
Π
i
=
1
k
n
i
lcm(n_1,n_2,...,n_k)=\Pi_{i=1}^{k}n_i
lcm(n1,n2,...,nk)=Πi=1kni ,因此找到满足条件的最小解就可以得到解的通式。
中国剩余定理的算法流程如下:
- 计算所有 n i n_i ni 的乘积,记为 n n n
- 对于第
i
i
i 个方程:
a. 计算 m i = n n i m_i=\frac{n}{n_i} mi=nin
b. 计算 m i m_i mi 在模 n i n_i ni 意义下的逆元 m i − 1 m_i^{-1} mi−1
c. 计算 c i = m i m i − 1 c_i=m_im_i^{-1} ci=mimi−1 (不要对 n i n_i ni 取模) - 方程的唯一解为 x = Σ i = 1 k a i c i m o d n x=\Sigma_{i=1}^{k}a_ic_i\mod n x=Σi=1kaicimodn
三、正确性证明
中国剩余定理是构造性的算法,因此证明最终结果代入原方程均满足即可。
证明:
对于
i
=
1
,
2
,
3
,
.
.
.
,
k
i=1,2,3,...,k
i=1,2,3,...,k ,有:
① 当
i
≠
j
i\neq j
i=j 时,
a
j
c
j
=
a
j
m
j
m
j
−
1
=
a
j
n
n
i
a_jc_j=a_jm_jm_j^{-1}=a_j\frac{n}{n_i}
ajcj=ajmjmj−1=ajnin
此时
n
n
i
\frac{n}{n_i}
nin 一定是
n
j
n_j
nj 的倍数,因此
a
j
c
j
≡
0
(
m
o
d
n
i
)
a_jc_j\equiv 0\pmod {n_i}
ajcj≡0(modni) 。
② 当
i
=
j
i=j
i=j 时,
a
j
c
j
=
a
i
c
i
=
a
i
m
i
m
i
−
1
≡
a
i
(
m
o
d
n
i
)
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
=
1
k
a
j
c
j
=
Σ
j
=
1
,
j
≠
i
k
a
j
c
j
+
a
i
c
i
≡
a
i
(
m
o
d
n
i
)
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
≡
a
1
(
m
o
d
n
1
)
x
≡
a
2
(
m
o
d
n
2
)
\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
=
k
1
n
1
+
a
1
x
=
k
2
n
2
+
a
2
(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)
两式做差,得:
n
1
k
1
−
n
2
k
2
=
a
1
−
a
2
n_1k_1 - n_2k_2 = a_1-a_2
n1k1−n2k2=a1−a2
其中,
k
1
,
k
2
k_1,k_2
k1,k2 是未知量,上式是一个二元一次方程,使用扩展欧几里得算法求解得到一组特解
k
1
∗
,
k
2
∗
k_1^*,k_2^*
k1∗,k2∗ ,写出
k
1
,
k
2
k_1,k_2
k1,k2 的通解形式为:
{
k
1
=
k
1
∗
+
n
2
gcd
(
n
1
,
n
2
)
t
k
2
=
k
2
∗
+
n
1
gcd
(
n
1
,
n
2
)
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
=
(
k
1
∗
+
n
2
gcd
(
n
1
,
n
2
)
t
)
n
1
+
a
1
=
k
1
∗
n
1
+
n
1
n
2
gcd
(
n
1
,
n
2
)
t
+
a
1
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
这里记
k
1
∗
n
1
+
a
1
k_1^*n_1+a_1
k1∗n1+a1 为
x
0
x_0
x0 ,得:
x
=
x
0
+
t
∗
l
c
m
(
n
1
,
n
2
)
x=x_0+t*lcm(n_1,n_2)
x=x0+t∗lcm(n1,n2)
即
x
≡
x
0
(
m
o
d
l
c
m
(
n
1
,
n
2
)
)
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 的计算部分进行了封装,比较适合当模板。