【BZOJ3601】一个人的数论

题目链接

题意简述

求小于 n 且与 n 互质的数的 k 次方之和。

Sol

要求的东西:

\[\sum_{i=1}^n i^k [gcd(i,n)=1] \]

枚举 gcd 上个莫比乌斯函数:

\[\sum_{i=1}^n i^k \sum_{d|n,d|i} \mu(d) \]

交换求和顺序

\[\sum_{d|n} \mu(d) \sum_{i=1}^{\frac{n}{d}} (id)^k \]

这样就差不多没有办法继续了:

\[\sum_{d|n} \mu(d)*d^k \sum_{i=1}^{\frac{n}{d}} i^k \]

题目中给出的是 n 的分解式 , 这使得我们往积性函数上想。

发现前面的 \(\mu(d)\) 是积性函数 ，\(d^k\) 是个完全积性函数，那么就只剩下最后那个东西了。

这就是个自然质数幂嘛，它是个关于 \(x=\frac{n}{d}\) 的 \(k+1\)次多项式，这样子我们可以把这个多项式的系数直接高斯消元或着拉格朗日插值给弄出来，假设系数是 \(a_i\) ，那么就是:

\[\sum_{d|n} \mu(d)*d^k \sum_{i=0}^{k+1} a_i\bigg(\frac{n}{d}\bigg)^i \]

交换个求和顺序就好了:

\[\sum_{i=0}^{k+1}a_i\sum_{d|n} \mu(d)*d^k \bigg(\frac{n}{d}\bigg)^i \]

后面就是个狄利克雷卷积，它也是个积性函数了，直接计算好 \(n\) 是 \(p^k\) 时的结果乘起来就好了，因为莫比乌斯函数必须是在无平方因子时才不为 \(0\)，所以 \(d\) 只有 \(1\) 和 \(p\) 两种取值，直接算就行了。

code:

#include<bits/stdc++.h>
using namespace std;
template<class T>inline void init(T&x){
	x=0;char ch=getchar();bool t=0;
	for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
	if(t) x=-x;return;
}
typedef long long ll;
const int mod=1e9+7,phi=mod-1;
template<class T>inline void Inc(T&x,int y){x+=y;if(x>=mod) x-=mod;return;}
template<class T>inline void Dec(T&x,int y){x-=y;if(x < 0 ) x+=mod;return;}
template<class T>inline int fpow(int x,T k){int ret=1;for(;k;k>>=1,x=(ll)x*x%mod)if(k&1) ret=(ll)ret*x%mod;return ret;}
inline int Sum(int x,int y){x+=y;return x>=mod? (x-mod):x;}
inline int Dif(int x,int y){x-=y;return x < 0 ? (x+mod):x;}
const int N=1021;
int d,w;
namespace Lagerange{
	int coef[N],Y[N];
	typedef vector<int> Poly;
	inline Poly Mul(Poly A,Poly B){
		Poly Res;int szl=A.size(),szr=B.size();Res.resize(szl+szr-1);
		for(int i=0;i<szl;++i) for(int j=0;j<szr;++j) Inc(Res[i+j],(ll)A[i]*B[j]%mod);
		return Res;
	}
	Poly Divide(int l,int r,int I){
		if(l==r) {
			int len=r-l+1;Poly Res;Res.resize(len+1);
			if(l==I) Res[0]=1;
			else {int inv=fpow(Dif(I,l),mod-2);Res[0]=(ll)(mod-l)*inv%mod;Res[1]=inv;}
			return Res;
		}
		int mid=l+r>>1;Poly A=Divide(l,mid,I),B=Divide(mid+1,r,I);
		return Mul(A,B);
	}
	inline void Solve(int K){
		int LIM=K+2;Y[0]=0;
		for(int i=1;i<=LIM;++i) Y[i]=Sum(Y[i-1],fpow(i,K));
		for(int i=1;i<=LIM;++i) {
			Poly ret=Divide(1,LIM,i);
			int sz=ret.size();
			for(int j=0;j<sz;++j) Inc(coef[j],(ll)Y[i]*ret[j]%mod);
		}
		return;
	}
}using Lagerange::coef;
int P[N],alpha[N];
inline int Calc(int i){
	int ret=1;
	for(int j=1;j<=w;++j) ret=(ll)ret*Dif(fpow(P[j],(ll)alpha[j]*i%phi),fpow(P[j],((ll)(alpha[j]-1)*i%phi+d)%phi))%mod;
	return ret;
}
int main()
{
	init(d),init(w);
	Lagerange::Solve(d);
	int ans=0;
	for(int i=1;i<=w;++i) init(P[i]),init(alpha[i]);
	for(int i=0;i<=d+1;++i) Inc(ans,(ll)coef[i]*Calc(i)%mod);
	cout<<ans<<endl;
}