BZOJ 3601 一个人的数论(容斥+高斯消元/拉格朗日插值)

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统计答案
}