BZOJ 3601 一个人的数论(容斥+高斯消元/拉格朗日插值)
题意：求\(\sum _1^{n}[gcd(i,n)=1]i^d\)
其中\(n=\Pi_1^mp_i^{c_i}\)，\(m\leq 1000,p_i,c_i\leq 10^9\)
???丧心病狂???
我们采取容斥来计算，但是首先要把这个\(i^d\)求和搞定
由于神奇的数学归纳法，我们知道，\(F(n)=\sum_{i=1}^ni^d\)可以用一个\(d+1\)次的多项式表示
即$F(n)=\sum _{i=1}^{d+1}w_i \cdot n^i $
带入\(d+1\)个值，求出\(w_i\)，当然我们可以高斯消元/拉格朗日插值
// 暴力高斯消元
namespace Gauss{
ll a[111][111];
void Solve() {
ll tmp=0;
rep(i,0,d) {
rep(j,0,d) a[i][j]=qpow(i+1,j+1);
tmp=(tmp+qpow(i+1,d))%P;
a[i][d+1]=tmp; // 预处理点值式
}
rep(i,0,d) {
rep(j,i,d) if(a[j][i]) {
swap(a[i],a[j]);
break;
}
if(!a[i][i]) continue;
ll Inv=qpow(a[i][i],P-2);
rep(j,i+1,d) {
ll t=a[j][i]*Inv%P;
rep(k,i,d+1) a[j][k]=(a[j][k]-a[i][k]*t%P+P)%P;
}
}
drep(i,d,0) {
w[i+1]=a[i][d+1]*qpow(a[i][i],P-2)%P;
rep(j,0,i-1) a[j][d+1]=(a[j][d+1]-a[j][i]*w[i+1]%P+P)%P;
}
}
}
考虑\(2^m\)枚举容斥，考虑枚举当前的数包含的因子集合为\(S\)，答案就是
\[\sum _{S\in p} (-1)^{|S|}\sum_{i=1}^{\frac{n}{\Pi S_j}} (i \Pi S_j)^d
\]
\[=\sum _{S\in p} (-1)^{|S|}(\Pi S_j)^d F(\frac{n}{\Pi S_j})
\]
n=1,d=rd(),m=rd();
rep(i,1,m) {
a[i]=rd(),b[i]=rd();
n=1ll*n*qpow(a[i],b[i])%P;
}
Gauss::Solve();
rep(i,0,(1<<m)-1) {
ll t=1,res=0,c=0;
rep(j,1,m) if(i&(1<<(j-1))) t=t*a[j]%P,c^=1;
ll x=n*qpow(t,P-2)%P;
rep(j,1,d+1) res=(res+qpow(x,j)*w[j])%P;
res=res*qpow(t,d)%P;
if(c) ans=(ans-res)%P;
else ans=(ans+res)%P;
}
ans=(ans%P+P)%P;
printf("%lld\n",ans),ans=0;
考虑\(F(n)\)的每一项\(w_in^i\)对于答案的贡献
\[\sum _{S\in p} (-1)^{|S|}(\Pi S_j)^d w_i(\frac{n}{\Pi S_j})^i
\]
这个值可以压进转移中直接完成计算
rep(i,1,d+1) {
dp[0]=1;
rep(j,1,m) {
dp[j]=dp[j-1]*qpow(qpow(a[j],b[j]),i)%P; //不选这个因数
dp[j]=(dp[j]-dp[j-1]*qpow(qpow(a[j],b[j]-1),i)%P*qpow(a[j],d))%P;// 选这个因数,系数取反
}
ans=(ans+dp[m]*w[i])%P; // 对于每个wi统计答案
}