有积性函数(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;
}