# bzoj 4407 于神之怒加强版

1 2

3 3

20

### HINT

1<=N,M,K<=5000000,1<=T<=2000

### 题解

\large \begin{aligned} &\;\;\;\;\; \sum_{i = 1}^n\sum_{j = 1}^m (i,j)^k \\ &= \sum_d d^k \sum_{i = 1}^n\sum_{j = 1}^m [(i,j)=b] \\ &= \sum_d d^k \sum_{i = 1}^{n \over d}\sum_{j = 1}^{m \over d}[(i,j)=1] \\ &= \sum_d d^k \sum_{i = 1}^{n \over d}\sum_{j = 1}^{m \over d}\sum_{p|(i,j)} \mu(p) \\ &= \sum_d d^k \sum_p \mu(p) {n\over dp} {m\over dp} \\ &= \sum_d d^k \sum_p \mu(p) {n\over T} {m\over T} \\ &= \sum_{T = 1}^{min(n,m)} {n\over T} {m\over T} \sum_{d|T}d^k \mu({T\over d}) \end{aligned}

$$f(T)=\sum_{d|T}d^k\mu({T\over d})$$

#include <cstdio>
#include <algorithm>
using namespace std;
const int N = 5000000, mod = 1000000007;
bool p[N + 10];
int T, pr[N + 10], n, m;
long long k, f[N + 10], ans;
long long fast_pow(long long a, int p)
{
long long ans=1;
for (; p; p >>= 1, a = a * a % mod)
if (p & 1) ans = ans * a % mod;
return ans;
}
void init() {
f[1] = 1;
for (int i = 2; i <= N; ++i) {
if (!p[i])
pr[++pr[0]] = i, f[i] = (fast_pow(i, k) - 1 + mod) % mod;
for (int j = 1; j <= pr[0] && i * pr[j] <= N; ++j) {
p[i * pr[j]] = 1;
if (i % pr[j] == 0) {
f[i * pr[j]] = f[i] * fast_pow(pr[j], k) % mod;
break;
}
f[i * pr[j]] = f[i] * f[pr[j]] % mod;
}
}
for (int i = 2; i <= N; ++i) f[i] = (f[i] + f[i - 1]) % mod;
}
main() {
scanf("%d%d", &T, &k);
init();
while (T--) {
scanf("%d%d", &n, &m);
if (n > m) swap(n, m);
ans = 0;
for (int i = 1, las; i <= n; i = las + 1) {
las = min(n / (n / i), m / (m / i));
ans += (f[las] - f[i - 1]) * (n / i) % mod * (m / i) % mod;
ans = (ans % mod + mod) % mod;
}
printf("%lld\n", ans);
}
}