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}^nsum_{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}
floor}sum_{j=1}^{lfloor frac{n}{d}
floor}[gcd(i,j)=1]
]
[=sum_{d=1}^{n}f^k(d)(2sum_{i=1}^{lfloorfrac{n}{d}
floor}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'(lfloorfrac{x}{p_i^{e}}
floor,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'(lfloorfrac{x}{p_i^{e}}
floor,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
", ans);
return 0;
}