(sum_{i=1}^nsum_{j=1}^nmathrm{sgcd}(i,j)^k=sum_{p=1}^ns(p)^ksum_{i=1}^nsum_{j=1}^n[gcd(i,j)=p]=sum_{p=1}^ns(p)^k(-1+2sum_{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)) (nge 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}^nsum_{j=0}^k{^k_j}{ichoose j}j!=sum_{j=0}^k{^k_j}j!sum_{i=1}^n{ichoose j}=sum_{j=0}^k{^k_j}j!{n+1choose 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
", 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
", res);
return 0;
}