记\(f(x)=\)\(x\)的次大因数，那么\(sgcd(i,j)=f(gcd(i,j))\)。
下面来推式子：
\[ \begin{aligned} \sum_{i=1}^n\sum_{j=1}^nsgcd(i,j)^k&=\sum_{i=1}^n\sum_{j=1}^nf(gcd(i,j))^k&记d=gcd(i,j)\\ &=\sum_{d=1}^nf(d)^k\sum_{\frac i d=1}^{\lfloor \frac n d\rfloor}\sum_{\frac j d=1}^{\lfloor \frac n d\rfloor}[gcd(\frac i d,\frac j d)=1]&f(1)=0\\ &=\sum_{d=2}^nf(d)^k\sum_{i=1}^{\lfloor \frac n d\rfloor}\sum_{j=1}^{\lfloor \frac n d\rfloor}[gcd(i,j)=1]\\ &=\sum_{d=2}^nf(d)^k(2\sum_{i=1}^{\lfloor \frac n d\rfloor}\varphi(i)-1) \end{aligned} \]
最后一步的括号是用欧拉函数的定义直接替换出来的。
我们发现，可以按\(\lfloor \frac n d \rfloor\)的取值数论分段，因为括号显然只受\(\lfloor \frac n d\rfloor\)的取值影响，关键是如何求\(f(x)^k\)的前缀和？
记\(S_x=\sum_{i=2}^xf(x)^k\)，
令min_25中的\(g\)数组以\(s(x)=x^k\)的定义计算\(g\)，考虑由\(g_{i,j-1}\)推到\(g_{i,j}\)的时候，减去的是什么？是\(s(p_j)*(g_{\lfloor\frac n {p_j}\rfloor,j-1}-g_{p_j-1,j-1})\)，后面括号的部分是什么呢？恰好是那些最小质因子为\(p_j\)且除去\(p_j\)后剩余部分最小质因数不小于\(p_j\)的数的\(k\)次方之和，我们发现这些数的\(f\)值之和就是后面的括号。又因为每一个数至多被如此枚举到1次，所以对于每一个\(i\)，我们把括号的值累加起来，这就是那些\([2,i]\)中非质数的\(f\)值之和；再加上\([2,i]\)中的质数的\(f\)值之和，也就是质数个数，我们就求得了\(S_i\)。
min_25真牛逼。
所以就做完了，\(g\)的初始化需要用到自然数幂求和，考虑到\(k\)比较小，可以用第二类斯特林数求解。
\[ \begin{aligned} Sum_k(n)&=\sum_{i=1}^ni^k\\ &=\sum_{j=1}^k\begin{Bmatrix}k\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1} \end{aligned} \]
Code
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <map>
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef map<ll,ui> mlu;
const int SQRTN=32000,LEN=1000001;
bool vis[LEN];
ui p[LEN],pcnt,pphi[LEN],s[51][51],pk[LEN];
ll n,k,sqrtn,m;
ll a[SQRTN*2],cnt,pos1[SQRTN],pos2[SQRTN];
ui g[SQRTN*2],sum[SQRTN*2],g1[SQRTN*2];
mlu record;
ui ksm(ui x,int y){
ui res=1;
for(;y;x=x*x,y>>=1)
if(y&1) res=res*x;
return res;
}
void prework(){
pphi[1]=1;
for(int i=2;i<LEN;i++){
if(!vis[i]){
p[++pcnt]=i;
pphi[i]=i-1;
}
for(int j=1;j<=pcnt&&i*p[j]<LEN;j++){
int x=i*p[j];
vis[x]=true;
if(i%p[j]==0){
pphi[x]=pphi[i]*p[j];
break;
}
pphi[x]=pphi[i]*(p[j]-1);
}
}
for(int i=2;i<LEN;i++) pphi[i]+=pphi[i-1];
s[0][0]=s[0][1]=1;
for(int i=1;i<=50;i++){
s[i][1]=1;
for(int j=2;j<=i;j++)
s[i][j]=s[i-1][j]*(ui)j+s[i-1][j-1];
}
}
inline int gp(ll x){
return x<=sqrtn?pos1[x]:pos2[n/x];
}
void Diz(){
for(ll i=1,j;i<=n;i=j+1){
j=n/(n/i);
a[++cnt]=n/i;
}
reverse(a+1,a+1+cnt);
for(int i=1;i<=cnt;i++)
if(a[i]<=sqrtn) pos1[a[i]]=i;
else pos2[n/a[i]]=i;
}
ui sumup(ll l,ll r){
ll a=l+r,b=r-l+1;
if(a&1) b/=2; else a/=2;
return (ui)a*(ui)b;
}
ui getPhi(ll n){
if(n<LEN) return pphi[n];
mlu::iterator pos=record.find(n);
if(pos!=record.end()) return pos->second;
ui res=sumup(1,n);
for(ll i=2,j;i<=n;i=j+1){
j=n/(n/i);
res-=getPhi(n/i)*(j-i+1);
}
record[n]=res;
return res;
}
ui getSumk(ll n){
ui res=0,t;
for(ll j=1;j<=k;j++){
t=1;
int md=(n-j+1)%(j+1);
for(ll i=n-j+1;i<=n+1;i++,md++)
if(!md||md==(j+1)) t*=(ui)i/(j+1);
else t*=(ui)i;
res+=t*s[k][j];
}
return res;
}
void calc_g(){
for(int i=2;i<=cnt;i++) g[i]=getSumk(a[i])-1,g1[i]=a[i]-1,sum[i]=0;
for(int j=1;j<=m;j++)
for(int i=cnt;i>=2&&a[i]>=p[j]*p[j];i--){
g1[i]-=g1[gp(a[i]/p[j])]-g1[gp(p[j]-1)];
ui delta=(g[gp(a[i]/p[j])]-g[gp(p[j]-1)]);
g[i]-=pk[j]*delta;
sum[i]+=delta;
}
}
ui calc(ll n){
return sum[gp(n)]+g1[gp(n)];
}
int main
prework();
scanf("%lld%lld",&n,&k);
sqrtn=(ll)sqrt(n);
m=upper_bound(p+1,p+1+pcnt,sqrtn)-p-1;
for(int i=1;i<=m;i++) pk[i]=ksm(p[i],k);
Diz();
calc_g();
ui ans=0,last=0,tmp;
for(ll i=2,j;i<=n;i=j+1){
j=n/(n/i);
tmp=calc(j);
ans+=(getPhi(n/i)*2-1)*(tmp-last);
last=tmp;
}
printf("%lu\n",ans);
return 0;
}