拓展中国剩余定理ExCRT

更新日志

2024/12/11：开工。

2024/12/12：添加模板与例题代码。

---

概念

同中国剩余定理，但模数不满足两两互质。

求解。

思路

我们先考虑两个方程如何解决。

\[\begin{cases}x\equiv r_1\pmod {m_1}\\x\equiv r_2\pmod {m_2}\end{cases}\\ \Rightarrow x=m_1p+r_1=m_2q+r_2\\ \Leftrightarrow m_1p-m_2q=r_2-r_1 \]

其中 \(m_1,m_2,r_1,r_2\) 均为已知量，所以这就是一个普通的线性不定方程，使用 \(\mathrm{exgcd}\) 解决即可。

注意判断是否有解，详见 \(\mathrm{exgcd}\)。

然后我们考虑多个方程组，两两合并即可。

具体的，我们将解决的两个方程组合并为一个新的同余方程。一个可行的例子如下：

\[x\equiv m_1p+r_1\pmod{\mathrm{lcm}(m_1,m_2)} \]

显然等价于原方程组。

模板

直接给出较为常用的模板吧。

在计算过程中的 \(x\) 均为不定方程最小正整数解，全部开了 long long。

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 excrt(){
	rep(i,2,n){
		ll x,y,k;
		ll d=abs(exgcd(m[i-1],m[i],x,y));
		if((r[i]-r[i-1])%d)return -1;
		x=(r[i]-r[i-1])/d*x;
		k=abs(m[i]/d);x=(x%k+k)%k;
		m[i]=m[i-1]/d*m[i];
		r[i]=(m[i-1]*x+r[i-1])%m[i];
	}
	return (r[n]%m[n]+m[n])%m[n];
}

小优化

不难发现不定方程的 x 计算过程中可能会爆 long long，两种解决方式。

龟速乘

那么，代价呢？复杂度多出 \(O(\log n)\)。

ll smul(ll a,ll b,ll m){
	a=(a%m+m)%m;
	b=(b%m+m)%m;
	ll res=0;
	while(b){
		if(b&1)(res+=a)%=m;
		(a+=a)%=m;
		b>>=1;
	}
	return res;
}
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 excrt(){
	rep(i,2,n){
		ll x,y,k;
		ll d=abs(exgcd(m[i-1],m[i],x,y));
		if((r[i]-r[i-1])%d)return -1;
		x=smul((r[i]-r[i-1])/d,x,abs(m[i]/d));
		m[i]=m[i-1]/d*m[i];
		r[i]=(m[i-1]*x+r[i-1])%m[i];
	}
	return (r[n]%m[n]+m[n])%m[n];
}

__int128

简单粗暴。

int n;
ll r[N],m[N];
ll exgcd(ll a,ll b,i128 &x,i128 &y){
	if(!b){
		x=1,y=0;
		return a;
	}
	ll d=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return d;
}
ll excrt(){
	rep(i,2,n){
		i128 x,y,k;
		ll d=abs(exgcd(m[i-1],m[i],x,y));
		if((r[i]-r[i-1])%d)return -1;
		x=(r[i]-r[i-1])/d*x;
		k=abs(m[i]/d);x=(x%k+k)%k;
		m[i]=m[i-1]/d*m[i];
		r[i]=(m[i-1]*x+r[i-1])%m[i];
	}
	return (r[n]%m[n]+m[n])%m[n];
}

例题

LG4777

代码

龟速乘版本

#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=1e5+5;

ll smul(ll a,ll b,ll m){
	a=(a%m+m)%m;
	b=(b%m+m)%m;
	ll res=0;
	while(b){
		if(b&1)(res+=a)%=m;
		(a+=a)%=m;
		b>>=1;
	}
	return res;
}
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 excrt(){
	rep(i,2,n){
		ll x,y,k;
		ll d=abs(exgcd(m[i-1],m[i],x,y));
		if((r[i]-r[i-1])%d)return -1;
		x=smul((r[i]-r[i-1])/d,x,abs(m[i]/d));
		m[i]=m[i-1]/d*m[i];
		r[i]=(m[i-1]*x+r[i-1])%m[i];
	}
	return (r[n]%m[n]+m[n])%m[n];
}

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<<excrt();
    return 0;
}

__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=1e5+5;

int n;
ll r[N],m[N];
ll exgcd(ll a,ll b,i128 &x,i128 &y){
	if(!b){
		x=1,y=0;
		return a;
	}
	ll d=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return d;
}
ll excrt(){
	rep(i,2,n){
		i128 x,y,k;
		ll d=abs(exgcd(m[i-1],m[i],x,y));
		if((r[i]-r[i-1])%d)return -1;
		x=(r[i]-r[i-1])/d*x;
		k=abs(m[i]/d);x=(x%k+k)%k;
		m[i]=m[i-1]/d*m[i];
		r[i]=(m[i-1]*x+r[i-1])%m[i];
	}
	return (r[n]%m[n]+m[n])%m[n];
}

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<<excrt();
    return 0;
}