Solution
這個的話。。看到跟\(gcd\)相關這種東西。。那肯定是要先反演一波啦
\(sgcd(i,j)\)這個東西看起來很麻煩不過其實就是\(gcd(i,j)\)的次大約數
為了方便表示我們用\(f(x)\)表示\(x\)的次大約數，那么原來的式子可以寫成：
然后倒數第二部到最后一步是直接由\(\phi\)的定義得到的嗯（然后我就特別弱智地用\(\mu\)推了好久還要推錯了。。）
然后現在我們就是要求出：
這個東西
后半部分的話是經典的杜教篩求歐拉函數前綴和，然后因為\(i\)的上界是\(\lfloor \frac{n}{d}\rfloor\)所以求出前綴和以后分塊搞一波就好了
關於杜教篩求歐拉函數前綴和。。具體一點的話。。。戳這里好了
然后前半部分的求法十分神仙（瘋狂%%%），為了防止變量名稱弄混，統一一下接下來的部分題目中給出的指數為\(K\)，而用來枚舉的變量為\(k\)
首先定義
然后我們令min_25中的那個\(g\)數組按照如下的定義計算：
然后我們思考一下在計算過程中，\(g_{i,j}=g_{i,j-1}-s(P_j)\cdot(g_{\lfloor\frac{i}{P_j}\rfloor,j-1}-g_{P_j-1,j-1})\)這一步中，括號中的那一部分的含義
其實就是\(\sum\limits_{x}s(x)[minp(x)=P_j且minp(\frac{x}{P_j})>=P_j]\)，然后這個其實就是這堆滿足條件的\(x\)的\(f(x)^k\)之和
然后又因為min_25篩中，每個數只會被考慮一次，所以我們直接多開一個數組把這堆東西加起來就好了
但是這並沒有算上\(k\in Prime\)的情況，所以此外我們還要單獨算一下素數情況下的取值，其實也就是素數的個數（因為\(f(prime)=1\)），這個也是直接用min_25篩一下就好了
最后還剩下一個小問題，\(g_{i,0}\)在初始化的時候要快速求出\(\sum\limits_{j=1}^{i}j^k\)
注意到\(k\)很小，那么我們用第二類斯特林數處理一下：
然后就很愉快滴解決啦
一些實現的細節（你別說個人感覺因為那個模數不是很友善所以寫起來有點麻煩。。）
1、取模的話可以自然溢出，\(unsigned\ int\)了解一下ovo
2、求下降階乘冪的時候，因為這題的模數不是質數，所以不能快樂求逆元，但是因為\(k\)比較小所以直接暴力一波，一個一個數乘，這堆數中必定有一個是\(j+1\)的倍數，所以可以直接先除再乘就好了，注意判斷的時候每次都取模判斷的話會愉快T掉，所以應該在一開始的時候先把余數算出來，然后根據余數的加法定理直接判就好了
代碼的話大概長這個樣子
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<cmath>
#define ll long long
#define ui unsigned int
using namespace std;
const int MAXN=1e5+10;
map<ll,ui>Rec;
ui P[MAXN],loc1[MAXN],loc2[MAXN];
ui g[MAXN*2],rec[MAXN*2],sphi[MAXN],pcnt[MAXN*2],pw[MAXN*2],G[MAXN*2];
ui S[60][60];
bool vis[MAXN];
ll n,m,K,cnt,cntval,sq;
ui ans;
void prework(int n);
void init_loc(ll n);
void get_g();
int Pos(ll x){return x<=sq?loc1[x]:loc2[n/x];}
ui Sum(ui l,ui r);
ui SumPhi(ll n);
ui calc_sumk(ll n);
ui ksm(ui x,ll y);
ui f(ll x){return G[Pos(x)]+pcnt[Pos(x)];}
int main(){
#ifndef ONLINE_JUDGE
freopen("a.in","r",stdin);
#endif
scanf("%lld%lld\n",&n,&K);
prework(MAXN);
sq=sqrt(n);
m=upper_bound(P+1,P+1+cnt,sq)-P-1;
init_loc(n);
get_g();
ans=0;
ui pre=0,now;
for (ll i=2,pos;i<=n;i=pos+1){
pos=n/(n/i);
now=f(i);
ans+=(SumPhi(n/i)*2-1)*(now-pre);
pre=now;
}
printf("%lu\n",ans);
}
ui ksm(ui x,ll y){
ui ret=1,base=x;
for (;y;y>>=1,base=base*base)
if (y&1) ret=ret*base;
return ret;
}
void prework(int n){
cnt=0;
sphi[1]=1;
for (int i=2;i<=n;++i){
if (!vis[i])
P[++cnt]=i,sphi[i]=i-1;
for (int j=1;j<=cnt&&P[j]*i<=n;++j){
vis[i*P[j]]=true;
if (i%P[j]==0){
sphi[i*P[j]]=sphi[i]*P[j];
break;
}
else
sphi[i*P[j]]=sphi[i]*sphi[P[j]];
}
}
for (int i=1;i<=n;++i) sphi[i]=sphi[i-1]+sphi[i];
S[0][0]=1;
for (int i=1;i<=K;++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];
}
for (int i=1;i<=cnt;++i) pw[i]=ksm(P[i],K);
}
void init_loc(ll n){
cntval=0;
for (ll i=1,pos;i<=n;i=pos+1){
rec[++cntval]=n/i;
pos=n/(n/i);
}
reverse(rec+1,rec+1+cntval);
for (int i=1;i<=cntval;++i)
if (rec[i]<=sq) loc1[rec[i]]=i;
else loc2[n/rec[i]]=i;
}
ui Sum(ui l,ui r){
ui tmp1=l+r,tmp2=r-l+1;
if (tmp1%2==0) tmp1/=2;
else tmp2/=2;
return tmp1*tmp2;
}
ui SumPhi(ll n){
if (n<MAXN) return sphi[n];
map<ll,ui>::iterator tmp=Rec.find(n);
if (tmp!=Rec.end()) return tmp->second;
ui ret=Sum(1,n);
for (ll i=2,pos;i<=n;i=pos+1){
pos=n/(n/i);
ret-=SumPhi(n/i)*(pos-i+1);
}
Rec[n]=ret;
return ret;
}
ui calc_sumk(ll n){
ui ret=0,tmpval;
int tmp;
for (int i=1;i<=K;++i){
tmpval=1;
tmp=(n-i+1)%(i+1);
for (ll j=n-i+1;j<=n+1;++j,++tmp){
if (tmp>=i+1) tmp-=i+1;
if (!tmp) tmpval*=(ui)j/(i+1);
else tmpval*=(ui)j;
}
ret+=S[K][i]*tmpval;
}
return ret;
}
void get_g(){
for (int i=2;i<=cntval;++i) g[i]=calc_sumk(rec[i])-1,pcnt[i]=rec[i]-1,G[i]=0;
for (int j=1;j<=m;++j){
for (int i=cntval;i>=1&&rec[i]>=1LL*P[j]*P[j];--i){
g[i]-=pw[j]*(g[Pos(rec[i]/P[j])]-g[Pos(P[j]-1)]);
pcnt[i]-=pcnt[Pos(rec[i]/P[j])]-pcnt[Pos(P[j]-1)];
G[i]+=g[Pos(rec[i]/P[j])]-g[Pos(P[j]-1)];
}
}
}