中国剩余定理 学习笔记

目录

• CRT
• exCRT

中国剩余定理 (Chinese Remainder Theorem) 简称CRT。

CRT

CRT可以求解一个这样的方程组的最小解

$$\left\{ \begin{array}{**lr**} x \equiv b_1 \pmod {a_1} \\ x \equiv b_2 \pmod {a_2} \\ x \equiv b_3 \pmod {a_3} \\ \vdots \\ x \equiv b_k \pmod {a_k} \\ \end{array} \right.$$

这里保证数组 $a$ 两两互质。

CRT的原理就是得到 $c_i$ 使 $c_i \equiv 1 \pmod {a_i}$ 并且 $c_i \equiv 0 \pmod {a_j},i\not= j$
显然答案就是 $\sum^k_{i=0} b_ic_i\bmod \left(\prod^{n}_{i=1}a_i\right)$

算法流程如下：

1. 计算 $a=\prod a_i$
2. 对于第 $i$ 个方程：
   a. 计算 $m_i=\frac{a}{a_i}$
   b. 计算 $m_i$ 在模 $a_i$ 的逆元 $m_i^{-1}$
   c. 计算 $c_i=m_i\times m_i^{-1}$
3. 方程的解为 $x \equiv \sum^k_{i=1}b_ic_i \pmod a$

证明略 （显然啊）

代码(洛谷P1495)：

#include<cstdio>
#define maxn 100039
using namespace std;
//#define debug
typedef long long ll;
typedef long long Type;
inline Type read(){
	Type sum=0;
	int flag=0;
	char c=getchar();
	while((c<'0'||c>'9')&&c!='-') c=getchar();
	if(c=='-') c=getchar(),flag=1;
	while('0'<=c&&c<='9'){
		sum=(sum<<1)+(sum<<3)+(c^48);
		c=getchar();
	}
	if(flag) return -sum;
	return sum;
}
ll n,a[maxn],b[maxn];
ll sn,m,c,ans;
ll gcd(ll a,ll b){
	if(a%b==0) return b;
	return gcd(b,a%b);
}
inline ll lcm(ll a,ll b){ return a*b/gcd(a,b); }
void exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){ x=1; y=0; return; }
	exgcd(b,a%b,x,y);
	ll t=x;
	x=y; y=t-a/b*y;
	return;
}
ll js(ll a1,ll b1,ll a2,ll b2){
	ll p,q;
	exgcd(a1,-a2,p,q);
	if((b2-b1)%gcd(a1,-a2)) return -1;
	q*=(b2-b1)/gcd(a1,-a2);
	p*=(b2-b1)/gcd(a1,-a2);
	//if(a1*p+b1 != a2*q+b2) printf("Oops!\n");
	//else printf("OK\n");
	int r=lcm(a1,a2);
	return ((a2*q+b2)%r+r)%r;
}//x=b_i(mod a_i)
ll excrt(){
	ll ta,tb;
	ta=a[1],tb=b[1];
	for(int i=2;i<=n;i++){
		tb=js(ta,tb,a[i],b[i]);
		if(tb==-1) return -1;
		ta=lcm(ta,a[i]);
	}
	return tb;
}
int main(){
    //freopen("1.in","r",stdin);
    //freopen("my.out","w",stdout);
    n=read(); sn=1;
    for(int i=1;i<=n;i++)
    	a[i]=read(),b[i]=read();//x=b_i(mod a_i)
	printf("%lld",excrt());
	return 0;
}

exCRT

如果数组 $a$ 不是两两互质呢？
我们尝试修复一下CRT的方法让它变成exCRT，但是我们发现，CRT的核心都说过了，它是求出 $c_i$ 使 $c_i \equiv 1 \pmod {a_i}$ 并且 $c_i \equiv 0 \pmod {a_j},i\not= j$ ，但是我们发现，我们有时不能找到这样的 $c_i$ ，所以CRT原来的思路是不能实现的。
看来只能用另外一种方法了。
考虑两个方程的情况。

$$\left\{ \begin{array}{**lr**} x \equiv b_1 \pmod {a_1} \\ x \equiv b_2 \pmod {a_2} \\ \end{array} \right.$$

令 $x=a_1q+b_1=a2_p+b_2$
移项，得到 $a_1q+a_2(-p)=b_2-b_1$
显然只要exgcd解出 $p,q$ 的一组解就可以解决问题了，然后显然合并后的方程是：

$$x \equiv a_1q+b_1 \pmod {\operatorname{lcm}\left(a_1,a_2\right)}$$

然后就可以过 洛谷P4777了，当然要注意精度问题，加上防爆乘，不然会炸精度。。。

#include<cstdio>
#define maxn 100039
using namespace std;
//#define debug
typedef long long ll;
typedef long long Type;
inline Type read(){
	Type sum=0;
	int flag=0;
	char c=getchar();
	while((c<'0'||c>'9')&&c!='-') c=getchar();
	if(c=='-') c=getchar(),flag=1;
	while('0'<=c&&c<='9'){
		sum=(sum<<1)+(sum<<3)+(c^48);
		c=getchar();
	}
	if(flag) return -sum;
	return sum;
}
ll mol(ll a,ll b,ll p){
	int flag=0;
	if(a<0) a=-a,flag^=1;
	if(b<0) b=-b,flag^=1;
	ll ans=0,tmp=a;
	while(b){
		if(b&1) ans=(ans+tmp)%p;
		b>>=1; tmp=(tmp<<1)%p;
	}
	if(flag) return (-ans%p+p)%p;
	return ans%p;
}
ll gcd(ll a,ll b){
	if(a%b==0) return b;
	return gcd(b,a%b);
}
inline ll lcm(ll a,ll b){ return a/gcd(a,b)*b; }
void exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){ x=1; y=0; return; }
	exgcd(b,a%b,x,y);
	ll t=x;
	x=y; y=t-(a/b)*y;
	return;
}
int n;
ll a[maxn],b[maxn];
ll excrt(){
	ll m=a[1],ans=b[1];
	ll p,q;
	for(int i=2;i<=n;i++){
		exgcd(m,a[i],p,q);
		//if((b[i]-ans)%gcd(m,a[i])!=0) return -1;
		p=mol(p,(b[i]-ans)/gcd(m,a[i]),a[i]/gcd(m,a[i]));
		ans=(mol(p,m,lcm(a[i],m))+ans)%lcm(a[i],m);
		m=lcm(a[i],m);
	}
	return ans%m;
}//x=b(mod a)
int main(){
	//freopen("P4777_13.in","r",stdin);
	//freopen(".in","w",stdout);
	n=read();
	for(int i=1;i<=n;i++)
		a[i]=read(),b[i]=read();
	printf("%lld",excrt());
	return 0;
}