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

传送门

---

中国剩余定理(CRT)

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

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

怎么做呢？
可以仿照拉格朗日插值法，构造一组通解：

\[ans=\sum_{i=1}^{n}b_i\prod_{j=1,j\ne i}^{n}(a[j]*a[j]^{-1}) \]

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

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

所以

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

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

总结一下思路

预处理出 \(\prod a_i\) 和 \(\psi(1)\text{ 到 }\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;
}