扩展中国剩余定理(EXCRT)学习笔记

EXCRT 其实和 CRT 半毛钱关系都没有……

我们在 CRT 中需要用到求逆元的操作：

$\operatorname{inv}(\dfrac{M}{b_i},b_i)$

但是如果 $b_i$ 不满足两两互质，我们就无法求出逆元，这样 CRT 就用不了了。

为了解决这个问题，求出 $b_i$ 不一定两两互质时同余方程组

$\begin{cases}x\equiv a_1\pmod{b_1}\\x\equiv a_2\pmod{b_2}\\x\equiv a_3\pmod{b_3}\\\qquad\dots\dots\\x\equiv a_n\pmod{b_n}\end{cases}$

的解，EXCRT 就诞生了。

首先考虑两个方程的情况：

$\begin{cases}x\equiv a_1\pmod{b_1}\\x\equiv a_2\pmod{b_2}\end{cases}$

可以转化为这样：

$\begin{cases}x=a_1+k_1\cdot b_1\\x=a_2+k_2\cdot b_2\end{cases}$

即：

$a_1+k_1\cdot b_1=a_2+k_2\cdot b_2$

移项：

$k_1\cdot b_1-k_2\cdot b_2=a_2-a_1$

这个柿子很明显能使用 exgcd 求解。

很明显，若 $a_2-a_1\mod{\gcd(b_1,b_2)} \ne 0$ 那么无解。

设 $y_1,y_2$ 满足 $y_1\cdot b_1+y_2\cdot b_2=\gcd(b_1,b_2)$，则 $k_1=y1\cdot \dfrac{a_2-a_1}{\gcd(b_1,b_2)},k_2=-y_2\cdot \dfrac{a_2-a_1}{\gcd(b_1,b_2)}$。

所以我们可以求出这个方程组的一个解为 $x=a_1+y1\cdot \dfrac{a_2-a_1}{\gcd(b_1,b_2)}\cdot b_1$。

很明显最小正整数解就是 $x_{min}=x\mod{\operatorname{lcm}(b_1,b_2)}$，但由于 $x$ 有可能是负数，所以取模时需要特殊处理。

接下来 EXCRT 的核心思想来了：

我们可以把 $x=a_1+y1\cdot \dfrac{a_2-a_1}{\gcd(b_1,b_2)}\cdot b_1\mod{\operatorname{lcm}(b_1,b_2)}$ 改写成同余方程，即：

$x\equiv a_1+y1\cdot \dfrac{a_2-a_1}{\gcd(b_1,b_2)}\cdot b_1\mod{\operatorname{lcm}(b_1,b_2)}\pmod{\operatorname{lcm}(b_1,b_2)}$

这时，我们就成功把两个同余方程合并为了一个！

所以只要一直合并下去，就能得到最终的解了。

模板题代码：（需要高精所以开了 __int128）

#include <iostream>
#include <cstdio>

using namespace std;

int n;
__int128 a[100005],b[100005];

inline __int128 read()
{
	__int128 s=0;
	int w=1,ch=getchar();
	while(ch<'0'||ch>'9') ch=='-'?w=-1,ch=getchar():ch=getchar();
	while(ch>='0'&&ch<='9') s=(s<<1)+(s<<3)+(ch^48),ch=getchar();
	return s;
}

void writen(__int128 x)
{
	if(x<0) putchar('-'),writen(-x);
	x>9?writen(x/10),putchar(x%10|48):putchar(x|48);
}

inline __int128 lcm(__int128 a,__int128 b)
{
	__int128 tmpa=a,tmpb=b;
	__int128 t=a%b;
	while(t!=0)
	{
		a=b;
		b=t;
		t=a%b;
	}
	return tmpa*tmpb/b;
}

__int128 exgcd(__int128 a,__int128 b,__int128 &x,__int128 &y)
{
	if(a<b)
	{
		return exgcd(b,a,y,x);
	}
	if(b==0)
	{
		x=1;
		y=0;
		return a;
	}
	__int128 tmpx,tmpy;
	__int128 res=exgcd(b,a%b,tmpx,tmpy);
	x=tmpy;
	y=tmpx-a/b*tmpy;
	return res;
}

int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		b[i]=read();
		a[i]=read();
	}
	a[1]=(a[1]%b[1]+b[1])%b[1];
	for(int i=2;i<=n;i++)
	{
		__int128 rig=a[i]-a[i-1];
		__int128 y1,y2;
		__int128 g=exgcd(b[i-1],b[i],y1,y2);
		if(rig%g!=0)
		{
			puts("Impossible!");
			break;
		}
		y1*=rig/g;
		a[i]=a[i-1]+b[i-1]*y1;
		b[i]=lcm(b[i-1],b[i]);
		a[i]=(a[i]%b[i]+b[i])%b[i];
	}
	writen(a[n]);
	putchar('\n');
	return 0;
}