CRT&EXCRT 中国剩余定理及其扩展

中国剩余定理详解

中国剩余定理(孙子定理)


有这样一个问题:
{x≡a1(mod m1)x≡a2(mod m2)……x≡an(mod mn)\begin{cases} x\equiv a_1(mod~m_1)\\ x\equiv a_2(mod~m_2)\\ ……\\ x\equiv a_n(mod~m_n) \end{cases}xa1(mod m1)xa2(mod m2)xan(mod mn)
其中,m两两互质,求满足所有同余方程的解。

求解:
  • 废话不多说,高中数学书上有讲过一个很对的算法。比如当前对于x≡a1(mod m1)x\equiv a_1(mod~m_1)xa1(mod m1)来说,我们找一个数t1t_1t1,让它是m2−mnm_2-m_nm2mn的公倍数,并且让它在模m1m_1m1的意义下等于1。然后我们把它乘上a1a_1a1,在模意义下它就等于a1a_1a1了。其他的方程都这样做。我们再把所有的ai∗tia_i*t_iaiti加起来,这个和就是答案(当然可能还需要取模)。
  • 为什么是正确的?我们依次看看每一个同余方程条件是否满足就好了。对于第iii个方程,除了我们的ai∗tia_i*t_iaiti,其他的a∗ta*tat都是mim_imi的倍数,所以同余方程显然成立,所以我们的算法正确性显然。
  • 因此,x=∑i=1nai∗ti% Mx=\sum_{i=1}^na_i*t_i\%~Mx=i=1naiti% M,其中MMM是所有数的lcmlcmlcm
校oj1898 曹冲称象
  • 模板题,但wa了居然半天没找到错误……
  • 数据只保证两两互质,但没说是质数,拿费马小定理求逆元就死的明明白白了。
  • exgcdexgcdexgcd求的xxx可能是负数,记得取模。
  • exgcdexgcdexgcd都不会写了,自闭了……
CodingCodingCoding
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=100;
int n,m[N],b[N];
ll M=1,ans;
ll exgcd(ll a,ll b,ll &x,ll &y){
	if(b==0){
		x=1,y=0;return a;
	}
	ll d=exgcd(b,a%b,x,y);
	ll tmp=x;x=y;y=tmp-a/b*x;
	return d;
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i){
		scanf("%d%d",&m[i],&b[i]);
		M*=m[i];
	}
	for(int i=1;i<=n;++i){
		ll Mi,t,x,y;
		Mi=M/m[i];
		exgcd(Mi,m[i],x,y);
		ans=(ans+1LL*x*Mi*b[i])%M;
	}
	cout<<(ans+M)%M<<endl;
	return 0;
}

扩展版

  • 还是上面的问题,这回不保证mmm两两互质了。我们先假设只有两个方程。用exgcdexgcdexgcd合并同余方程就好了。
