分析:
考场上写min25筛T到飞起,插值的时候脑抽没优化用的\(O(K^2)\),结果跟暴力差不多速度
由于\(f(p^k)\)是\(O(1)\)可求的,答案:
\(\sum_{i=1}^{n}f(i)=G_k(n,|P|)\)
所以无脑筛就好了。。。
答案进行了优化,引入了一个叫powerful number的奥妙重重的东西
这里可以去拜读zzq大神的博客(在线膜拜)
我们要求的是一个积性函数\(F(p)=p^k\)
我们先考虑另外一个积性函数\(G(x)=x^k\)
注意\(p\)指的是指数,而\(x\)指正整数
显然当\(k\)很小时,\(G(x)\)可以使用拉格朗日插值法\(O(k)\)求出
接下来先介绍powerful number(多次方很多的数(大雾)
powerful number的定义是每个质因子次数都大于1的数。首先,小于\(n\)的powerful number个数是\(O(\sqrt{n})\)的,这是因为所有powerful number的质因子都含次方,开根后对应的数是\(O(\sqrt{n})\)的
于是可以通过爆搜求出全部powerful number
我们现在想要找到一个函数\(H\),使得\(F=G*H\)(这里指狄利克雷卷积),所以:
\(F(i)=\sum_{d|i}G(d)H(\frac{i}{d})\)
\(\sum_{i=1}^{n}F(i)=\sum_{ij\leq n}G(i)H(j)=\sum_{i=1}^{n}H(i)\sum_{j=1}^{\frac{n}{i}}G(j)\)
因为\(F(p)=G(p)=p^k\)
所以\(H(1)=1\)且\(H(p)=0\)
由于\(H\)为积性函数,所以\(H\)除1以外不为0的位置只可能是powerful number
于是枚举powerful number,算出对应位置的\(H\)就好了
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<iostream>
#include<map>
#include<string>
#define maxn 10000005
#define INF 0x3f3f3f3f
#define MOD 1000000007
using namespace std;
inline long long getint()
{
long long num=0,flag=1;char c;
while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;
while(c>='0'&&c<='9')num=num*10+c-48,c=getchar();
return num*flag;
}
long long n,S;
int K;
int pri[maxn],np[maxn],cnt;
int prik[maxn],sumprik[maxn];
int pre[maxn],sub[maxn];
int fac[maxn],inv[maxn],y[maxn];
int pos1[maxn],pos2[maxn],tot,sum[maxn];
long long num[maxn];
inline int add(int x,int y){return x+y<MOD?x+y:x+y-MOD;}
inline int ksm(int num,int k)
{
int ret=1;
for(;k;k>>=1,num=1ll*num*num%MOD)if(k&1)ret=1ll*ret*num%MOD;
return ret;
}
inline void init(int N)
{
for(int i=2;i<=N;i++)
{
if(!np[i])pri[++cnt]=i,prik[cnt]=ksm(pri[cnt],K),sumprik[cnt]=add(sumprik[cnt-1],prik[cnt]);
for(int j=1;j<=cnt&&i*pri[j]<=N;j++){np[i*pri[j]]=1;if(i%pri[j]==0)break;}
}
}
inline int &id(long long x){return x<=S?pos1[x]:pos2[n/x];}
inline void dfs(long long x,int cur,int p)
{
sum[id(n/x)]=add(sum[id(n/x)],p);
for(int i=cur+1;i<=cnt&&x<=n/pri[i]/pri[i];i++)
{
long long tmp=pri[i];
do{
tmp*=pri[i];
dfs(x*tmp,i,1ll*add(prik[i],MOD-1ll*prik[i]*prik[i]%MOD)*p%MOD);
}while(tmp<=n/x/pri[i]);
}
}
inline int Sumk(int m)
{
int N=K+2,ans=0;
pre[0]=sub[N+1]=1;
for(int i=1;i<=N;i++)pre[i]=1ll*pre[i-1]*(m-i)%MOD;
for(int i=N;i>=1;i--)sub[i]=1ll*sub[i+1]*(m-i)%MOD;
for(int i=1;i<=N;i++)
{
int tmp=1ll*pre[i-1]*sub[i+1]%MOD,p=1ll*((N-i)&1?MOD-1:1)*inv[i-1]%MOD*inv[N-i]%MOD;
ans=add(ans,1ll*tmp*y[i]%MOD*p%MOD);
}
return ans;
}
int main()
{
n=getint(),K=getint();
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=1;i<=K+2;i++)y[i]=add(y[i-1],ksm(i,K));
for(int i=2;i<=K+2;i++)fac[i]=1ll*fac[i-1]*i%MOD;
for(int i=2;i<=K+2;i++)inv[i]=1ll*(MOD-MOD/i)*inv[MOD%i]%MOD;
for(int i=2;i<=K+2;i++)inv[i]=1ll*inv[i]*inv[i-1]%MOD;
init(S=sqrt(n+1));
for(long long i=1;i<=n;i=n/(n/i)+1)num[id(n/i)=++tot]=n/i;
dfs(1,0,1);
int ans=0;
for(int i=1;i<=tot;i++)if(sum[i])ans=add(ans,1ll*sum[i]*Sumk(num[i]%MOD)%MOD);
printf("%d\n",(ans%MOD+MOD)%MOD);
}