中国剩余定理 学习笔记

Posted 2021-05-20 jiangtaizhe001

目录

• 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\\)

算法流程如下：

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;
}