CodingCodingCoding
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e5+10;
int n;ll m[N],a[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,x,y);
	ll tmp=x;x=y;y=tmp-a/b*y;
	return d;
}
ll mul(ll a,ll b,ll mod){
	ll res=0;
	for(;b;b>>=1){
		if(b&1) res=(res+a)%mod;
		a=(a+a)%mod;
	}
	return res;
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;++i) scanf("%lld%lld",&m[i],&a[i]);
	ll k1,k2,m1,m2,a1,a2,c,d,x,y,x0;
	m1=m[1],a1=a[1];
	for(int i=2;i<=n;++i){
		m2=m[i],a2=a[i];
		c=a2-a1;c=(c%m[i]+m[i])%m[i];
		d=exgcd(m1,m2,x,y);
		if(c%d){
			printf("-1\n");
			return 0;
		}
		ll lc=m1/d*m2;
		k1=(mul(x,c/d,m2)+m2)%m2;
		x0=mul(k1,m1,lc)+a1;x0=(x0+lc)%(lc);
		m1=lc;a1=x0;
	}
	cout<<x0<<endl;
	return 0;
}
[NOI2018]屠龙勇士
  • 每次攻击时攻击力是确定的。设攻击力为kkk,恢复力为ppp,生命值aaa,则有kx=py+akx=py+akx=py+a,类似于kx≡a(mod p)kx\equiv a(mod~p)kxa(mod p)。我们中国剩余定理合并同余方程的时候,xxx的系数必须是1,所以想办法消掉kkk。求逆元显然不合适,因为有没有逆元还不一定。
  • 可以求解一下这个同余方程啊。kx+py=a,x0=x′ad+kpdkx+py=a,x_0=x&#x27;\frac{a}{d}+k\frac{p}{d}kx+py=a,x0=xda+kdp,其中x0x_0x0是方程的通解。那么可以转化为同余方程x0≡x′ad(mod pd)x_0\equiv x&#x27;\frac{a}{d}(mod~\frac{p}{d})x0xda(mod dp)。此时就可以合并了。
  • 细节问题1:最终的xxx一定要保证把龙杀死。也就是x&gt;maxaatkx&gt;max{\frac{a}{atk}}x>maxatka,注意向上取整。
  • 细节问题2:一个是如果生命值和攻击力都是恢复力的倍数,显然你攻击次数任意,只要满足xxx把龙杀死就好,不参与合并。
  • 细节问题3:如果生命值不是恢复力的倍数,而攻击力是,那么一定无解,直接输出-1。
  • 细节问题4:其实也不是什么细节,就是自己有点sb。exgcdexgcdexgcd求方程的解的时候,a,ba,bab的顺序和对应的x,yx,yxy的关系一定要搞好,否则错的你不明不白的。
CodingCodingCoding
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e5+10;
ll t,n,m,atk[N],p[N],a[N];
multiset<ll>s;
multiset<ll>::iterator it;
ll in(){
	char ch=getchar();ll num=0,f=1;
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){num=num*10+(ch^48);ch=getchar();}
	return num*f;
}
ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1;y=0;return a;}
	ll d=exgcd(b,a%b,x,y);
	ll tmp=x;x=y;y=tmp-a/b*x;
	return d;
}
ll mul(ll a,ll b,ll mod){
	ll res=0;
	for(;b;b>>=1){
		if(b&1) res=(res+a)%mod;
		a=(a+a)%mod;
	}
	return res;
}
int main(){
	t=in();
	E:while(t--){
		n=in(),m=in();ll temp;
		for(int i=1;i<=n;++i) a[i]=in();
		for(int i=1;i<=n;++i) p[i]=in();
		for(int i=1;i<=n;++i) atk[i]=in();
		s.clear();
		for(int i=1;i<=m;++i){temp=in();s.insert(temp);}
		ll k,d,tmp,lcm,x,y,c1=0,m1=1,c,m,M=0;
		for(int i=1;i<=n;++i){
			it=s.upper_bound(a[i]);
			if(it!=s.begin()) it--;
			k=*it,s.erase(it);
			s.insert(atk[i]);
			M=max(M,(a[i]-1)/k+1);
			k%=p[i];a[i]%=p[i];
			if(!k&&!a[i]) continue;
			if(!k&&a[i]){printf("-1\n");goto E;}
			d=exgcd(k,p[i],x,y);
			if(a[i]%d){printf("-1\n");goto E;}
			m=p[i]/d;c=(mul((x%m+m)%m,a[i]/d,m)+m)%m;
			d=exgcd(m1,m,x,y);tmp=c-c1;lcm=m1/d*m;tmp=(tmp+lcm)%lcm;
			if(tmp%d){printf("-1\n");goto E;}
			x=(mul(x,tmp/d,m/d)+m/d)%(m/d);
			c1=(mul(m1,x,lcm)+c1)%lcm;m1=lcm;
			//m1=m1/d*m;
			//c1=(c1+mul(mul(m1/m,((c-c1)%m1+m1)%m1,m1),(x%m1+m1)%m1,m1))%m1;
		}
		printf("%lld\n",c1>=M?c1:c1+m1*((M-c1-1)/m1+1));
	}
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值