[BZOJ 4407] 于神之怒加强版

题目大意

给定 \(n\)\(m\)\(k\),求 \[ \sum_{i = 1}^{n} \sum_{j = 1}^{m} gcd(i, \; j)^k \bmod 1,000,000,007 \] 多组询问(\(k\)都相同)。

\(1 \leqslant T \leqslant 2,000\)

\(1 \leqslant n, \; m, \; k \leqslant 5,000,000\)

题目链接

BZOJ 4407

题解

\[ \begin{align} &\sum_{i = 1}^{n} \sum_{j = 1}^{m} gcd(i, \; j)^k \\ =& \sum_{d = 1}^{n} d^k \sum_{i = 1}^{n} \sum_{j = 1}^{m} [gcd(i, \; j) = d] \\ = &\sum_{d = 1}^{n} d^k \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(i) \lfloor \frac{n}{di} \rfloor \lfloor \frac{m}{di} \rfloor \\ \end{align} \]

此时,如果每次枚举 \(d\),每次询问复杂度是 \(O(n \log k)\) 的(\(O(\log k)\) 来自快速幂,\(O(n)\) 来自 \(\sum_{i = 1}^{n} \sqrt{\frac{n}{i}}\),画了个图像发现渐近线是 \(y = 2x\)),会 TLE。不过,似乎有一些卡着时限 AC 的人是这个做法(猜的)。

改变枚举对象为 \(T = di\)\[ \begin{align} &\sum_{d = 1}^{n} d^k \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(i) \lfloor \frac{n}{di} \rfloor \lfloor \frac{m}{di} \rfloor \\ = &\sum_{T = 1}^{n} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \sum_{d | T} d^k \mu(\lfloor \frac{T}{d} \rfloor) \end{align} \] 记后面的和式为 \(f(T)\),发现它有积性(积性函数的狄利克雷卷积为积性函数)。

\(T = \prod p_i^{x_i}\),有: \[ \begin{align} f(T) = &\sum_{d | T} d^k \mu(\frac{T}{d}) \\ = &\prod f(p_i^{x_i}) \\ = &\prod p_i^{x_I} \mu(1) + p_i^{k(x_i - 1)} \mu(p_i) \\ = &\prod p_i^{k(x_i - 1)} (p_i^k - 1) \end{align} \] 可以得到,当 \(T = T' \times p^2\) 时,有 \(f(T) = f(T') \times (f(p) + 1)\)

于是乎,我们就可以用线性筛求 \(f(x)\) 函数了。询问时,分块计算。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <cstdio>
#include <algorithm>
const int MAXN = 5000005;
const int MOD = 1000000007;
long long f[MAXN];
int prime[MAXN], primeCnt;
bool notPrime[MAXN];
long long pow(long long a, long long n) {
long long res = 1;
for (; n; n >>= 1, a = a * a % MOD) if (n & 1) res = res * a % MOD;
return res;
}
int k;
void linearShaker() {
notPrime[0] = notPrime[1] = true;
f[1] = 1;
for (int i = 2; i < MAXN; i++) {
if (!notPrime[i]) {
prime[++primeCnt] = i;
f[i] = (pow(i, k) - 1 + MOD) % MOD;
}
for (int j = 1; j <= primeCnt && i * prime[j] < MAXN; j++) {
notPrime[i * prime[j]] = true;
if (i % prime[j] == 0) {
f[i * prime[j]] = f[i] * (f[prime[j]] + 1) % MOD;
break;
} else f[i * prime[j]] = f[i] * f[prime[j]] % MOD;
}
}
for (int i = 2; i < MAXN; i++) f[i] = (f[i] + f[i - 1]) % MOD;
}
long long calc(int n, int m) {
if (n > m) std::swap(n, m);
long long res = 0;
for (int i = 1, last; i <= n; i = last + 1) {
last = std::min(n / (n / i), m / (m / i));
(res += (f[last] - f[i - 1]) * (n / i) % MOD * (m / i) % MOD) %= MOD;
}
return (res % MOD + MOD) % MOD;
}
int main() {
int T;
scanf("%d %d", &T, &k);
linearShaker();
while (T--) {
int n, m;
scanf("%d %d", &n, &m);
printf("%lld\n", calc(n, m));
}
return 0;
}