PE639 Summing a multiplicative function

有积性函数(f(i))，其中(f_k(p^e)=p^k,e>0)。

(S_k(n)=sum_{i=1}^nf_k(i))

求(sum_{k=1}^{50} S_k(10^{12}))

---

学习狄利克雷生成函数的时候碰见的例题：https://zhuanlan.zhihu.com/p/50817119

很容易写出生成函数：

[F(x)=prod_p(1+frac{p^k}{p^x}+frac{p^k}{p^{2x}}+dots) ]

推一下式子：

[F(x)=prod(1+frac{p^k}{p^x}frac{1}{1-p^{-x}})\=prod(frac{1+p^{-x}+p^{k-x}}{1-p^{-x}})\=prod(frac{1+p^{-x}+p^{k-x}}{1-p^{-x}}frac{1-p^{k-x}}{1-p^{k-x}})\=(sumfrac{1}{1-p^{k-x}})prod(1+(p^k-p^{2k})(p^{-2x}+p^{-3x}+dots))\=zeta(x-k)prod(1+(p^k-p^{2k})(p^{-2x}+p^{-3x}+dots))\=(sumfrac{n^k}{n^x})prod(1+(p^k-p^{2k})(p^{-2x}+p^{-3x}+dots)) ]

可以发现右半边的式子中只有powerful number处有值。于是枚举右半边式子中有值的位置，它的贡献就是系数乘上左半边对应的前缀和。powerful number有(O(sqrt n))个。前缀和可以(O(k))算。

所以理论上时间复杂度是(O(k^2sqrt n))，用家里超慢的机子跑了108s。

---

using namespace std;
#include <bits/stdc++.h>
#define ll long long
#define mo 1000000007
#define K 55
#define M 1000005
ll qpow(ll x,ll y=mo-2){
	ll r=1;
	for (;y;y>>=1,x=x*x%mo)
		if (y&1)
			r=r*x%mo;
	return r;
}
int p[M],np,c[M];
bool inp[M];
int mu[M];
void initp(int n){
	mu[1]=1;
	for (int i=2;i<=n;++i){
		if (!inp[i])
			p[++np]=i,mu[i]=-1;
		for (int j=1;j<=np && i*p[j]<=n;++j){
			inp[i*p[j]]=1;
			if (i%p[j]==0)
				break;
			mu[i*p[j]]=-mu[i];
		}
	}
}
ll n,m;
int maxk,k;
ll S[K][K],inv[K],sc[K];
int presum(ll n){
	int s=0,r=(n+1)%mo,t=r;
	for (int j=1;j<=k;++j){
		--r;
		t=(ll)t*r%mo;
		s=(s+sc[j]*t)%mo;
//		s=(s+S[k][j]*t%mo*inv[j+1])%mo;
	}
	return s;
}
ll ans;
void dfs(int i,ll v,ll w){
	(ans+=presum(n/v)*w)%=mo;
	for (;i<=np && (__int128)v*p[i]*p[i]<=n;++i){
		ll v_=v*p[i],w_=w*c[i]%mo;
		while ((__int128)v_*p[i]<=n)
			dfs(i+1,v_*=p[i],w_);
	}
}
ll doit(int _k){
	k=_k;
	for (int j=1;j<=k;++j)
		sc[j]=S[k][j]*inv[j+1]%mo;
	for (int i=1;i<=np;++i){
		ll tmp=qpow(p[i],k);
		c[i]=(tmp-tmp*tmp%mo+mo)%mo;
	}
	dfs(1,1,1);
}
int main(){
	freopen("in.txt","r",stdin);
	scanf("%lld%d",&n,&maxk);
	m=sqrt(n);
	while ((m+1)*(m+1)<=n) ++m;
	initp(m);
	S[0][0]=1;
	for (int i=1;i<=maxk;++i)
		for (int j=1;j<=i;++j)
			S[i][j]=(S[i-1][j-1]+j*S[i-1][j])%mo;
	inv[1]=1;
	for (int i=2;i<=maxk+1;++i)
		inv[i]=(mo-mo/i)*inv[mo%i]%mo;
//	doit(maxk);
	for (int i=1;i<=maxk;++i){
		doit(i);
		printf("%d
",i);
	}
	ans=(ans+mo)%mo;
	printf("%lld
",ans);
	//ans=797866893
	return 0;
}