P1495 【模板】中国剩余定理(CRT)/曹冲养猪

传送门


中国剩余定理(CRT)

我的第一反应是小学奥数题——韩信点兵。
转化成数学语言,就是给你 (n) 个关于 (x) 的同余方程(保证 (a_i) 互质),求最小整数解。

[egin{cases} xequiv b_1pmod {a_1}\ xequiv b_2pmod {a_2}\ xequiv b_3pmod {a_3}\ cdotscdotscdotscdots\ end{cases}]

怎么做呢?
可以仿照拉格朗日插值法,构造一组通解:

[ans=sum_{i=1}^{n}b_iprod_{j=1,j e i}^{n}(a[j]*a[j]^{-1}) ]

其中 (a[j]^{-1}) 表示的是 (a[j]) 在模 (a[i]) 下的逆元。
因为 (a[i]) 之间两两互质,所以我们可以放心的用欧拉定理

[a^{psi(m)}equiv 1pmod m ]

所以

[a^{psi(m)-1}equiv frac{1}{a}pmod m ]

可见 (a[j]^{psi(a[i])-1}) 即是逆元。

总结一下思路

预处理出 (prod a_i)(psi(1) ext{ 到 }psi(max\_a)),带入公式中,利用快速幂求出逆元,把答案累加起来。

注意事项

  • 数据有锅,第一个点的 (a) 高达 (1093258),注意数组开的范围
  • 模数最大为 (10^{18}) ,所以需要使用快(gui)速乘

AC代码

#include<iostream>
using namespace std;
bool vis[2000005];
int n,cnt,prime[2000005];
long long a[15],b[15],prod_a=1,phi[2000005],maxa,ans;
inline long long ksc(long long x,long long y,long long mod)
{
    return (x*y-(long long)((long double)x/mod*y)*mod+mod)%mod;     
}
long long q_p(long long a,long long b){
	if(b==0) return 1;
	if(b==1) return a%prod_a;
	long long res=q_p(a,b/2);
	if(b&1) return ksc(ksc(res,res,prod_a),a,prod_a);
	return ksc(res,res,prod_a);
}
int main()
{
    cin>>n;
    for(int i=1;i<=n;i++) cin>>a[i]>>b[i];
    for(int i=1;i<=n;i++) prod_a*=a[i],maxa=max(maxa,a[i]);
    phi[1]=1;
    for(int i=2;i<=maxa;i++){
    	if(!vis[i]) prime[++cnt]=i,phi[i]=i-1;
    	for(int j=1;j<=cnt&&i*prime[j]<=maxa;j++){
    		vis[i*prime[j]]=1;
    		if(i%prime[j]){
    			phi[i*prime[j]]=phi[i]*phi[prime[j]];
			}else{
				phi[i*prime[j]]=phi[i]*prime[j];
				break;
			}
		}
	}
    for(int i=1;i<=n;i++){
    	ans=(ans+ksc(ksc(prod_a/a[i],q_p(prod_a/a[i],phi[a[i]]-1),prod_a),b[i],prod_a))%prod_a;
	} 
	cout<<ans;
    return 0;
}
原文地址:https://www.cnblogs.com/yinyuqin/p/14869428.html