前言
中国剩余定理(也叫孙子定理)并不是很复杂,由于最近用到了,以前学的时候还不写博客,所以现在补一下
中国剩余定理(CRT)
问题
给出
n
n
n个同余方程
x
≡
a
1
(
m
o
d
p
1
)
x
≡
a
2
(
m
o
d
p
2
)
⋅
⋅
⋅
⋅
⋅
⋅
⋅
⋅
x
≡
a
n
(
m
o
d
p
n
)
\begin{aligned} x&\equiv a_1\pmod{p_1}\\ x&\equiv a_2\pmod{p_2}\\ ··&\ \ \ \ \ ···\ \ \ \ \ \ \ \ ···\\ x&\equiv a_n\pmod{p_n}\\ \end{aligned}
xx⋅⋅x≡a1(modp1)≡a2(modp2) ⋅⋅⋅ ⋅⋅⋅≡an(modpn)
保证所有的
p
p
p两两互质
求最小自然数解
x
x
x
做法
我们令
P
=
∏
i
=
1
n
p
i
P=\prod_{i=1}^n p_i
P=i=1∏npi
容易发现,
x
<
P
x<P
x<P
证明,若现在求出一个
x
x
x满足条件且
x
≥
P
x\ge P
x≥P
那么
(
x
m
o
d
  
P
)
(x\mod P)
(xmodP)一定也满足条件且
(
x
m
o
d
  
p
)
<
P
(x\mod p)<P
(xmodp)<P
我们考虑对于一个限制
x
≡
a
i
(
m
o
d
p
i
)
x\equiv a_i\pmod{p_i}
x≡ai(modpi),我们要进行一种操作,使得其它
n
−
1
n-1
n−1个
p
p
p的模意义下值不变,模
p
i
p_i
pi意义下的值从
0
0
0变为
a
i
a_i
ai
考虑如果我们加上的数是
X
∗
P
p
i
X*\frac{P}{p_i}
X∗piP形式的,那么我们就可以满足“其它
n
−
1
n-1
n−1个
p
p
p的模意义下值不变”,由于保证了所有
p
p
p两两互质,所以
P
p
i
̸
≡
0
(
m
o
d
p
i
)
\frac{P}{p_i}\not\equiv0\pmod{p_i}
piP̸≡0(modpi)
设
P
p
i
≡
v
i
(
m
o
d
p
i
)
\frac{P}{p_i}\equiv v_i\pmod{p_i}
piP≡vi(modpi)
那么只要加上
(
a
i
v
i
−
1
(
m
o
d
P
p
i
)
)
P
p
i
(a_iv_i^{-1}\pmod{\frac{P}{p_i}})\frac{P}{p_i}
(aivi−1(modpiP))piP就可以满足这条限制
最后求出
x
x
x后再模
P
P
P就好了
例题
竟然有模板题。。。
洛谷P3868 [TJOI2009]猜数字
有个要注意的地方是这里的逆元应用扩展欧几里得(exGCD)求,另外由于模数过大,需要用到慢速乘(复杂度
Θ
(
l
o
g
n
)
\Theta(logn)
Θ(logn)),可以用杜教乘优化(复杂度
Θ
(
1
)
\Theta(1)
Θ(1),我的这篇博客里有提到)
代码
ac代码(其实也是板子)
#include<cstdio>
#include<cctype>
#include<algorithm>
namespace fast_IO
{
const int IN_LEN=10000000,OUT_LEN=10000000;
char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;
inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}
inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}
inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
}
using namespace fast_IO;
#define getchar() getchar_()
#define putchar(x) putchar_((x))
#define rg register
typedef long long LL;
template <typename T> inline T max(const T a,const T b){return a>b?a:b;}
template <typename T> inline T min(const T a,const T b){return a<b?a:b;}
template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;}
template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;}
template <typename T> inline T abs(const T a){return a>0?a:-a;}
template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);}
template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;}
template <typename T> inline T square(const T x){return x*x;};
template <typename T> inline void read(T&x)
{
char cu=getchar();x=0;bool fla=0;
while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
if(fla)x=-x;
}
template <typename T> inline void printe(const T x)
{
if(x>=10)printe(x/10);
putchar(x%10+'0');
}
template <typename T> inline void print(const T x)
{
if(x<0)putchar('-'),printe(-x);
else printe(x);
}
LL EXgcd(const LL a,const LL b,LL &x,LL &y)
{
if(!b)
{
x=1,y=0;
return a;
}
const LL res=EXgcd(b,a%b,y,x);
y-=a/b*x;
return res;
}
inline LL msc(LL a,LL b,LL mod)
{
LL v=(a*b-(LL)((long double)a/mod*b+1e-8)*mod);
return v<0?v+mod:v;
}
int n,a[11],p[11];
LL CRT()
{
LL P=1,sum=0;
for(rg int i=1;i<=n;i++)P*=p[i];
for(rg int i=1;i<=n;i++)
{
const LL m=P/p[i];
LL x,y;
EXgcd(p[i],m,x,y);
sum=(sum+msc(msc(y,m,P),a[i],P))%P;
}
return sum;
}
int main()
{
read(n);
for(rg int i=1;i<=n;i++)read(a[i]);
for(rg int i=1;i<=n;i++)
{
read(p[i]);
a[i]=(a[i]%p[i]+p[i])%p[i];
}
print(CRT());
return flush(),0;
}
扩展中国剩余定理(exCRT)
问题
给出
n
n
n个同余方程
x
≡
a
1
(
m
o
d
p
1
)
x
≡
a
2
(
m
o
d
p
2
)
⋅
⋅
⋅
⋅
⋅
⋅
⋅
⋅
x
≡
a
n
(
m
o
d
p
n
)
\begin{aligned} x&\equiv a_1\pmod{p_1}\\ x&\equiv a_2\pmod{p_2}\\ ··&\ \ \ \ \ ···\ \ \ \ \ \ \ \ ···\\ x&\equiv a_n\pmod{p_n}\\ \end{aligned}
xx⋅⋅x≡a1(modp1)≡a2(modp2) ⋅⋅⋅ ⋅⋅⋅≡an(modpn)
求最小自然数解
x
x
x
(与上面的问题的不同在于这里不保证所有的
p
p
p两两互质)
做法
考虑我们现在已经符合了前
k
−
1
k-1
k−1个限制了
设
P
=
l
c
m
(
p
1
,
p
2
⋅
⋅
⋅
p
k
−
1
)
P=lcm(p_1,p_2···p_k-1)
P=lcm(p1,p2⋅⋅⋅pk−1),容易发现,满足前
k
−
1
k-1
k−1个限制的最小的自然数
x
<
P
x<P
x<P
证明同理,也就是现在答案是
x
+
t
P
x+tP
x+tP形式的
我们要满足第
k
k
k个限制,其实很方便,就是找到一个
t
t
t满足
x
+
t
P
≡
a
k
(
m
o
d
p
k
)
x+tP\equiv a_k\pmod{p_k}
x+tP≡ak(modpk)
移项
t
P
≡
a
k
−
x
(
m
o
d
p
k
)
tP\equiv a_k-x\pmod{p_k}
tP≡ak−x(modpk)
转化
t
P
+
q
p
k
≡
a
k
−
x
tP+qp_k\equiv a_k-x
tP+qpk≡ak−x
然后直接拓展欧几里得求解
若无解输出无解
否则直接将答案更新
例题
洛谷P4777 【模板】扩展中国剩余定理(EXCRT)
这道题保证有解,所以不用判无解了,保证了解的范围,所以容易发现不用高精度
代码
贴出ac代码,代码里有句注释,是用来判无解的
为了卡常才把它注释掉的
#include<cstdio>
#include<cctype>
#include<algorithm>
namespace fast_IO
{
const int IN_LEN=10000000,OUT_LEN=10000000;
char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;
inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}
inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}
inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
}
using namespace fast_IO;
#define getchar() getchar_()
#define putchar(x) putchar_((x))
#define rg register
typedef long long LL;
template <typename T> inline T max(const T a,const T b){return a>b?a:b;}
template <typename T> inline T min(const T a,const T b){return a<b?a:b;}
template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;}
template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;}
template <typename T> inline T abs(const T a){return a>0?a:-a;}
template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);}
template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;}
template <typename T> inline T square(const T x){return x*x;};
template <typename T> inline void read(T&x)
{
char cu=getchar();x=0;bool fla=0;
while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
if(fla)x=-x;
}
template <typename T> inline void printe(const T x)
{
if(x>=10)printe(x/10);
putchar(x%10+'0');
}
template <typename T> inline void print(const T x)
{
if(x<0)putchar('-'),printe(-x);
else printe(x);
}
LL EXgcd(const LL a,const LL b,LL &x,LL &y)
{
if(!b)
{
x=1,y=0;
return a;
}
const LL res=EXgcd(b,a%b,y,x);
y-=a/b*x;
return res;
}
inline LL msc(LL a,LL b,LL mod)
{
LL v=(a*b-(LL)((long double)a/mod*b+1e-8)*mod);
return v<0?v+mod:v;
}
LL n,a[100001],p[100001];
LL EXCRT()
{
LL P=1,sum=0;
for(rg int i=1;i<=n;i++)
{
const LL o=((a[i]-sum)%p[i]+p[i])%p[i];
//if(o%gcd(P,p[i])!=0)return -1;
const LL xs=o/gcd(P,p[i]);
LL x,y;
EXgcd(p[i],P,x,y);
const LL M=P;P=P/gcd(P,p[i])*p[i];
sum=(sum+msc(msc(y,M,P),xs,P))%P;
}
return sum;
}
int main()
{
read(n);
for(rg int i=1;i<=n;i++)
{
read(p[i]),read(a[i]);
a[i]=(a[i]%p[i]+p[i])%p[i];
}
print(EXCRT());
return flush(),0;
}
总结
并不是很难的两个算法,学起来很方便
后记:感觉exCRT严格优于CRT,有没有大佬指出CRT有什么用啊