Sol
设 \(f(d)\) 表示 \(d\) 所有约数中第二大的，\(low_d\) 表示 \(d\) 的最小质因子
\[f(d)=\frac{d}{low_d}
\]
\]
那么
\[\sum_{i=1}^{n}\sum_{j=1}^{n}sgcd^k(i,j)
\]
\]
\[=\sum_{i=1}^n\sum_{j=1}^{n}f^k(gcd(i,j))
\]
\]
\[=\sum_{d=1}^{n}f^k(d)\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=d]
\]
\]
\[=\sum_{d=1}^{n}f^k(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor \frac{n}{d}\rfloor}[gcd(i,j)=1]
\]
\]
\[=\sum_{d=1}^{n}f^k(d)(2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\varphi(i)-1)
\]
\]
可以数论分块，后面的 \(\varphi\) 直接杜教筛
考虑计算
\[\sum_{d=2}^{n}f^k(d)=\sum_{d=2}^{n}\frac{d}{low_d}^k
\]
\]
设 \(p_j\) 表示第 \(j\) 个质数
\(g(x,i)\) 表示 \(2\) 到 \(x\) 之间最小质因子大于等于 \(p_i\) 的或者质数的 \(f\) 的 \(k\) 次方和
\(g'(x,i)\) 表示 \(2\) 到 \(x\) 之间最小质因子大于等于 \(p_i\) 的或者质数的 \(k\) 次方和
设 \(s(i)\) 表示小于等于 \(p_i\) 的质数的 \(k\) 次方和
那么就是要求 \(g(n, 1)\)
\(s\) 直接 \(min25\) 筛
\[g(x,i)=g(x,i+1)+\sum_{e=1}^{p_i^{e+1}\le x}p_i^{k(e-1)}(g'(\lfloor\frac{x}{p_i^{e}}\rfloor,i+1)-s(i)+p_i^{k})
\]
\]
\[g'(x,i)=g'(x,i+1)+\sum_{e=1}^{p_i^{e+1}\le x}p_i^{ke}(g'(\lfloor\frac{x}{p_i^{e}}\rfloor,i+1)-s(i)+p_i^{k})
\]
\]
这里要用到自然幂数和求和，用第二类斯特林数就好了
# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned int uint;
inline uint Pow(uint x, int y) {
register uint ret = 1;
for (; y; y >>= 1, x = x * x) if (y & 1) ret = ret * x;
return ret;
}
const int maxn(1e6 + 5);
int pr[maxn], tot, k, d, id1[maxn], id2[maxn], cnt, phi[maxn];
uint n, s[55][55], f[maxn], g[maxn], val[maxn], sk[maxn];
uint s1[maxn], s2[maxn], sphi[maxn], sp[maxn];
bitset <maxn> ispr;
inline void Sieve(int mx) {
register int i, j;
for (phi[1] = 1, ispr[1] = 1, i = 2; i <= mx; ++i) {
if (!ispr[i]) pr[++tot] = i, sk[tot] = sk[tot - 1] + Pow(i, k), phi[i] = i - 1;
for (j = 1; j <= tot && pr[j] * i <= mx; ++j) {
ispr[i * pr[j]] = 1;
if (i % pr[j]) phi[i * pr[j]] = phi[i] * (pr[j] - 1);
else {
phi[i * pr[j]] = phi[i] * pr[j];
break;
}
}
}
for (i = 1; i <= mx; ++i) sphi[i] = sphi[i - 1] + phi[i];
}
# define ID(x) ((x) <= d ? id1[x] : id2[n / (x)])
inline uint Sum(uint x) {
register uint i, j, v = 0, t, r, tmp = k <= x ? k : x;
for (i = 1; i <= tmp; ++i) {
t = i + 1, r = s[k][i];
for (j = x - i + 1; j <= x + 1; ++j)
if (t > 1 && j % t == 0) r *= j / t, t = 1;
else r = r * j;
v += r;
}
return v;
}
uint Sumphi(uint x) {
if (x <= d) return sphi[x];
if (sp[ID(x)]) return sp[ID(x)];
register uint ans = (x & 1) ? ((x + 1) >> 1) * x : (x >> 1) * (x + 1), i, j;
for (i = 2; i <= x; i = j + 1) j = x / (x / i), ans -= Sumphi(x / i) * (j - i + 1);
return sp[ID(x)] = ans;
}
int main() {
register uint i, j, e, ans = 0, lst = 0, cur, r, v, tmp, now;
scanf("%u%d", &n, &k), Sieve(d = sqrt(n));
for (i = 1; i <= 50; ++i)
for (s[i][1] = 1, j = 2; j <= i; ++j)
s[i][j] = s[i - 1][j] * j + s[i - 1][j - 1];
for (i = 1; i <= n; i = j + 1) {
val[++cnt] = n / i, j = n / (n / i);
val[cnt] <= d ? id1[val[cnt]] = cnt : id2[n / val[cnt]] = cnt;
f[cnt] = val[cnt] - 1, g[cnt] = Sum(val[cnt]) - 1;
}
for (i = 1; i <= tot && pr[i] * pr[i] <= n; ++i)
for (j = 1; j <= cnt && pr[i] * pr[i] <= val[j]; ++j) {
f[j] -= f[ID(val[j] / pr[i])] - i + 1;
g[j] -= (sk[i] - sk[i - 1]) * (g[ID(val[j] / pr[i])] - sk[i - 1]);
}
for (i = 1; i <= cnt; ++i) s1[i] = f[i], s2[i] = g[i];
for (r = 1; r <= tot && pr[r] * pr[r] <= n; ++r);
for (i = r - 1; i; --i)
for (tmp = sk[i] - sk[i - 1], j = 1; j <= cnt && pr[i] * pr[i] <= val[j]; ++j)
for (cur = e = 1, v = pr[i]; v <= val[j] / pr[i]; ++e, v *= pr[i], cur *= tmp)
now = (s2[ID(val[j] / v)] - sk[i] + tmp) * cur, s1[j] += now, s2[j] += tmp * now;
for (i = 1; i <= n; i = j + 1) {
j = n / (n / i), cur = s1[ID(j)];
ans += (cur - lst) * (Sumphi(n / i) * 2 - 1), lst = cur;
}
printf("%u\n", ans);
return 0;
}