中国剩余定理CRT
更新日志
2024/12/09:开工。2024/12/10:完成证明、添加模板和例题
概念
就不写引入部分了,直接进入正题。
中国剩余定理用来求解形如下式的方程的一组特解:
\[\begin{cases}
x\equiv r_1\pmod {m_1}\\
x\equiv r_2\pmod {m_2}\\
x\equiv r_3\pmod {m_3}\\
\dots\\
x\equiv r_n\pmod {m_n}
\end{cases}
\]
其中 \(m\) 两两互质。
如果不互质,见拓展中国剩余定理。
思路
考虑构造法,构造出一组可行解即可。
事实上,先人已经为我们构造出了一种简单的 \(O(n)\) 方案,所以我们直接记即可。
为了简洁,这里记 \(M\gets\prod\limits_{j=1}^nm_i\)。
\[x=\sum_{i=1}^n[r_i\cdot\frac{M}{m_i}\cdot(\frac{M}{m_i})^{-1}]
\]
其中 \((\frac{M}{m_i})^{-1}\) 表示 \(\frac{M}{m_i}\) 在模 \(m_i\) 意义下的逆元。
证明
我们证明 \(\forall i\in[1,n],x\equiv r_i\pmod {m_i}\) 即可。
\[\because\forall j\ne i,\frac{M}{m_j}\equiv 0\pmod {m_i}
\]
\[\begin{align*}
\forall i\in[1,n],x&\equiv \sum_{j=1}^n[r_j\cdot\frac{M}{m_j}\cdot(\frac{M}{m_j})^{-1}]&\pmod {m_i}\\
x&\equiv r_i\cdot\frac{M}{m_i}\cdot(\frac{M}{m_i})^{-1}&\pmod {m_i}\\
x&\equiv r_i&\pmod{m_i}
\end{align*}
\]
证毕。
常见要求
求通解
不难发现 \(\forall i\in[1,n],M\equiv0\pmod{m_i}\)。
故而通解为 \(x_0+kM,k\in Z\)。
求最小正整数解
我们根据通解知道 \(x_0\) 加减任意个 \(M\) 仍然为解。
直接 \(x_0\bmod M\) 即可。若 \(\le0\),\(+M\) 就行了。
模板
int n;
int r[N],m[N];
int exgcd(int a,int b,int &x,int &y){
if(!b){
x=1,y=0;
return a;
}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
int CRT(){
int M=1;rep(i,1,n)M*=m[i];
int res=0,inv,tmp;
rep(i,1,n){
exgcd(M/m[i],m[i],inv,tmp);
res+=r[i]*(M/m[i])*inv;
}
return res;
}
请注意数据范围,如数据范围很大,开启long long
且计算答案时实时%M
,保险起见乘法过程中强制转换__int128
。
下面是一个修改示例:
int n;
ll r[N],m[N];
ll exgcd(ll a,ll b,ll &x,ll &y){
if(!b){
x=1,y=0;
return a;
}
ll d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
ll CRT(){
ll M=1;rep(i,1,n)M*=m[i];
ll res=0,inv,tmp;
rep(i,1,n){
exgcd(M/m[i],m[i],inv,tmp);
res=((res+(__int128)r[i]*(M/m[i])*inv%M)%M+M)%M;
}
return res;
}
例题
代码
前注:非题解,不做详细讲解
求中间有一段开了__int128
。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef __int128 i128;
typedef double db;
typedef long double ld;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;
template <typename Type>
using vec=vector<Type>;
template <typename Type>
using grheap=priority_queue<Type>;
template <typename Type>
using lrheap=priority_queue<Type,vector<Type>,greater<Type> >;
#define fir first
#define sec second
#define pub push_back
#define pob pop_back
#define puf push_front
#define pof pop_front
#define chmax(a,b) a=max(a,b)
#define chmin(a,b) a=min(a,b)
#define rep(i,x,y) for(int i=x;i<=y;i++)
#define per(i,x,y) for(int i=x;i>=y;i--)
const int inf=0x3f3f3f3f;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int mod=1e9+7/*998244353*/;
const int N=15;
int n;
ll r[N],m[N];
ll exgcd(ll a,ll b,ll &x,ll &y){
if(!b){
x=1,y=0;
return a;
}
ll d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
ll CRT(){
ll M=1;rep(i,1,n)M*=m[i];
ll res=0,inv,tmp;
rep(i,1,n){
exgcd(M/m[i],m[i],inv,tmp);
res=((res+(__int128)r[i]*(M/m[i])*inv%M)%M+M)%M;
}
return res;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
rep(i,1,n)cin>>m[i]>>r[i];
cout<<CRT();
return 0;
}