题目链接
sol
[ans=sum_{i=1}^nsum_{x_1,x_2,ldots,x_m=1}^ioperatorname{lcm}(gcd(i,x_1),gcd(i,x_2),ldots,gcd(i,x_m))\
=sum_{i=1}^nsum_{x_1,x_2,ldots,x_m=1}^ifrac{i}{gcd(frac{i}{gcd(i,x_1)},frac{i}{gcd(i,x_2)},ldots,frac{i}{gcd(i,x_m)})}\]
然后考虑(y_j=frac i{gcd(i,x_j)})，容易发现它出现的次数等于(varphi(y_j))。（左式简单变换即可
即求
[ans=sum_{i=1}^nsum_{y_1,y_2,ldots,y_mmid i}varphi(y_1)varphi(y_2)cdotsvarphi(y_m)frac{i}{gcd(y_1,y_2,ldots,y_m)}
]
令(f(n)=sum_{y_1,y_2,ldots,y_mmid n}varphi(y_1)varphi(y_2)cdotsvarphi(y_m)frac{n}{gcd(y_1,y_2,ldots,y_m)})
然后这好像就是个积性函数？讲真我看不出来有谁会证评论教教我好不好呀(QwQ)。
下面我们假装它是积性函数。
实际上我们只要会算(f(p^k))就行了，和他们一样，我们暂时用(d)来代替(k)。
[f(p^k)=sum p^{max(x1,x2,...,xd)}prodvarphi(p^{k-x_i})
]
设(mx)为最大指数
[g(mx)=sum_{i=0}^{mx}(varphi(x^{k-i}))^d\
=egin{cases}
p^{kd}&,mx=1\
(p^k-p^{mx-1})^d&,2leq mxleq k
end{cases}]
则
[f(p^k)=sum_{mx=0}^kp^{mx}(g(mx)-g(mx-1))\
f(p^{k+1})=p^df(p^k)-p^k(p^{(k+1)d}-(p^{k+1}-1)^d)+p^{k+1}(p^{(k+1)d}-(p^{k+1}-1)^d)\
=p^df(p^k)+p^k(p-1)(p^{(k+1)d}-(p^{k+1}-1)^d)]
就是考虑(g(mx))中的(k)的变化，暴力推。
现在我们可以做(nle10^7)了
关于(n>10^7)，有(k<100)于是我们想到(Min\_25)筛。
[f(p)=p^{d+1}-(p-1)^{d+1}\
=sum_{i=0}^{d}(-1)^{d-i}inom{d+1}{i}p^i]
(f(p^k))暴力算。
然后就可以(Min\_25)筛了。
关于这份代码，我一开始写斯特林数处理自然数幂和，莫名(WA)了，于是蒯了yww的插值做法，恩。
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define gt getchar()
#define ll long long
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
inline int in()
{
int k=0;char ch=gt;
while(ch<'-')ch=gt;
while(ch>'-')k=k*10+ch-'0',ch=gt;
return k;
}
const int YL=1e9+7;
inline int MO(const int &x){return x>=YL?x-YL:x;}
inline int OD(const int &x){return x<0?x+YL:x;}
inline int ksm(int a,ll k){int r=1;while(k){if(k&1)r=1ll*r*a%YL;a=1ll*a*a%YL,k>>=1;}return r;}
int n,k,d;
namespace w1
{
const int N=1e7+5;bool np[N];
int pr[N],tot,pw[N],f[N],c[N];
void work()
{
f[1]=pw[1]=1;
for(int i=2;i<=n;++i)
{
if(!np[i])pr[++tot]=i,pw[i]=ksm(i,k);
for(int j=1;i*pr[j]<=n;++j)
{
np[i*pr[j]]=1;pw[i*pr[j]]=1ll*pw[i]*pw[pr[j]]%YL;
if(i%pr[j]==0)break;
}
}
for(int i=2;i<=n;++i)
{
if(!np[i])c[i]=i,f[i]=OD(ksm(i,k+1)-ksm(i-1,k+1));
for(int j=1;i*pr[j]<=n;++j)
{
int v=i*pr[j];
if(i%pr[j]==0)
{
c[v]=c[i]*pr[j];
if(c[v]==v)
f[v]=(1ll*pw[pr[j]]*f[i]%YL+1ll*i*(pr[j]-1)%YL*OD(pw[v]-pw[v-1])%YL)%YL;
else f[v]=1ll*f[c[v]]*f[v/c[v]]%YL;break;
}
c[v]=pr[j],f[v]=1ll*f[i]*f[pr[j]]%YL;
}
}
ll ans=0;for(int i=1;i<=n;++i)ans+=f[i];
printf("%lld
",ans%YL);
}
}
namespace w2
{
const int N=2e5+5,M=N-5;bool np[N];
int pr[N],tot,pw[N],ans[N],g[N],id1[N],id2[N],sp[N];
int fnv[N],fac[N],s[N],w[N],m,S[105][105],sum[N];
inline int C(int x,int y){return y>x?0:1ll*fac[x]*fnv[y]%YL*fnv[x-y]%YL;}
void init()
{
fac[0]=1;
for(int i=1;i<=M;++i)fac[i]=1ll*fac[i-1]*i%YL;fnv[M]=ksm(fac[M],YL-2);
for(int i=M;i>=1;--i)fnv[i-1]=1ll*fnv[i]*i%YL;
for(int i=2;i<=M;++i)
{
if(!np[i])pr[++tot]=i,sum[tot]=MO(sum[tot-1]+OD(ksm(i,k+1)-ksm(i-1,k+1)));
for(int j=1;i*pr[j]<=M;++j)
{np[i*pr[j]]=1;if(i%pr[j]==0)break;}
}
}
void pre_(int x)
{
int cnt=0;s[1]=1;
for(int i=2;i<=M;++i)
{
if(!np[i])pw[i]=ksm(i,x),++cnt,sp[cnt]=MO(sp[cnt-1]+pw[i]);
for(int j=1;i*pr[j]<=M;++j)
{pw[i*pr[j]]=1ll*pw[i]*pw[pr[j]]%YL;if(i%pr[j]==0)break;}
s[i]=MO(s[i-1]+pw[i]);
}
}
int pre[N],suf[N];
int calc(int x,int y)
{
if(x<=y+2)return s[x];
pre[0]=1;
for(int i=1;i<=y+2;i++)pre[i]=1ll*pre[i-1]*(x-i)%YL;
suf[y+3]=1;
for(int i=y+2;i>=1;i--)suf[i]=1ll*suf[i+1]*(x-i)%YL;
ll res=0;
for(int i=1;i<=y+2;i++)res=(res+1ll*s[i]*pre[i-1]%YL*suf[i+1]%YL*fnv[i-1]%YL*fnv[y+2-i]%YL*((y+2-i)&1?YL-1:1))%YL;
return (res%YL+YL)%YL;
}
#define ID(x) ((x)<=d?id1[x]:id2[n/(x)])
void deal(int x)
{
pre_(x);
for(int i=1;i<=m;++i)g[i]=((w[i]<=d)?s[w[i]]:calc(w[i],x))-1;
for(int i=1;pr[i]<=n/pr[i];++i)
for(int j=1;j<=m&&pr[i]<=w[j]/pr[i];++j)
g[j]=OD(g[j]-1ll*pw[pr[i]]*OD(g[ID(w[j]/pr[i])]-sp[i-1])%YL);
}
inline int G(int z,int x,int y)
{
if(z<0)return 0;if(z==y)return ksm(x,1ll*y*k);
return ksm(OD(ksm(x,y)-ksm(x,y-z-1)),k);
}
inline int Get(int x,int y)
{
int res=0;
for(int i=0;i<=y;++i)
res=(1ll*ksm(x,i)*OD(G(i,x,y)-G(i-1,x,y))+res)%YL;
return res;
}
int solve(int x,int y)
{
if(x<=1||pr[y]>x)return 0;
int res=OD(ans[ID(x)]-sum[y-1]);
for(int i=y;i<=tot&&1ll*pr[i]*pr[i]<=x;++i)
for(int e=1,p=pr[i];1ll*p*pr[i]<=x;++e,p*=pr[i])
res=MO(res+(1ll*solve(x/p,i+1)*Get(pr[i],e)+Get(pr[i],e+1))%YL);
return res;
}
void work()
{
d=sqrt(n);init();pw[1]=1;
for(int i=1;i<=n;i=n/(n/i)+1)w[++m]=n/i,ID(w[m])=m;
for(int i=0;i<=k;++i)
{
deal(i);int v=1ll*C(k+1,i)*((k-i)&1?YL-1:1)%YL;
for(int j=1;j<=m;++j)ans[j]=(1ll*v*g[j]+ans[j])%YL;
}
int ans=solve(n,1)+1;printf("%d
",ans);
}
}
int main()
{
n=in(),k=in();
if(n<=1e7)w1::work();
else w2::work();
return 0;
}