「BZOJ 4916」神犇和蒟蒻-ExtendedEratosthenesSieve

求 $A = \sum\limits_{i = 1} ^ n \mu(i ^ 2), B = \sum\limits_{i = 1} ^ n \varphi(i ^ 2)$

链接

BZOJ 4916

题解

显然 $A = 1$,然后我们只用求 $B$。

考虑扩展埃拉托色尼筛法,这里直接不推式子暴力筛

令 $f(n) = \varphi(n ^ 2)$,则 $f(1) = 1$
当 $n = p$,$p$ 为质数时,$f(p) = (p - 1)p = p ^ 2 - p$
当 $n = p ^ e$ 时,$f(p ^ e) = (p - 1)p ^ {e - 1}$

于是我们需要 $p ^ 0, p ^ 1, p ^ 2$ 的前缀和,$f(p)$ 可以通过 $p ^ 2 - p$ 计算出,枚举指数时先算出 $p ^ 2$,维护 $p ^ e$ 的答案,初始为 $p ^ 2 - p$,然后每次乘以 $p ^ 2$ 就可以保证复杂度,然后就做完了。

代码

显然可以优化,但是这是不推式子直接暴力筛,然而跑得飞快

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
101
102
103
104
105
106
107
108
109
110
/**
* Copyright (c) 2017-2018, xehoth
* All rights reserved.
* 「BZOJ 4916」神犇和蒟蒻 24-01-2018
* Extended Eratosthenes Sieve
* @author xehoth
*/
#include <bits/stdc++.h>

int d, M, n;
const int MOD = 1e9 + 7;
const int INV6 = 166666668;

std::vector<int> pre[3], suc[3], primes;

inline int add(const int x, const int v) {
return x + v >= MOD ? x + v - MOD : x + v;
}

inline int dec(const int x, const int v) {
return x - v < 0 ? x - v + MOD : x - v;
}

int rec(int res, int last, int mul) {
int t = dec((res > M ? suc[2][n / res] : pre[2][res]),
pre[2][primes[last] - 1]);
int ret = (unsigned long long)t * mul % MOD;
for (int i = last, p; i < (int)primes.size(); i++) {
p = primes[i];
if (p * p > res) break;
const int p2 = (unsigned long long)p * p % MOD;
for (int q = p, nres = res, nmul = p * (p - 1ull) % MOD * mul % MOD;
(unsigned long long)p * q <= res; q *= p) {
ret = add(ret, rec(nres /= p, i + 1, nmul));
nmul = (unsigned long long)nmul * p2 % MOD;
ret = add(ret, nmul);
}
}
return ret;
}

inline int solve(const int n) {
M = sqrt(n);
primes.reserve(M);
pre[0].resize(M + 1);
pre[1].resize(M + 1);
pre[2].resize(M + 1);
suc[0].resize(M + 1);
suc[1].resize(M + 1);
suc[2].resize(M + 1);
for (int i = 1, t; i <= M; i++) {
pre[0][i] = i - 1;
pre[1][i] = (i * (i + 1ull) / 2 - 1 + MOD) % MOD;
pre[2][i] =
(i * (i + 1ull) % MOD * (2 * i + 1ull) % MOD * INV6 + MOD - 1) %
MOD;
t = n / i;
suc[0][i] = t - 1;
suc[1][i] = (t * (t + 1ull) / 2 - 1 + MOD) % MOD;
suc[2][i] =
(t * (t + 1ull) % MOD * (2 * t + 1ull) % MOD * INV6 + MOD - 1) %
MOD;
}
for (int p = 2, end; p <= M; p++) {
if (pre[0][p] == pre[0][p - 1]) continue;
primes.push_back(p);
const int q = p * p, m = n / p;
const int pcnt = pre[0][p - 1], psum = pre[1][p - 1],
psum2 = pre[2][p - 1], p2 = (unsigned long long)p * p % MOD;
end = std::min(M, n / q);
for (int i = 1, w = M / p; i <= w; i++) {
suc[0][i] = dec(suc[0][i], dec(suc[0][i * p], pcnt));
suc[1][i] =
dec(suc[1][i],
(unsigned long long)dec(suc[1][i * p], psum) * p % MOD);
suc[2][i] =
dec(suc[2][i],
(unsigned long long)dec(suc[2][i * p], psum2) * p2 % MOD);
}
for (int i = M / p + 1; i <= end; i++) {
suc[0][i] = dec(suc[0][i], dec(pre[0][m / i], pcnt));
suc[1][i] =
dec(suc[1][i],
(unsigned long long)dec(pre[1][m / i], psum) * p % MOD);
suc[2][i] =
dec(suc[2][i],
(unsigned long long)dec(pre[2][m / i], psum2) * p2 % MOD);
}
for (int i = M; i >= q; i--) {
pre[0][i] = dec(pre[0][i], dec(pre[0][i / p], pcnt));
pre[1][i] =
dec(pre[1][i],
(unsigned long long)dec(pre[1][i / p], psum) * p % MOD);
pre[2][i] =
dec(pre[2][i],
(unsigned long long)dec(pre[2][i / p], psum2) * p2 % MOD);
}
}
primes.push_back(M + 1);
for (int i = 1; i <= M; i++) {
pre[2][i] = dec(pre[2][i], pre[1][i]);
suc[2][i] = dec(suc[2][i], suc[1][i]);
}
return n > 1 ? 1 + rec(n, 0, 1) : 1;
}

int main() {
std::cin >> n;
std::cout << "1\n" << solve(n);
}
#

Comments

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×