设 \(f(n)\) 表示 \(n\) 的次大因数。
\[\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^nf(\gcd(i,j))^k\\ =&\sum_{d=1}^nf(d)^k\sum_{i=1}^{(n/d)}\sum_{j=1}^{(n/d)}[\gcd(i,j)=1]\\ =&\sum_{d=1}^nf(d)^k\left(2\sum_{i=1}^{(n/d)}\varphi(i)-1\right) \end{aligned} \]（其中第二步推到第三步是根据 \(\varphi\) 的定义，括号内 \(-1\) 是因为 \(i=j=1\) 的情况我们算了两次）
我们数论分块，然后后面括号里的 \(\varphi\) 的前缀和可以用杜教筛解决。
对于前面的 \(f(d)^k\) 的前缀和，我们分为质数和合数两部分。
对于质数，\(f(p)=1\)，于是询问的就是区间质数个数，这个可以用 min_25 筛解决。
对于合数，观察 min_25 筛中推导 \(g_k(n,i)=\sum\limits_{x=1}^n[x\in \mathbb{P}\lor \operatorname{lpf}(x)>p_i]x^k\) 的过程：\(g_k(n,i)\) 由 \(g_k(n,i-1)\) 减去最小质因数为 \(p_i\) 的合数的贡献得到。所以我们只需要加上求 \(g_k(n,i)\) 过程中用 \(g_k(n,i-1)\) 减去的东西即可。
但是注意初始化 \(g_k(n,0)=\sum\limits_{x=2}^nx^k\) 时，不能直接用拉格朗日插值来求自然数幂求和，因为模数 \(2^{32}\) 不支持取逆。
那么就需要用到第二类斯特林数了：
\[\sum_{i=0}^ni^k=\sum_{j=0}^k\begin{Bmatrix}k\\j\end{Bmatrix}\dfrac{(n+1)^{\underline{j+1}}}{j+1} \]我们 \(O(k^2)\) 预处理第二类斯特林数，然后每次 \(O(k^2)\) 查询。
#include<bits/stdc++.h>
#define K 55
#define ll unsigned int
using namespace std;
inline ll poww(ll a,int b)
{
int ans=1;
while(b)
{
if(b&1) ans=ans*a;
a=a*a;
b>>=1;
}
return ans;
}
inline int read()
{
int x=0,f=1;
char ch=getchar();
while(ch<'0'||ch>'9')
{
if(ch=='-') f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9')
{
x=(x<<1)+(x<<3)+(ch^'0');
ch=getchar();
}
return x*f;
}
int n,k;
namespace s1
{
const int maxn=100000,N=maxn+10;
const int SN=10010;
int cnt,p[N];
bool notprime[N];
ll phi[N],sphi[SN];
void init()
{
phi[1]=1;//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
for(int i=2;i<=maxn;i++)
{
if(!notprime[i])
{
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1,x;j<=cnt&&(x=i*p[j])<=maxn;j++)
{
notprime[x]=1;
if(!(i%p[j]))
{
phi[x]=phi[i]*p[j];
break;
}
phi[x]=phi[i]*phi[p[j]];
}
}
for(int i=2;i<=maxn;i++)
phi[i]=phi[i-1]+phi[i];
}
int get(int x)
{
return n/x;
}
ll S(int n)
{
if(n<=maxn) return phi[n];
if(sphi[get(n)]) return sphi[get(n)];
ll ans=1ll*n*(n+1)/2;
for(int l=2,r;l<=n;l=r+1)
{
r=n/(n/l);
ans-=S(n/l)*(r-l+1);
}
return sphi[get(n)]=ans;
}
ll query(int n)
{
return S(n);
}
}
namespace s2
{
const int N=32000;
int sn,tot,b[N<<1],id1[N],id2[N];
int cnt,p[N];
bool notprime[N];
ll s[K][K];
ll powk[N],gpk[N];
ll g0[N<<1],gk[N<<1];
ll ans1[N<<1];
void init()
{
for(int i=2;i<=sn;i++)
{
if(!notprime[i])
{
p[++cnt]=i,powk[cnt]=poww(i,k);
gpk[cnt]=gpk[cnt-1]+powk[cnt];
}
for(int j=1,x;j<=cnt&&(x=i*p[j])<=sn;j++)
{
notprime[x]=1;
if(!(i%p[j])) break;
}
}
}
ll getsum(int n)
{
ll ans=0;
for(int j=0;j<=k;j++)
{
ll tmp=1;
for(int l=0;l<=j;l++)
{
if((n+1-l)%(j+1)) tmp*=(n+1-l);
else tmp*=(n+1-l)/(j+1);
}
ans+=s[k][j]*tmp;
}
return ans;
}
int get(int x)
{
return x<=sn?id1[x]:id2[n/x];
}
void work()
{
sn=sqrt(n);
init();
s[0][0]=1;
for(int i=1;i<=k;i++)
for(int j=1;j<=k;j++)
s[i][j]=s[i-1][j-1]+s[i-1][j]*j;
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
int x=n/l;
b[++tot]=x;
if(x<=sn) id1[x]=tot;
else id2[n/x]=tot;
g0[tot]=x-1;
gk[tot]=getsum(x)-1;
}
for(int i=1;i<=cnt;i++)
{
int pi2=p[i]*p[i];
for(int j=1;j<=tot&&pi2<=b[j];j++)
{
int v=get(b[j]/p[i]);
g0[j]=g0[j]-(g0[v]-(i-1));
gk[j]=gk[j]-powk[i]*(gk[v]-gpk[i-1]);
ans1[j]+=gk[v]-gpk[i-1];
}
}
}
ll query(int n)
{
return ans1[get(n)]+g0[get(n)];
}
}
int main()
{
n=read(),k=read();
s1::init();
s2::work();
ll ans=0,lst=0;
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ll now=s2::query(r);
ans+=(now-lst)*(s1::query(n/l)*2-1);
lst=now;
}
cout<<ans<<endl;
return 0;
}
/*
10 1
*/
/*
15 1
*/
/*
1000000 50
*/