CRT
用于求解线性同余方程组 ,其中模数两两互质 .
引理:
1.若 \(x \equiv 0 \pmod p\) 且 \(y \equiv 0 \pmod p\) ,则有 \(x+y \equiv 0 \pmod p\)
2.若 \(x \equiv b \pmod p\) 且 ,\(y \equiv 0 \pmod p\) ,则有 \(x+y \equiv b \pmod p (0\leq b<p)\)
则整个方程组可以写为
\(b_1\times \begin{bmatrix}1\\0\\0\end{bmatrix} + \cdots+b_i\times \begin{bmatrix}0\\1\\0\end{bmatrix}\) + \(\cdots\) +\(b_n\times \begin{bmatrix}0\\0\\1\end{bmatrix}\)
也就是说,想要求出最终方程组的解就只需要求出上面每个方程的解。
考虑每个式子中,因为每两个模数之间两两互质,设 \(N=\sum_{i=1}^nm_i\)
则 \(gcd(\frac{N}{m_i} , m_i)=1\) ,\(\frac{N}{m_i}\) 即为其他模数的 lcm。
也为其他 \(n-1\) 个方程的公共解,则第 \(i\) 个方程的解可以写成 \(x_i=\frac{N}{m_i} \times y\)
有 \(\frac{N}{m_i}\times y \equiv 1\pmod {m_i}\)
等价于 \(\frac{N}{m_i}\times y+m_i\times z=1\) ,可以用扩展欧几里得求解。
求出 \(x_i\) 后 ,有 \(b_i\times x_i \equiv b_i \pmod {m_i}\)
最后整个方程组的解就为 \(X = \sum_{i=1}^n b_i \times x_i \pmod N\)
Code:
#include<stdio.h>
#define N 100
long long W[N],B[N],nn,T;
inline void exgcd(long long a,long long b,long long &x,long long &y){
if(!b){
x=1;y=0;
}else{
exgcd(b,a%b,x,y);
long long t=x;
x=y;
y=t-a/b*y;
}
}
inline long long China(long long k){
long long x,y,a=0,m,n=1;
for(long long i=0;i<k;i++) n*=W[i];
for(long long i=0;i<k;i++){
m=n/W[i];
exgcd(W[i],m,x,y);
a=a+y*B[i]*m;
}
a%=n;
while(a<=0) a+=n;
while(a>=n) a-=n;
return a;
}
signed main(){
scanf("%lld",&nn);
for(long long i=0;i<nn;i++)
scanf("%lld",&B[i]);
for(long long i=0;i<nn;i++)
scanf("%lld",&W[i]);
for(long long i=0;i<nn;i++)
B[i]=(B[i]%W[i]+W[i])%W[i];
printf("%lld",China(nn));
}
EXCRT
其实中国剩余定理和它的扩展版本关系并不大。
模数不是两两互质的,也就是说先前求的 \(N\) 不是它们的 LCM,不能直接用。
我们考虑构造的方法,假设我们已经求出了前 \(k-1\) 个方程的和解,现在需要合并第 \(k\) 个方程。
设前 \(k-1\) 个方程的和解为 \(X\),则他们的通解可以写成 \(X + t \times N\) (其中 \(N\) 为前 \(k-1\) 个模数的 LCM)
则要合并第k个方程,即要解如下方程 $ X + t\times N \equiv b_k \pmod {m_k}$
此方程又等价于 $ t \times N + m_k \times y = b_k-X $
显然可用扩欧。求出 \(t\) 后 ,前 \(k\) 个方程的解就为 \(X^{'}=X+t\times N\),继续更新 \(N\) ,\(N=lcm(N,m_k)\)
Code:
#include<stdio.h>
#define N 100006
#define ll long long
ll a[N],b[N],ans,B,M,GCD,x,y;
int n;
template<class T>
inline void read(T &x){
x=0;T flag=1;char c=getchar();
while(c<'0'||c>'9'){if(c=='-') flag=-1;c=getchar();}
while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+c-48;c=getchar();}
x*=flag;
}
inline ll exgcd(ll a,ll b,ll &xx,ll &yy){
if(!b){
xx=1;yy=0;
return a;
}
ll gcd=exgcd(b,a%b,xx,yy);
ll tmp=xx;
xx=yy;
yy=tmp-(a/b)*yy;
return gcd;
}
inline ll mul(ll a,ll b,ll Mod){
ll ans=0;
while(b){
if(b&1) ans=(ans+a)%Mod;
a=(a<<1)%Mod;
b>>=1;
}
return ans;
}
int main(){
read(n);
for(int i=1;i<=n;i++)
read(a[i]),read(b[i]);
ans=b[1];
M=a[1];
for(int i=2;i<=n;i++){
B=(b[i]-ans%a[i]+a[i])%a[i];
GCD=exgcd(M,a[i],x,y);
x=mul(x,B/GCD,a[i]);
ans+=M*x;
M*=a[i]/GCD;
ans=(ans%M+M)%M;
}
printf("%lld",ans);
}
/*
3
11 6
25 9
33 17
*/