EXCRT

求一个形如(xequiv a_i (mod p_i)) 的线性方程组的解

我们考虑合并两个方程:

[egin{cases} x=a_1 (mod p_1)\ \ x=a_2 (mod p_2)\ end{cases} ]

先考虑将(x)设为(a_1),再不断的加(p_1)使得(x)在模(p_2)意义下为(a_2)

我们可以用(exgcd)求得一组(mp_1+np_2= gcd(p_1,p_2))的解,即(mp_1equiv gcd(p1,p2) (mod p2))

也就是说我们每次增加(mp_1)即可让(x)在模(p_1)意义下增加(gcd(p_1,p_2))

(d=gcd(p_1,p_2)),那么我们应该增加((a_2-a_1)/d) 次,即(x)应该加上((a_2-a_1)/d*mp_1),然后我们知道方程合并后的模数应该是(P=p_1*p_2/gcd(p_1,p_2)),所以再将结果对(P)取模即可。

代码:

LL exgcd(LL x,LL y,LL &a,LL &b){
	LL z;
	return y?(z=exgcd(y,x%y,b,a),b-=x/y*a,z):(a=1,b=0,x);
}
LL mult(LL a,LL b,LL p){
	LL ans=0;
	while(b){
		if(b&1)ans=(ans+a)%p;
		b>>=1; a=(a<<1)%p;
	}
	return ans;
}
struct CRT{
	LL a[N],p[N],P,x;
	bool excrt(LL a1,LL p1,LL a2,LL p2){
		LL m,n,d=exgcd(p1,p2,m,n); 
		LL a=(a2-a1%p2+p2)%p2;
		if(a%d!=0)return false;
		P=p1/d*p2; m=(m%p2+p2)%p2;
		x=(x+mult(a/d*p1,m,P))%P; //注意这里一般先计算a/d*p1,然后再乘m,因为题目一般保证了所有p_i的最小公倍数小于某个值,这样前面的运算就不会溢出,而后面的乘m可能会溢出,所以还要用快速乘。
		return true;
	}
	bool run(){
		x=a[1],P=p[1];
		for(int i=2;i<=n;i++)
			if(!excrt(x,P,a[i],p[i]))
				return false;
		return true;
	}
}crt;




[NOI2018]屠龙勇士

题目要把一个形如(Axequiv a (mod p))的方程变成(xequiv b(mod p))的形式。

肯定有无解的情况,我们求(d=gcd(A,p)),然后将方程每一项除以(d),如果(d)不是(a)的约数,那么方程无解。

然后我们得到了等价方程(A'xequiv a' (mod p')),其中(A'=A/d, a'=a/d, p'=p/d)

我们用(exgcd)求得(mA+np=gcd(A,p))的一组解,则这组解满足(mA'+np'=1),即(mA'equiv1 (mod p'))(m)即为(A')在模(p')意义下的逆元。

这样上面的方程就可以转化为(xequiv a'm (mod p'))

这样就是一个(xequiv b(mod p))的形式了。

代码:

bool flag=0;
for(int i=1;i<=n;i++){
	LL A1=A[i],a1=a[i],p1=p[i];
	LL m,n,d=exgcd(A1,p1,m,n); 
	if(a1%d!=0){flag=1; break;}
	A1/=d,a1/=d,p1/=d; m=(m%p1+p1)%p1;
	crt.a[i]=mult(a1,m,p1),crt.p[i]=p1;
}

完整代码

原文地址:https://www.cnblogs.com/lishuyu2003/p/13287249.html