扩展中国剩余定理(EXCRT)

扩展中国剩余定理(EXCRT)

中国剩余定理可以解形如下面的一元线性同余方程组：
$$\{\begin{array}{c}ans\equiv x_1(\mathrm{m}\mathrm{o}\mathrm{d}{a}_1)\\ ans\equiv x_2(\mathrm{m}\mathrm{o}\mathrm{d}{a}_2)\\ ans\equiv x_3(\mathrm{m}\mathrm{o}\mathrm{d}{a}_3)\\ .\\ .\\ .\\ ans\equiv x_n(\mathrm{m}\mathrm{o}\mathrm{d}{a}_n)\end{array}$$
其中a₁,a₂,a₃…a_n两两互素。

如果a₁,a₂,a₃…a_n不两两互素时，就不能利用中国剩余定理去解上述方程组了。因为不存在另外n-1个a_i的乘积关于剩下一个数的逆元。

我们先考虑将如下方程组合并成一个方程
$$\{\begin{array}{c}ans\equiv x_1(\mathrm{m}\mathrm{o}\mathrm{d}{a}_1)\\ ans\equiv x_2(\mathrm{m}\mathrm{o}\mathrm{d}{a}_2)\end{array}$$

该方程存在唯一解ans，使得$ans=m_1{a}_1+x_1=m_2{a}_2+x_2$成立。
所以有方程
$$m_1{a}_1+x_1=m_2{a}_2+x_2$$即
$$m_1{a}_1-m_2{a}_2=x_2-x_1$$
成立。
如果有$gcd(a_1,a_2)|(x_2-x_1)$，我们则可以通过扩展欧几里得算法求得整数解。
否则无解。

在求得一组解(m₁,m₂)后，可以将上述方程组转化为以下方程：
$ans\equiv X(\mathrm{m}\mathrm{o}\mathrm{d}lcm(a_1,a_2))$
其中$X=m_1{a}_1+b_1=m_2{a}_2+b_2$。

显然X是上述方程组的一个特解。
设$M=lcm(a_1,a_2)$,那么通解为$X^{\prime }=kM+X$

可以将上述结论推广至由k个同余方程组成的一元线性同余方程组。

算法流程如下：
设x_k为前k个一元线性同余方程合并后的一个特解，x_k+1为第k+1个同余方程的特解。a_k为前k个一元线性同余方程合并后的模数，即$lcm(a1,a2,...,a_k)$。a_k+1为第k+1个同余方程的模数。
$fork=1ton:$
1、解方程$a_{k}x-a_{k+1}y=gcd(a_k,a_{k+1})$。
2、若$gcd(a_k,a_{k+1})|(x_{k+1}-x_k)$，则可以求出$a_{k}x-a_{k+1}y=x_{k+1}-x_k$的一个特解$x^{\prime }=\frac{x_{k+1}-x_k}{gcd(a_k,a_{k+1})}x$；否则无解，算法结束。
3、合并前k+1个方程的后一个特解为$x_{k+1}=a_k{x}^{\prime }+x_k$,
4、合并前k+1个方程的后的模数为$lcm(a_k,a_{k+1})$

知道对于方程$ax+by=gcd(a,b)$,由一对特解(x,y)，那么由通解X,Y，$X=k\cdot \frac{b}{gcd(a,b)}+x,Y=k\cdot \frac{a}{gcd(a,b)}+y$。
在上述算法流程中，如果$\frac{x_{k+1}-x_k}{gcd(a_k,a_{k+1})}$较大，直接与x相乘可能会造成溢出。我们可以把x’限制在$[0\frac{a_{k+1}}{gcd(a_k,a_{k+1})}]$中，通过慢速乘法(类似与快速幂)防止数据溢出。

模板代码:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
void exgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
	if(!b){
		x=1;y=0;d=a;return ;
	}
	else{
		exgcd(b,a%b,d,y,x);y-=(a/b)*x;
	}
}
ll mul(ll a,ll b,ll mod)
{
	ll ans=0;
	while(b)
	{
		if(b&1)ans=(ans+a)%mod;
		a=(a+a)%mod;
		b>>=1;
	}
	return ans;
}
int main(void)
{
	int n;
	scanf("%d",&n);
	ll b1=0,a1=1;
	ll b2,a2;
	while(n--)
	{
		scanf("%lld%lld",&a2,&b2);
		ll x,y,d,t;//求解x*a1+b1=b2(mod a2) 
		t=((b2-b1)%a2+a2)%a2;//x*a1=b2-b1(mod a2)
		exgcd(a1,a2,d,x,y);//求解x*a1+y*a2=gcd(a1,a2) 
		ll g=a2/d;
		x=mul(x,t/d,g);
		b1+=x*a1;
		a1*=a2/d;//lcm
		b1=(b1%a1+a1)%a1;
	}
	printf("%lld\n",b1);
}