Solution
这个的话。。看到跟(gcd)相关这种东西。。那肯定是要先反演一波啦
(sgcd(i,j))这个东西看起来很麻烦不过其实就是(gcd(i,j))的次大约数
为了方便表示我们用(f(x))表示(x)的次大约数，那么原来的式子可以写成：
然后倒数第二部到最后一步是直接由(phi)的定义得到的嗯（然后我就特别弱智地用(mu)推了好久还要推错了。。）
然后现在我们就是要求出：
这个东西
后半部分的话是经典的杜教筛求欧拉函数前缀和，然后因为(i)的上界是(lfloor frac{n}{d} floor)所以求出前缀和以后分块搞一波就好了
关于杜教筛求欧拉函数前缀和。。具体一点的话。。。戳这里好了
然后前半部分的求法十分神仙（疯狂%%%），为了防止变量名称弄混，统一一下接下来的部分题目中给出的指数为(K)，而用来枚举的变量为(k)
首先定义
然后我们令min_25中的那个(g)数组按照如下的定义计算：
然后我们思考一下在计算过程中，(g_{i,j}=g_{i,j-1}-s(P_j)cdot(g_{lfloorfrac{i}{P_j}
floor,j-1}-g_{P_j-1,j-1}))这一步中，括号中的那一部分的含义
其实就是(sumlimits_{x}s(x)[minp(x)=P_j且minp(frac{x}{P_j})>=P_j])，然后这个其实就是这堆满足条件的(x)的(f(x)^k)之和
然后又因为min_25筛中，每个数只会被考虑一次，所以我们直接多开一个数组把这堆东西加起来就好了
但是这并没有算上(kin Prime)的情况，所以此外我们还要单独算一下素数情况下的取值，其实也就是素数的个数（因为(f(prime)=1)），这个也是直接用min_25筛一下就好了
最后还剩下一个小问题，(g_{i,0})在初始化的时候要快速求出(sumlimits_{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,&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
",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)];
}
}
}