题目大意：
\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)^klcm(i, j)[gcd(i,j)\in prime]\pmod {1e9 + 7}
其中 n\le {10}^{10},\ k\le 100
分析：
方法一：
\begin{aligned}
f(n,k)&=\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)^klcm(i, j)[gcd(i,j)\in prime]\pmod {1e9 + 7}\cr
&=\sum_{p}\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=p]\cdot p^k\cdot\frac{ij}{p}\cr &=\sum_{p}p^{k-1}\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)=p]\cdot ij\cr &=\sum_{p}p^{k+1}\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{n}{p}\right\rfloor}[gcd(i,j)=p]\cdot ij\cr &=\sum_{p}p^{k+1}\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}i^2\cdot\varphi(i)
\end{aligned}
后面的柿子显然是个杜教筛；前面的质数点值之和用 min25 筛维护，用拉格朗日插值法得到 k+1 次方的前缀和，并且对杜教筛部分整除分块。
方法二：对前面的柿子交换求和顺序：
\begin{aligned}
f(n,k)&=\sum_{p}p^{k+1}\sum_{i=1}^{\left\lfloor\frac{n}{p}\right\rfloor}i^2\cdot\varphi(i)\cr
&=\sum_{i=1}^{n}i^2\cdot\varphi(i)\sum_{p=1}^{\left\lfloor\frac{n}{i}\right\rfloor}p^{k+1}[p\in prime]
\end{aligned}
于是可以对 min25 筛部分做数论分块，用杜教筛维护前面的前缀和。
由于min25筛做完之后天然自带数论分块，因此十分好写：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int P = 1e9 + 7;
const int N = 5e6 + 50;
const int M = 1e5 + 50;
int inv2, inv6;
ll n; int k;
int ksm(int a, int b) {
int ret = 1;
for(; b; b >>= 1, a = 1LL * a * a % P) {
if(b & 1) ret = 1LL * ret * a % P;
}
return ret;
}
int ni(int a) {
return ksm(a, P - 2);
}
int ifac[M];
void init_inv() {
inv2 = ni(2); inv6 = ni(6);
ifac[0] = ifac[1] = 1;
for (int i = 2; i < M; i++) {
ifac[i] = 1LL * (P - P / i) * ifac[P % i] % P;
}
for (int i = 2; i < M; i++) {
ifac[i] = 1LL * ifac[i] * ifac[i - 1] % P;
}
}
int p[N], vis[N], tot;
int phi[N], Sf[N];
unordered_map<ll, int> Sumf;
void init_phi() {
phi[1] = 1;
for (int i = 2; i < N; i++) {
if (!vis[i]) {
p[++tot] = i;
phi[i] = i - 1;
}
for (int j = 1, x; j <= tot && (x = p[j] * i) < N; j++) {
vis[x] = 1;
if (i % p[j]) {
phi[x] = phi[i] * phi[p[j]];
}
else {
phi[x] = phi[i] * p[j]; break;
}
}
}
for (int i = 1; i < N; i++) {
phi[i] = 1LL * phi[i] * i % P * i % P;
Sf[i] = (Sf[i - 1] + phi[i]) % P;
}
}
int Scube(ll x) { // 三次方前缀和
x %= P;
int tmp = 1LL * x * (x + 1) / 2 % P;
return 1LL * tmp * tmp % P;
}
int Ssquare(ll x) { // 二次方前缀和
x %= P;
return 1LL * x * (x + 1) % P * (2 * x + 1) % P * inv6 % P;
}
int S(ll x) {
if (x < N) return Sf[x];
if (Sumf.count(x)) return Sumf[x];
int ret = Scube(x);
for (ll l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
int tmp = 1LL * S(x / l) * (Ssquare(r) - Ssquare(l - 1) + P) % P;
ret = (ret - tmp + P) % P;
}
return Sumf[x] = ret;
}
struct qm {
vector<int> y;
int a[N], b[N];
void init(int k) {
y.resize(1); y[0] = 0;
for (int i = 1; i <= k + 1; i++) {
y.push_back((y.back() + ksm(i, k)) % P);
}
}
int Lagrange(ll X) {
int n = (int)y.size() - 1, ans = 0; X %= P;
a[0] = (int)X; b[n + 1] = 1;
for (int i = 1; i <= n; i++) {
a[i] = 1LL * a[i - 1] * (X - i) % P;
}
for (int i = n; i >= 0; i--) {
b[i] = 1LL * b[i + 1] * (X - i) % P;
}
for (int i = 0; i <= n; i++)
(ans += 1LL * y[i] * (i == 0 ? 1 : a[i - 1]) % P * b[i + 1] % P
* ifac[i] % P * (((n - i)&1) ? -1 : 1) * ifac[n - i] % P) %= P;
return (ans + P) % P;
}
int solve(ll n) {
return Lagrange(n);
}
}I;
struct Min25 {
ll id1[M << 1], id2[M << 1], val[M << 1], g[M << 1];
ll n, lim, sp[N], powp[N];
void init() {
for (int i = 1; i <= tot; i++) {
powp[i] = ksm(p[i], k + 1);
sp[i] = (sp[i - 1] + powp[i]) % P;
}
}
void build(ll _n) {
n = _n; lim = sqrt(n);
int m = 0;
for(ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
val[++m] = n / l;
if(val[m] <= lim) id1[val[m]] = m;
else id2[n/val[m]] = m;
g[m] = (I.solve(val[m]) - 1 + P) % P;
}
for(int j = 1; j <= tot; j++) {
for(int i = 1; 1LL * p[j] * p[j] <= val[i]; i++) {
ll tmp = val[i] / p[j];
if(tmp <= lim) tmp = id1[tmp];
else tmp = id2[n/tmp];
g[i] -= 1LL * powp[j] * (g[tmp] - sp[j-1] + P) % P;
if(g[i] < 0) g[i] += P;
}
}
}
ll solve1(ll l) {
if (n / (n / l) <= lim) return g[id1[n / (n / l)]];
return g[id2[n / l]];
}
ll solve2(ll l) { /* n / l 以内质数点值之和 */
return n / l <= lim ? g[id1[n / l]] : g[id2[n / (n / l)]];
}
}Sn;
ll work1() {
ll ret = 0;
for (ll l = 1, r, lst = 0; l <= n; l = r + 1) {
r = n / (n / l);
ll tmp = Sn.solve1(l);
(ret += (tmp - lst + P) * S(n / l) % P) %= P;
lst = tmp;
}
return ret;
}
ll work2() {
ll ret = 0;
for (ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
(ret += Sn.solve2(l) * (S(r) - S(l - 1) + P) % P) %= P;
}
return ret;
}
int main(int argc, const char * argv[]) {
init_inv(); init_phi();
int T; cin >> T;
while (T--) {
cin >> n >> k;
I.init(k + 1);
Sn.init(); Sn.build(n);
cout << work1() << endl;
//cout << work2() << endl;
}
return 0;
}