\(\sum_{i=1}^n\sum_{j=1}^n\mathrm{sgcd}(i,j)^k=\sum_{p=1}^ns(p)^k\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=p]=\sum_{p=1}^ns(p)^k(-1+2\sum_{i=1}^{n/p}\varphi(i))\)
由于 \(n\) 的范围是 \(10^9\) ，对于后面的我们最多只有根号种取值，根据套路，可以杜教筛/Min_25筛一波。
至于前面的东西，我们可以考虑Min_25筛的过程：
Min_25筛我们设 \(g(n,j)\) 为2~n内素数或最小质因子\(\ge p_j\)的数的\(k\)次方的和
考虑转移：
\(g(n,j)=g(n,j-1)-p_j^k(g(a/p_j,j-1)-g(p_{j-1},j-1)) (n\ge p_j^2)\)
我们发现后面减去的 \(p_j^k(g(a/p_j,j-1)-g(p_{j-1},j-1))\) 就是最小质因子恰为 \(p_j\) 的合数的k次方的和
那么 \(g(a/p_j,j-1)-g(p_{j-1},j-1)\) 就是[1, n]内这最小质因子为 \(p_j\) 的合数的 \(s(x)^k\) 的和。
我们对于每个n/x向下取整开一个数，记录每个 j 的这个值，然后差分一下就是区间的答案了。（注意加上质数的贡献，也就是区间质数个数）
预处理Min_25筛的数组时，由于模数是合数，所以不能拉格朗日插值，要用第二类斯特林数：
\(\sum_{i=1}^ni^k=\sum_{i=1}^n\sum_{j=0}^k\{^k_j\}{i\choose j}j!=\sum_{j=0}^k\{^k_j\}j!\sum_{i=1}^n{i\choose j}=\sum_{j=0}^k\{^k_j\}j!{n+1\choose j+1}=\sum_{j=0}^k\{^k_j\}\frac{n+1^{\underline{j+1}}}{j+1}\)
phi用Min_25筛的代码：(观察txc巨佬的)
#include <cstdio>
#include <cmath>
#define int unsigned
using namespace std;
int n, m, k, sqt, tot;
int a[233333], g0[233333], g1[233333], gk[233333], gphi[233333], ans[233333], prime[233333], s[233][233];
int getid(int x) { return x <= sqt ? x : m - n / x + 1; }
int qsum(int n)
{
int ans = 0;
for (int i = 0; i <= k; i++)
{
int tmp = 1;
for (int j = 0; j <= i; j++) //(n - j + 1)
if ((n - j + 1) % (i + 1) == 0) tmp *= ((n - j + 1) / (i + 1));
else tmp *= (n - j + 1);
ans += tmp * s[k][i];
}
return ans;
}
signed main()
{
scanf("%u%u", &n, &k), sqt = sqrt(n);
s[0][0] = 1;
for (int i = 1; i <= k; i++) for (int j = 1; j <= i; j++) s[i][j] = s[i - 1][j] * j + s[i - 1][j - 1];
for (int i = 1; i <= n; i = a[m] + 1)
{
a[++m] = n / (n / i);
gk[m] = qsum(a[m]) - 1;
g1[m] = a[m] * (long long)(a[m] + 1) / 2 - 1;
g0[m] = a[m] - 1;
}
for (int i = 1; i <= sqt; i++) if (g0[i] != g0[i - 1])
{
prime[++tot] = i;
int tmp = 1;
for (int t = 1; t <= k; t++) tmp *= i;
for (int j = m; a[j] >= i * i; j--)
{
int id = getid(a[j] / i);
gk[j] -= tmp * (gk[id] - gk[i - 1]);
ans[j] += gk[id] - gk[i - 1];
g1[j] -= i * (g1[id] - g1[i - 1]);
g0[j] -= g0[id] - g0[i - 1];
}
}
for (int i = 1; i <= m; i++) gphi[i] = (g1[i] -= g0[i]);
for (int i = tot; i >= 1; i--)
for (int j = m; a[j] >= prime[i] * prime[i]; j--)
for (int x = prime[i], f = x - 1; x * (long long)prime[i] <= a[j]; x *= prime[i], f *= prime[i])
gphi[j] += f * (gphi[getid(a[j] / x)] - g1[prime[i]] + prime[i]);
int res = 0;
for (int i = 2; i <= m; i++) res += (ans[i] - ans[i - 1] + g0[i] - g0[i - 1]) * (2 * gphi[m - i + 1] + 1);
printf("%u\n", res);
return 0;
}
phi用杜教筛写的代码：
#include <cmath>
#include <cstdio>
#include <cstring>
using namespace std;
#define int unsigned
int n, k, sqt;
int a[233333], g0[233333], gk[233333], prime[233333], ans[233333], tot, m;
int s[60][60];
int phi[233333], mem[233333];
bool vis[233333];
int qid(int x) { return x <= sqt ? x : m - n / x + 1; }
int qsum(int x, int k)
{
int ans = 0;
for (int i = 0; i <= k; i++)
{
int tmp = 1;
for (int j = 0; j <= i; j++)
{
if ((x + 1 - j) % (i + 1) == 0) tmp *= (x + 1 - j) / (i + 1);
else tmp *= x + 1 - j;
}
ans += tmp * s[k][i];
}
return ans;
}
int chuphi(int id)
{
if (a[id] <= sqt) { return phi[a[id]]; }
if (mem[id] != 0xffffffff) return mem[id];
int n = a[id], ans = n * (n + 1) / 2;
for (int i = 2; i <= n; i = n / (n / i) + 1)
ans -= chuphi(qid(n / i)) * (n / (n / i) - i + 1);
return mem[id] = ans;
}
signed main()
{
memset(mem, 0xff, sizeof(mem)); phi[1] = 1;
scanf("%u%u", &n, &k), sqt = sqrt(n);
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] + j * s[i - 1][j];
for (int i = 1; i <= n; i = a[m] + 1)
{
a[++m] = n / (n / i);
g0[m] = a[m] - 1;
gk[m] = qsum(a[m], k) - 1;
}
for (int i = 1; i <= sqt; i++)
{
if (g0[i] != g0[i - 1])
{
prime[++tot] = i, phi[i] = i - 1;
int tmp = 1; for (int j = 1; j <= k; j++) tmp *= i;
for (int j = m; a[j] >= i * i; j--)
{
int id = qid(a[j] / i);
g0[j] -= g0[id] - g0[i - 1];
gk[j] -= tmp * (gk[id] - gk[i - 1]);
ans[j] += gk[id] - gk[i - 1];
}
}
for (int j = 1; j <= tot && i * prime[j] <= sqt; j++)
{
if (i % prime[j] == 0) { phi[i * prime[j]] = phi[i] * prime[j]; break; }
phi[i * prime[j]] = phi[i] * (prime[j] - 1);
}
phi[i] += phi[i - 1];
}
int res = 0;
for (int i = 1; i <= m; i++)
{
res += (ans[i] - ans[i - 1] + g0[i] - g0[i - 1]) * (2 * chuphi(m - i + 1) - 1);
}
printf("%u\n", res);
return 0;
}