中国剩余定理 学习笔记
中国剩余定理 (Chinese Remainder Theorem) 简称CRT。
CRT
CRT可以求解一个这样的方程组的最小解
\[\left\{
\begin{array}{**lr**}
x \equiv b_1 \pmod {a_1} \\
x \equiv b_2 \pmod {a_2} \\
x \equiv b_3 \pmod {a_3} \\
\vdots \\
x \equiv b_k \pmod {a_k} \\
\end{array}
\right.
\]
这里保证数组 \(a\) 两两互质。
CRT的原理就是得到 \(c_i\) 使 \(c_i \equiv 1 \pmod {a_i}\) 并且 \(c_i \equiv 0 \pmod {a_j},i\not= j\)
显然答案就是 \(\sum^k_{i=0} b_ic_i\bmod \left(\prod^{n}_{i=1}a_i\right)\)
算法流程如下：
- 计算 \(a=\prod a_i\)
- 对于第 \(i\) 个方程：
a. 计算 \(m_i=\frac{a}{a_i}\)
b. 计算 \(m_i\) 在模 \(a_i\) 的逆元 \(m_i^{-1}\)
c. 计算 \(c_i=m_i\times m_i^{-1}\) - 方程的解为 \(x \equiv \sum^k_{i=1}b_ic_i \pmod a\)
证明略 （显然啊）
代码(洛谷P1495)：
#include<cstdio>
#define maxn 100039
using namespace std;
//#define debug
typedef long long ll;
typedef long long Type;
inline Type read(){
Type sum=0;
int flag=0;
char c=getchar();
while((c<'0'||c>'9')&&c!='-') c=getchar();
if(c=='-') c=getchar(),flag=1;
while('0'<=c&&c<='9'){
sum=(sum<<1)+(sum<<3)+(c^48);
c=getchar();
}
if(flag) return -sum;
return sum;
}
ll n,a[maxn],b[maxn];
ll sn,m,c,ans;
ll gcd(ll a,ll b){
if(a%b==0) return b;
return gcd(b,a%b);
}
inline ll lcm(ll a,ll b){ return a*b/gcd(a,b); }
void exgcd(ll a,ll b,ll &x,ll &y){
if(!b){ x=1; y=0; return; }
exgcd(b,a%b,x,y);
ll t=x;
x=y; y=t-a/b*y;
return;
}
ll js(ll a1,ll b1,ll a2,ll b2){
ll p,q;
exgcd(a1,-a2,p,q);
if((b2-b1)%gcd(a1,-a2)) return -1;
q*=(b2-b1)/gcd(a1,-a2);
p*=(b2-b1)/gcd(a1,-a2);
//if(a1*p+b1 != a2*q+b2) printf("Oops!\n");
//else printf("OK\n");
int r=lcm(a1,a2);
return ((a2*q+b2)%r+r)%r;
}//x=b_i(mod a_i)
ll excrt(){
ll ta,tb;
ta=a[1],tb=b[1];
for(int i=2;i<=n;i++){
tb=js(ta,tb,a[i],b[i]);
if(tb==-1) return -1;
ta=lcm(ta,a[i]);
}
return tb;
}
int main(){
//freopen("1.in","r",stdin);
//freopen("my.out","w",stdout);
n=read(); sn=1;
for(int i=1;i<=n;i++)
a[i]=read(),b[i]=read();//x=b_i(mod a_i)
printf("%lld",excrt());
return 0;
}
exCRT
如果数组 \(a\) 不是两两互质呢？
我们尝试修复一下CRT的方法让它变成exCRT，但是我们发现，CRT的核心都说过了，它是求出 \(c_i\) 使 \(c_i \equiv 1 \pmod {a_i}\) 并且 \(c_i \equiv 0 \pmod {a_j},i\not= j\) ，但是我们发现，我们有时不能找到这样的 \(c_i\) ，所以CRT原来的思路是不能实现的。
看来只能用另外一种方法了。
考虑两个方程的情况。
\[\left\{
\begin{array}{**lr**}
x \equiv b_1 \pmod {a_1} \\
x \equiv b_2 \pmod {a_2} \\
\end{array}
\right.
\]
令 \(x=a_1q+b_1=a2_p+b_2\)
移项，得到 \(a_1q+a_2(-p)=b_2-b_1\)
显然只要exgcd解出 \(p,q\) 的一组解就可以解决问题了，然后显然合并后的方程是：
\[x \equiv a_1q+b_1 \pmod {\operatorname{lcm}\left(a_1,a_2\right)}
\]
然后就可以过 洛谷P4777了，当然要注意精度问题，加上防爆乘，不然会炸精度。。。
#include<cstdio>
#define maxn 100039
using namespace std;
//#define debug
typedef long long ll;
typedef long long Type;
inline Type read(){
Type sum=0;
int flag=0;
char c=getchar();
while((c<'0'||c>'9')&&c!='-') c=getchar();
if(c=='-') c=getchar(),flag=1;
while('0'<=c&&c<='9'){
sum=(sum<<1)+(sum<<3)+(c^48);
c=getchar();
}
if(flag) return -sum;
return sum;
}
ll mol(ll a,ll b,ll p){
int flag=0;
if(a<0) a=-a,flag^=1;
if(b<0) b=-b,flag^=1;
ll ans=0,tmp=a;
while(b){
if(b&1) ans=(ans+tmp)%p;
b>>=1; tmp=(tmp<<1)%p;
}
if(flag) return (-ans%p+p)%p;
return ans%p;
}
ll gcd(ll a,ll b){
if(a%b==0) return b;
return gcd(b,a%b);
}
inline ll lcm(ll a,ll b){ return a/gcd(a,b)*b; }
void exgcd(ll a,ll b,ll &x,ll &y){
if(!b){ x=1; y=0; return; }
exgcd(b,a%b,x,y);
ll t=x;
x=y; y=t-(a/b)*y;
return;
}
int n;
ll a[maxn],b[maxn];
ll excrt(){
ll m=a[1],ans=b[1];
ll p,q;
for(int i=2;i<=n;i++){
exgcd(m,a[i],p,q);
//if((b[i]-ans)%gcd(m,a[i])!=0) return -1;
p=mol(p,(b[i]-ans)/gcd(m,a[i]),a[i]/gcd(m,a[i]));
ans=(mol(p,m,lcm(a[i],m))+ans)%lcm(a[i],m);
m=lcm(a[i],m);
}
return ans%m;
}//x=b(mod a)
int main(){
//freopen("P4777_13.in","r",stdin);
//freopen(".in","w",stdout);
n=read();
for(int i=1;i<=n;i++)
a[i]=read(),b[i]=read();
printf("%lld",excrt());
return 0;
}
