搬运自远古的洛咕博客,故文风与现在有很大不同
CRT 中国剩余定理
中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组(其中 (n_1, n_2, cdots, n_k) 两两互质):
[egin{cases}
x &equiv a_1 pmod {n_1} \
x &equiv a_2 pmod {n_2} \
&vdots \
x &equiv a_k pmod {n_k} \
end{cases}
]
算法流程
- 计算所有模数的积 (n) ;
- 对于第 (i) 个方程:
- 计算 (m_i=frac{n}{n_i}) ;
- 计算 (m_i) 在模 (n_i) 意义下的逆元 (m_i^{-1}) ;
- 计算 (c_i=m_im_i^{-1}) ( 不要对 (n_i) 取模 )。
- 方程组的唯一解为: (a=sum_{i=1}^k a_ic_i pmod n) 。
证明:
对于任意一组同余方程 (i,j) ,其中 (i ot=j) 。因为 (m_i) 中含有 (n_j) ,所以有:
[a_ic_iequiv a_im_im_i^{-1}equiv 0pmod {n_j}
]
另外有:
[a_ic_iequiv a_ipmod {n_i}
]
则对于任意同余方程 (i) ,有:
[sum_{i=1}^ka_ic_iequiv a_ipmod {n_i}
]
证毕。
模版:
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define lxl long long
using namespace std;
lxl a[15],b[15];
int n;
inline lxl exgcd(lxl a,lxl b,lxl &x,lxl &y)
{
if(!b) {x=1,y=0;return a;}
lxl k=exgcd(b,a%b,x,y);
lxl z=x;x=y,y=z-a/b*y;
return k;
}
inline lxl china()
{
lxl M=1,ans=0;
for(int i=1;i<=n;i++)
M*=a[i];
for(int i=1;i<=n;i++)
{
lxl tx,y,Mi=M/a[i];
exgcd(Mi,a[i],tx,y);
ans=(ans+b[i]*Mi*tx)%M;
}
return (ans+M)%M;
}
int main()
{
//freopen("P1495.in","r",stdin);
scanf("%d",&n);
for(int i=1;i<=n;i++)
scanf("%lld%lld",&a[i],&b[i]);
printf("%lld
",china());
return 0;
}
应用
某些计数问题或数论问题出于加长代码、增加难度、或者是一些其他不可告人的原因,给出的模数: 不是质数 !
但是对其质因数分解会发现它没有平方因子,也就是该模数是由一些不重复的质数相乘得到。
那么我们可以分别对这些模数进行计算,最后用 CRT 合并答案。
exCRT 扩展中国剩余定理
说是中国剩余定理,但是好像和CRT关系不大。
扩展中国剩余定理就是合并线性同余方程式,解决 (m_i) 不互素的情况。
先考虑两个同余方程式:
[egin{cases}
x equiv a_1 ({
m {mod}} b_1) \
x equiv a_2 ({
m {mod}} b_2)
end{cases}
implies\
egin{cases}
x = a_1 +k_1 * b_1 \
x equiv a_2 + k_2 * b_2
end{cases}
implies\
a_1 +k_1 * b_1=a_2 + k_2 * b_2
implies\
k_1 * b_1+k_2 * b_2=a_2-a_1
]
如果不定方程有解,使用扩展欧几里德算法求出一个 (k_1) 的特解 (k_0),则 (x) 的一个特解 (x_0=a_1+k_0 * b_1)。通解
[k_i=k_0+u* frac {b_2}{{
m{gcd}}(b_1,b_2)},u in {
m{Z}}
]
则:
[egin{aligned}
x&=k_i * b_1 +a_1\
&=u * frac {b_1b_2}{{
m{gcd}}(b_1,b_2)}+a_1+k_0 * b_1\
&=u * {
m {lcm}}(b_1,b_2)+x_0\
x& equiv x_0 ({
m {mod}} {
m {lcm}} (b_1,b_2))\
end{aligned}
]
于是就把两个同余方程式合并成了一个。同理,将所有方程式按此方法合并,答案为最后的 (x_0)。
模版:
#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#define lxl long long
#define maxn 100005
using namespace std;
inline lxl times(lxl a,lxl b,lxl mod)
{
lxl ans=0;
while(b>0)
{
if(b%2) ans=(ans+a)%mod;
a=(a+a)%mod;
b>>=1;
}
return ans;
}
inline lxl exgcd(lxl a,lxl b,lxl &x,lxl &y)
{
if(!b) {x=1,y=0;return a;}
lxl k=exgcd(b,a%b,x,y);
lxl z=x;x=y,y=z-a/b*y;
return k;
}
int n;
lxl a[maxn],b[maxn];
int main()
{
//freopen("P4777.in","r",stdin);
scanf("%d",&n);
for(int i=1;i<=n;i++) scanf("%lld%lld",&a[i],&b[i]);
lxl M=a[1],ans=b[1];
for(int i=2;i<=n;i++)
{
lxl a1=M,a2=a[i],bi=(b[i]-ans%a2+a2)%a2,x,y;
lxl g=exgcd(a1,a2,x,y);
a2/=g,bi/=g;
x=times(x,bi,a2);
ans+=M*x;
M*=a[i]/g;
ans=(ans+M)%M;
}
printf("%lld
",(ans+M)%M);
return 0;
}