Chinese Remainder Theorem(CRT)
又称孙子定理,是用来解决一类同余方程组问题的方法。
这个要用到欧几里得算法相关知识,不会走这里GO
我们从一个简单的例子入手:
现有如下式子( p1,p2 p 1 , p 2 为质数):
①{x≡a1(mod p1)x≡a2(mod p2) ① { x ≡ a 1 ( m o d p 1 ) x ≡ a 2 ( m o d p 2 )
按照同余式,上式可以写成:
②{x=a1+k1p1x=a2+k2p2 ② { x = a 1 + k 1 p 1 x = a 2 + k 2 p 2 其中 k1,k2∈Z k 1 , k 2 ∈ Z
通过联立起来,将其写成 a1+k1p1=a2+k2p2 a 1 + k 1 p 1 = a 2 + k 2 p 2 , k1p1−k2p2=a2−a1 k 1 p 1 − k 2 p 2 = a 2 − a 1
又因为 gcd(p1,p2)=1 g c d ( p 1 , p 2 ) = 1 ,所以此方程肯定有解。
对应 ax+by=c a x + b y = c 式子,那么 ⎧⎩⎨a=p1b=−p2c=a2−a1 { a = p 1 b = − p 2 c = a 2 − a 1 ,所以根据这种形式,我们用 exgcd e x g c d 可以求出一组解 (k^1,k^2) ( k ^ 1 , k ^ 2 ) ,然后我们将其带回原方程组 ② ② 可得,然后通过 exgcd e x g c d 的知识我们可以由这组 (k^1,k^2) ( k ^ 1 , k ^ 2 ) 特解,推导出其它的解: {k1=p2t+k^1k2=p1t+k^2,k1,k2∈Z { k 1 = p 2 t + k ^ 1 k 2 = p 1 t + k ^ 2 , k 1 , k 2 ∈ Z 。
带入原方程求解得:
x=a1+k1p1=a1+p1p2t+k^1=x0+p1p2t
x
=
a
1
+
k
1
p
1
=
a
1
+
p
1
p
2
t
+
k
^
1
=
x
0
+
p
1
p
2
t
等价于
x≡x0(mod p1p2)
x
≡
x
0
(
m
o
d
p
1
p
2
)
。
我们就得到了方程的答案,同时将模数合并了起来。
所以对于多个方程构成的方程组,我们每次选择两个合并出一个新的,然后再用新的去和另一个合并(注意:模数相乘后小心爆范围 long long l o n g l o n g )。
那么对于如下的模方程组:
所有的 mi m i 均互质,所以用 exgcd e x g c d 求出特解即可解决。(LRJ的蓝书上也有讲解)。
模板代码:
\\[TJOI2009猜数字]
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=21;
int n;
ll a[M],b[M];
void exgcd(ll a,ll b,ll &x,ll &y){
if(!b){x=1;y=0;}
else {exgcd(b,a%b,y,x);y-=x*(a/b);}
}
ll mul(ll a,ll b,ll mod){
ll c=0;
for(;b>0;b>>=1,a=(a+a)%mod)if(b&1) c=(c+a)%mod;
return c;
}
ll CRT(ll *a,ll *b,int n){
ll M=1,d=0,x,y;
for(int i=1;i<=n;i++) M*=b[i];
for(int i=1;i<=n;i++){
ll w=M/b[i];
exgcd(b[i],w,x,y);y=(y%b[i]+b[i])%b[i];//exgcd求解
d=(d+mul(mul(w,y,M),(a[i]%b[i]+b[i])%b[i],M))%M;
}
return (d+M)%M;
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)scanf("%lld",&a[i]);
for(int i=1;i<=n;i++)scanf("%lld",&b[i]);
ll ans=CRT(a,b,n);
printf("%lld\n",ans);
return 0;
}
扩展中国剩余定理
我们现在解决的一类模方程的求法,但是这是在模数互质的条件下。大多数时候模数都不会互质,那么应该如何求解呢?
我们仍然从简单的入手
看如下一个方程组:
同样的我们将其变形得到 a1+k1m1=a2+k2m2 a 1 + k 1 m 1 = a 2 + k 2 m 2 ,如果该方程要有整数解,那么根据裴蜀定理,要满足 gcd(m1,m2)|(a2−a1) g c d ( m 1 , m 2 ) | ( a 2 − a 1 ) 。
这里我们先假设用
exgcd
e
x
g
c
d
求出了一组解特解
(k1,k2)
(
k
1
,
k
2
)
,那么通解有是什么呢?
接下来令
g=gcd(m1,m2)
g
=
g
c
d
(
m
1
,
m
2
)
,又因为
g|(a2−a1)
g
|
(
a
2
−
a
1
)
(因为我们先假设它有解),所以原方程可以写成
回带入原方程,我们就可以得到:
lcm为最小公倍数。
实质上就等于解: x≡x0 (mod lcm(m1,m2)) x ≡ x 0 ( m o d l c m ( m 1 , m 2 ) ) 这个方程。
举个例子,解如下方程:
先自己试着解一下再看下面解答。
先联立 ①② ① ② 得 ④x≡ 8(mod 15) ④ x ≡ 8 ( m o d 15 ) ,再联立 ③④ ③ ④ 得到最终答案 x≡23 (mod 105) x ≡ 23 ( m o d 105 ) ,那么答案就为 x=23 x = 23 ,反带回去发现确实成立。
【模板IN】
模板代码:
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=201;
ll xx(ll a,ll b,ll mod){
if(a>b) swap(a,b);ll c=0;
for(ll i=0;(1ll<<i)<=b;i++){if(b&(1ll<<i)){c=(c+(1ll<<i)*a)%mod;}}
return c;
}
void exgcd(ll a,ll b,ll &d,ll &x,ll &y){
if(!b){d=a;x=1;y=0;}
else {exgcd(b,a%b,d,y,x);y-=x*(a/b);}
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll inv(ll a,ll b){
ll g,x,y;
exgcd(a,b,g,x,y);
if(g!=1) return -1;
return (x%b+b)%b;
}
bool merge(ll a1,ll b1,ll a2,ll b2,ll &a,ll &b){
ll g=gcd(b1,b2);
ll c=a2-a1;
if(c%g) return 0;//判断有无解
c=(c%b2+b2)%b2;
b1/=g;b2/=g;
c/=g;c=xx(c,inv(b1,b2),b2);
b=b1*b2*g;
c=xx(c,xx(b1,g,b),b);
c+=a1;
a=(c%b+b)%b;
return 1;
}
ll CRT(ll *a,ll *b,int n){
ll aa=a[1],bb=b[1];
for(int i=2;i<=n;i++){
ll da=a[i],db=b[i];
ll nexa,nexb;
//合并
if(!merge(aa,bb,da,db,nexa,nexb)) return -1;
aa=nexa;bb=nexb;
}
return (aa%bb+bb)%bb;
}
ll a[M],b[M];
int n;
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){scanf("%lld%lld",&b[i],&a[i]);}
ll ans=CRT(a,b,n);
printf("%lld\n",ans);
return 0;
}
//模板2
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define LL long long
using namespace std;
const int N=210;
int n;
LL x,y,lcm,equ[N][2];
LL multi(LL a,LL b,LL p)
{
a=(a%p+p)%p;b=(b%p+p)%p;LL ans=0;
for(;a;a>>=1,b=(b*2)%p)if(a&1)ans=(ans+b)%p;
return ans;
}
LL exgcd(LL a,LL b,LL &x,LL &y)
{
if(!b)
{
x=1,y=0;
return a;
}
LL val=exgcd(b,a%b,x,y);
LL t=x;x=y;y=t-a/b*y;return val;
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++)
scanf("%lld%lld",&equ[i][1],&equ[i][0]);
for(int i=1;i<n;i++)
{
LL val=exgcd(equ[i][1],equ[i+1][1],x,y);
lcm=equ[i][1]/val*equ[i+1][1];
if((equ[i+1][0]-equ[i][0])%val)
{
printf("No Answer\n");
return 0;
}
val=multi(y,(equ[i+1][0]-equ[i][0])/val,lcm);
equ[i+1][0]=(multi(equ[i+1][1],-val,lcm)+equ[i+1][0]+lcm)%lcm;
equ[i+1][1]=lcm;
}
printf("%lld\n",(equ[n][0]%equ[n][1]+equ[n][1])%equ[n][1]);
return 0;
}
一个额外的小技巧:
例如NOI2018的屠龙勇士。
当同余方程组有系数,如下式子:
如果 a,b a , b 都互质,我们用扩展欧几里得定理就可以求得逆元,然后把系数除过去就可以计算了。
但是当 a,b a , b 不一定互质怎么办?
我们观察式子, ax≡b(mod m) a x ≡ b ( m o d m ) ,可以将其写成 ax=km+b a x = k m + b ,我们记 g=gcd(a,m) g = g c d ( a , m ) (gcd表示求最大公约数),那么左右可以写成 (gw1)x=k(gw2)+b ( g w 1 ) x = k ( g w 2 ) + b ,由于左边为 g g 的倍数,右边那么同理也应该为的倍数,所以显然 b b 也应该为的倍数,若不为则该方程组无整数解,有的话将 a,b,m a , b , m 同时除以 g g <script type="math/tex" id="MathJax-Element-155">g</script>,就得到了系数互质的方程,用原来的方法解就可以了。
End
Thanks For Reading!