Loading

扩展中国剩余定理(EXCRT)

P4777 【模板】扩展中国剩余定理(EXCRT)

老师说NOIP会考,那还是复习&总结一下吧

求方程的最小正整数解:

\[\begin{cases} x\equiv a_1 \pmod {b_1}\\ x\equiv a_2 \pmod {b_2}\\ \cdots\\ x\equiv a_n \pmod {b_n}\\ \end{cases} \]

假设现在已经求出了前 \(i-1\) 个方程的最小正整数解 \(res\) ,记 \(lcm=\operatorname{lcm}(b_1,b_2,\cdots,b_{i-1})\)

当加入第 \(i\) 个方程的时候,我们需要求出一个整数 \(x\) 满足 \(res+x\times lcm \equiv a_i \pmod {b_i}\)

移项可得 \(x\times lcm + y \times b_i = a_i - res\)

这是一个扩展欧几里得的形式,方程组无解当且仅当存在一个扩展欧几里得方程无解,即 \((a_i - res) \% \gcd(lcm, b_i) \not = 0\)

否则根据小奥知识,\(res' = res + x \times lcm \% lcm'\)\(lcm'\) 是指加入 \(b_i\) 之后的 \(lcm\)

但是直接拿 \(x\) 去乘 \(lcm\) 很可能爆 long long

我们其实可以把 \(x\) 通过取模限制在 \(\dfrac{b_i}{\gcd(b_i, lcm)}\) 以内的,这样只要出题人保证所有数的 \(lcm\) 最后不爆 long long 就行了

具体方法如下:

首先先看一个化简版的(变量与上面无关)

\[xy \% yz = y \times (x \% z) \]

这个其实很好理解,就是余式两边同时乘上 \(y\) 结果不变

那么回到本题

\[x \times lcm \% lcm' = x \times lcm \% \dfrac{lcm \times b_i}{\gcd(lcm, b_i)}=lcm \times (x \% \dfrac{b_i}{\gcd(lcm, b_i)}) \]

直接取模就行了

当然在实现的时候,由 \(exgcd\) 特性我们会先解出 \(x_0,y_0\) 满足 \(x_0 \times lcm + y_0 \times b_i = \gcd(lcm, b_i)\) ,我们的 \(x\) 应该等于 \(x_0 \times \dfrac{a_i - res}{\gcd(lcm, b_i)}\) ,这一步要龟速乘或者光速乘,可以在这一步取模。

但是由于 \(exgcd\) 的缘故 \(EXCRT\) 必然带 \(\log\) ,用光速乘意义不大。

另外,这东西严格吊打 CRT 啊,到现在为止还没见过 CRT 能做 EXCRT 不能做的事。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
#define x first
#define y second
#define sz(v) (int)v.size()
#define pb(x) push_back(x)
#define mkp(x,y) make_pair(x,y)
#define int long long//------------------Warning!!!-------------------------------------------
//#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
//char buf[1<<21],*p1=buf,*p2=buf;
inline int read(){
	int x=0,f=1;char c=getchar();
	while(!isdigit(c)){if(c=='-')f=0;c=getchar();}
	while(isdigit(c))x=x*10+c-'0',c=getchar();
	return f?x:-x;
}
#define N 100005
int n,a[N],b[N];
void exgcd(int&x,int&y,int&d,int a,int b){
	if(!b)return x=1,y=0,d=a,void();
	exgcd(y,x,d,b,a%b),y-=a/b*x;
}
int mul(int x,int y,int mod){
//	int res=x*y-(long long)((long double)x/mod*y+0.5)*mod;
//	return res<0?res+mod:res;//光速乘
	int res=0;
	while(y){
		if(y&1)(res+=x)%=mod;
		y>>=1,(x<<=1)%=mod;
	}
	return res;
}//龟速乘
int excrt(){
	int res=a[1],lcm=b[1];
	for(int i=2;i<=n;++i){
		int x,y,d,t=((a[i]-res)%b[i]+b[i])%b[i];
		exgcd(x,y,d,lcm,b[i]);
		if(t%d)return -1;
		res+=lcm*mul(x,t/d,b[i]/d);
		lcm=lcm/d*b[i];
		res=(res+lcm)%lcm;
	}
	return res;
}
signed main(){
	n=read();
	for(int i=1;i<=n;++i)b[i]=read(),a[i]=read();
	printf("%lld\n",excrt());
	return 0;
}
posted @ 2020-11-01 09:11  zzctommy  阅读(88)  评论(0编辑  收藏  举报