根据Crash的数字表格，很容易可以将式子化简为
$$
\begin{aligned} Ans &= \sum\limits_{i = 1}^n \sum\limits_{j = 1} ij(i, j) \\ &= \sum\limits_{d = 1}^n d^3 \sum\limits_{k = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} \mu(k) k^2 \left( \sum\limits_{i = 1}^{\left\lfloor\frac{n}{kd}\right\rfloor} i \right)^2 \end{aligned}
$$
感觉 $d, k$ 放在一起式子无法继续化简，主要是有 $kd$ 存在，故令 $T = kd$ ，则有
$$
Ans = \sum\limits_{T = 1}^n \left( \sum\limits_{i = 1}^{\left\lfloor\frac{n}{T}\right\rfloor} i \right)^2 T^2 \sum\limits_{d | T} d \mu(\frac{T}{d})
$$
那么考虑整除分块，现在需要处理的是后半部分
通过观察（看题解）可以发现， $\sum\limits_{d | T} d \mu(\frac{T}{d})$ 可以看成 $\mu * id$ ，故可替换成 $\phi(T)$ ，则有
$$
Ans = \sum\limits_{T = 1}^n \left( \sum\limits_{i = 1}^{\left\lfloor\frac{n}{T}\right\rfloor} i \right)^2 T^2 \phi(T)
$$
现在考虑将 $T^2 \phi(T)$ 部分用杜教筛解决
$$
h(n) = \sum\limits_{d | n} d^2 \phi(d) g(\frac{n}{d})
$$
为了消除 $d^2$ ，令 $g(\frac{n}{d}) = n^2$ ，则有
$$
h(n) = n^3
$$
故得
$$
S(n) = \sum\limits_{i = 1}^n i^3 - \sum\limits_{d = 2}^n d^2 S(\left\lfloor\frac{n}{d}\right\rfloor)
$$
又（通过看题解）有一个知识点
$$
\sum\limits_{i = 1}^n i^3 = \left( \sum\limits_{i = 1}^n i \right)
$$
那么就可以直接杜教筛了
至于复杂度，将最外围的整除分块与杜教筛看为一体，故复杂度为 $O (n^{\frac{2}{3}})$
代码
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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
| #include <iostream>
#include <cstdio>
#include <cstring>
#include <tr1/unordered_map>
using namespace std;
typedef long long LL;
const int MAXN = 6e06 + 10;
LL MOD, N;
LL power (LL x, LL p) {
LL cnt = 1;
while (p) {
if (p & 1)
cnt = cnt * x % MOD;
x = x * x % MOD;
p >>= 1;
}
return cnt;
}
LL inv2, inv6;
int prime[MAXN];
int vis[MAXN]= {0};
int pcnt = 0;
LL phi[MAXN]= {0};
LL sumphi[MAXN]= {0};
int MAX = 6e06;
void linear_sieve () {
phi[1] = 1;
for (int i = 2; i <= MAX; i ++) {
if (! vis[i]) {
prime[++ pcnt] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= pcnt && i * prime[j] <= MAX; j ++) {
vis[i * prime[j]] = 1;
if (! (i % prime[j])) {
phi[i * prime[j]] = phi[i] * 1ll * prime[j] % MOD;
break;
}
phi[i * prime[j]] = phi[i] * 1ll * (prime[j] - 1) % MOD;
}
}
for (int i = 1; i <= MAX; i ++)
sumphi[i] = (sumphi[i - 1] + 1ll * i % MOD * 1ll * i % MOD * phi[i] % MOD) % MOD;
}
tr1::unordered_map<LL, LL> maphi;
inline LL sqr (LL x) {
return x * x % MOD;
}
inline LL eqm (LL n) {
return (n + 1) % MOD * (n % MOD) % MOD * (2 * n % MOD + 1) % MOD * inv6 % MOD;
}
inline LL oseqm (LL n) {
return n % MOD * ((n + 1) % MOD) % MOD * inv2 % MOD;
}
LL phi_sieve (LL n) {
if (n <= MAX)
return sumphi[n];
if (maphi[n])
return maphi[n];
LL total = sqr (oseqm (n));
for (LL l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
total = (total - (eqm (r) - eqm (l - 1) + MOD) % MOD * phi_sieve (n / l) % MOD + MOD) % MOD;
}
return maphi[n] = total;
}
LL Solve () {
LL ans = 0;
for (LL l = 1, r; l <= N; l = r + 1) {
r = N / (N / l);
ans = (ans + sqr (oseqm (N / l)) * ((phi_sieve (r) - phi_sieve (l - 1) + MOD) % MOD) % MOD) % MOD;
}
return ans;
}
int main () {
scanf ("%lld%lld", & MOD, & N);
inv2 = power (2ll, MOD - 2), inv6 = power (6ll, MOD - 2);
MAX = (int) min (1ll * MAX, N), linear_sieve ();
LL ans = Solve ();
cout << ans << endl;
return 0;
}
|
