Link
设(f(n))表示(n)的次大因子，那么(sgcd(i,j)=f(gcd(i,j)))。
[egin{aligned}
ans&=sumlimits_{i=1}^nsumlimits_{j=1}^nf(gcd(i,j))^k\
&=sumlimits_{i=1}^nf(i)^ksumlimits_{u=1}^nsumlimits_{v=1}^n[gcd(i,j)=i]\
&=sumlimits_{i=1}^nf(i)^ksumlimits_{u=1}^{lfloorfrac ni
floor}sumlimits_{v=1}^{lfloorfrac ni
floor}sumlimits_{d|gcd(u,v)}mu(d)\
&=sumlimits_{i=1}^nf(i)^ksumlimits_{d=1}^{lfloorfrac ni
floor}mu(d)lfloorfrac n{id}
floor^2
end{aligned}
]
设(g=f^k imesmu)（幂在点乘意义下），那么
[egin{aligned}
ans&=sumlimits_{i=1}^nf(i)^ksumlimits_{d=1}^{lfloorfrac nk
floor}mu(d)lfloorfrac n{id}
floor^2\
&=sumlimits_{i=1}^ng(i)lfloorfrac ni
floor^2
end{aligned}
]
对(lfloorfrac ni
floor^2)整除分块，那么我们要计算的就是(g)在所有(lfloorfrac ni
floor)的前缀和。
注意到(g imes I=f)，因此如果我们能够计算(f^k)在所有(lfloorfrac ni
floor)的前缀和的前缀和，就能够用杜教筛计算(sum g)了。
[sumlimits_{i=1}^nf(i)^k=pi(n)+sumlimits_{plesqrt n}sumlimits_{i=1}^{lfloorfrac np
floor}[min(i)ge p]i^k
]
其中(min(i))表示(i)的最小非平凡质因子，若(i=1vee iinmathbb P)则(min(i)=-infty)。
后面和Min_25筛计算(sum p^k)时求(h)的过程中被筛掉的部分比较类似，模拟即可。
同时可以顺便求出需要用到的(pi(n))。
求(h)的过程中需要求自然数幂和，考虑利用Stirling数处理。
[egin{aligned}
&sumlimits_{i=1}^ni^k\
=&sumlimits_{i=1}^nsumlimits_{j=0}^kleft{katop j
ight}i^{underline j}\
=&sumlimits_{j=0}^kleft{katop j
ight}frac{(n+1)^{underline{j+1}}}{j+1}
end{aligned}
]
(j+1)不一定有逆元，注意到显然有(j+1|(n+1)^{underline{j+1}})，因此暴力计算下降幂时除掉即可。
#include<cstdio>
using u32=unsigned;
const int N=100007,K=57;
int n,k,tot,lim,cnt,pr[N],pos[N];u32 S[K][K],pw[N],sum_pw[N],sum[N];
int read(){int x;scanf("%d",&x);return x;}
int getid(int x){return x<=lim? x:cnt+1-n/x;}
u32 pow(u32 p,int e){u32 q=1;for(;e;e>>=1,p*=p)if(e&1)q*=p;return q;}
u32 power_sum(int n)
{
u32 ans=0;
for(int i=1;i<=k&&i<=n;++i)
{
u32 now=1;int f=0;
for(int j=n+1;j>n-i;--j) if(!(j%(i+1))&&!f) now*=j/(i+1),f=1; else now*=j;
ans+=now*S[k][i];
}
return ans;
}
void init()
{
static int low[N];S[0][0]=1;
for(lim=1;(lim+1)*(lim+1)<=n;++lim);
for(int i=1;i<=k;++i) for(int j=1;j<=i;++j) S[i][j]=S[i-1][j-1]+S[i-1][j]*j;
for(int i=2;i<=lim;++i)
{
if(!low[i]) pr[++tot]=low[i]=i,sum_pw[tot]=sum_pw[tot-1]+(pw[tot]=pow(i,k));
for(int j=1;j<=tot&&i*pr[j]<=lim&&pr[j]<=low[i];++j) low[i*pr[j]]=pr[j];
}
}
void Min_25_Sieve()
{
static u32 pi[N],t[N];
for(int l=1,r;l<=n;l=r+1) r=n/(n/l),pos[++cnt]=r,pi[cnt]=r-1,t[cnt]=power_sum(r)-1;
for(int i=1;i<=tot&&pr[i]<=lim;++i) for(int j=cnt,p;j&&pr[i]*pr[i]<=pos[j];--j) p=getid(pos[j]/pr[i]),pi[j]-=pi[p]-(i-1),sum[j]+=t[p]-sum_pw[i-1],t[j]-=pw[i]*(t[p]-sum_pw[i-1]);
for(int i=1;i<=cnt;++i) sum[i]+=pi[i];
}
void Xudyh_Sieve()
{
for(int i=1;i<=cnt;++i) for(int l=2,r;l<=pos[i];l=r+1) r=pos[i]/(pos[i]/l),sum[i]-=(r-l+1)*sum[getid(pos[i]/l)];
}
void calc()
{
u32 ans=0;
for(int l=1,r,p;l<=n;l=r+1) p=getid(r=n/(n/l)),ans+=(sum[p]-sum[p-1])*(n/l)*(n/l);
printf("%u",ans);
}
int main()
{
n=read(),k=read();
init(),Min_25_Sieve(),Xudyh_Sieve(),calc();
}