倍数ゼータ・メビウス変換のメモ
Last updated
Apr 19, 2022
約数ゼータ・メビウス変換に続いて、倍数もメモします。
# 倍数ゼータ・メビウス変換
関数$f(n)$に対する倍数ゼータ変換は以下の定義です。$n | m$は「$m$は$n$の倍数」という意味で、$m$はmultipleのmです。
$$
F(n) = \sum_{ n | m } f(m)
$$
倍数ゼータ変換の逆操作が倍数メビウス変換です。
# 倍数ゼータ・メビウス変換の実装
高速ゼータ変換の約数版 O(N log(log(N))) - noshi91のメモを参考にすると早い実装が得られます。
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
| template <class T>
void multiple_transform(vector<T> &a) {
int n = a.size();
vector<bool> sieve(n, true);
for (int p = 2; p < n; ++p) {
if (sieve[p]) {
for (int k = (n - 1) / p; k != 0; --k) {
sieve[k * p] = false;
a[k] += a[k * p];
}
}
}
for (int i = 0; ++i != n;) {
a[i] += a[0];
}
}
template <class T>
void inverse_multiple_transform(vector<T> &a) {
int n = a.size();
vector<bool> sieve(n, true);
for (int i = 0; ++i != n;) {
a[i] -= a[0];
}
for (int p = 2; p < n; ++p) {
if (sieve[p]) {
for (int k = 1; k * p < n; ++k) {
sieve[k * p] = false;
a[k] -= a[k * p];
}
}
}
}
|
# gcdとの関連
$f(n)$の定義が「$\gcd$が$n$となる集合についての対応する値の総和」、つまり下のような式の形をしていたとしましょう。
$$
f(n) = \sum_{n = \gcd(S)} g(S)
$$
これの倍数ゼータ変換$F(n)$を考えると
$$
\begin{align}
F(n) = & \sum_{n | m} \sum_{m = \gcd(S)} g(S) \\\ = & \sum_{n | \gcd(S)} g(S) \\\ = & \sum_{n | S_0, n | S_1, \cdots} g(S) \\\ \end{align}
$$
となって、すべての要素が$n$を約数にもつような集合$S$についての計算になります。
これがかなりいい性質で、数$A$の約数列挙は$O(\sqrt A)$でできるので、
$F(n)$の$n$それぞれについての計算が高速にできる場合があります。
# AGC038 C - LCMs
AGC038 C - LCMs
$$
\begin{align}
& \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} \operatorname{lcm}(A_i, A_j) \\\ = & \sum_{i = 0}^{N - 1} \sum_{j = 0}^{i - 1} \frac{A_i A_j}{\gcd(A_i, A_j)} \\\ \end{align}
$$
ここで、$g = \gcd(A_i, A_j)$となる$i, j \\ (i < j)$それぞれについてまとめて計算することにします。
$$
\begin{align}
= & \sum_{g} \sum_{g = gcd(A_i, A_j), i < j} \frac{A_i A_j}{g} \\\ = & \sum_{g} \frac{1}{g} \sum_{g = gcd(A_i, A_j), i < j} A_i A_j \\\ \end{align}
$$
$f(g) = \sum_{g = gcd(A_i, A_j), i < j} A_i A_j$として、これを倍数ゼータ変換をした$F(n)$を考えましょう。
$$
\begin{align}
F(n) & = \sum_{n | m} f(m) \\\ & = \sum_{n | m} \sum_{m = \gcd(A_i, A_j), i < j} A_i A_j \\\ & = \sum_{n | \gcd(A_i, A_j), i < j} A_i A_j \\\ & = \sum_{n | A_i, n | A_j, i < j} A_i A_j \\\ & = \frac{1}{2} ((\sum_{n | A_i} A_i)^2 - \sum_{n | A_i} A_i^2)
\end{align}
$$
$\sum_{n | A_i} A_i$と$\sum_{n | A_i} A_i^2$は、各$i$について$A_i$の約数を列挙すれば計算できて、全体で$O(N \sqrt A)$で計算できます。あとは計算した$F(n)$を倍数メビウス変換にかければ$f(n)$が得られます。
提出コード
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
| constexpr i64 MAX = 1e6 + 1;
int main() {
i64 N;
cin.tie(nullptr);
std::ios::sync_with_stdio(false);
cin >> N;
vector<i64> A(N);
rep(i,0,N) cin >> A[i];
vector<fp> S(MAX);
vector<fp> ans(MAX);
rep(i,0,N) {
for(i64 d = 1; d * d <= A[i]; d++) {
if(A[i] % d == 0) {
S[d] += fp(A[i]);
ans[d] -= A[i] * A[i];
if(A[i] / d != d) {
S[A[i] / d] += fp(A[i]);
ans[A[i] / d] -= A[i] * A[i];
}
}
}
}
fp i2 = fp(2).inv();
rep(g,1,MAX) {
ans[g] += S[g] * S[g];
ans[g] *= i2;
}
inverse_multiple_transform(ans);
fp res;
rep(g,1,MAX) {
res += fp(g).inv() * ans[g];
}
cout << res << endl;
}
|
# ちなみに…
$\sum_{n | A_i} A_i$は倍数ゼータ変換の形をしているので、この計算を倍数ゼータ変換して求めることで高速に計算できます。これがいわゆる$\gcd$畳み込みです。
提出コード
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
| constexpr i64 MAX = 1e6 + 1;
int main() {
i64 N;
cin.tie(nullptr);
std::ios::sync_with_stdio(false);
cin >> N;
vector<i64> A(N);
rep(i,0,N) cin >> A[i];
vector<fp> S(MAX);
rep(i,0,N) {
S[A[i]] += A[i];
}
multiple_transform(S);
rep(g,1,MAX) {
S[g] = S[g] * S[g];
}
inverse_multiple_transform(S);
fp res;
rep(g,1,MAX) {
res += fp(g).inv() * S[g];
}
rep(i,0,N) {
res -= A[i];
}
cout << res / fp(2) << endl;
}
|
# ちなみに…2
約数を列挙する部分でosa_k法を用いるとめちゃ速くなります。
提出コード
osa_kの実装は
えびちゃんのエラトステネスの篩に基づく高速な素因数分解 - Qiitaを参考にしました。
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
| struct osa_k {
using int_type = int;
std::vector<int_type> min_fact;
// O(NlogN)
static std::vector<int_type> min_prime_factor(int n) {
std::vector<int_type> res(n);
std::iota(std::begin(res), std::end(res), 0);
for(int i = 2; i * i < n; i++) {
if(res[i] < i) continue;
for(int j = i * i; j < n; j += i) {
if(res[j] == j) res[j] = i;
}
}
return res;
}
void build(int n) {
min_fact = min_prime_factor(n);
}
// O(logN)
std::vector<std::pair<int_type, int>> prime_factors(int n) const {
std::vector<std::pair<int_type, int>> res;
while(n > 1) {
if(res.empty() || res.back().first != min_fact[n]) {
res.push_back({ min_fact[n], 0 });
}
res.back().second++;
n /= min_fact[n];
}
return res;
}
// The divisors are not sorted
// O(logN + |divisors|)
template<class F>
void enumerate_divisors(int n, F f) const {
std::vector<std::pair<int_type, int>> prime_facts = prime_factors(n);
if(prime_facts.empty()) {
f(1);
return;
}
std::vector<int> cnt(prime_facts.size());
std::vector<int> acc(prime_facts.size(), 1);
while(true){
f(acc.front());
int i = 0;
for(; i < prime_facts.size(); i++) {
if((cnt[i]++) == prime_facts[i].second) {
cnt[i] = 0;
}
else {
acc[i] *= prime_facts[i].first;
break;
}
}
if(i == prime_facts.size()) {
break;
}
while(i --> 0) {
acc[i] = acc[i + 1];
}
}
}
};
...
int main() {
i64 N;
cin.tie(nullptr);
std::ios::sync_with_stdio(false);
cin >> N;
vector<i64> A(N);
rep(i,0,N) cin >> A[i];
osa_k osa;
osa.build(MAX);
vector<fp> S(MAX);
vector<fp> ans(MAX);
rep(i,0,N) {
osa.enumerate_divisors(A[i], [&](int d) {
S[d] += fp(A[i]);
ans[d] -= fp(A[i]) * fp(A[i]);
});
}
fp i2 = fp(2).inv();
rep(g,1,MAX) {
ans[g] += S[g] * S[g];
ans[g] *= i2;
}
inverse_multiple_transform(ans);
fp res;
rep(g,1,MAX) {
res += fp(g).inv() * ans[g];
}
cout << res << endl;
}
|
# ABC248 G - GCD cost on the tree
ABC248 G - GCD cost on the tree
$P(i, j)$を$(i, j)$パス上の頂点の$A$の値の集合、$d(i, j)$を$(i, j)$パスの長さとすると、
$$
\begin{align}
\sum_{i < j} (d(i, j) + 1) \gcd(P(i, j))
\end{align}
$$
が答えです。ここから変形していきます。
$$
\begin{align}
= & \sum_{g} \sum_{g = \gcd(P(i, j)), i < j} (d(i, j) + 1) g \\\ = & \sum_{g} g \sum_{g = \gcd(P(i, j)), i < j} (d(i, j) + 1) \\\ \end{align}
$$
$f(g) = \sum_{g = \gcd(P(i, j)), i < j} (d(i, j) + 1)$として、これの倍数ゼータ変換$F(n)$を考えましょう。
$$
\begin{align}
F(n) = & \sum_{n | m} \sum_{m = \gcd(P(i, j)), i < j} (d(i, j) + 1) \\\ = & \sum_{n | \gcd(P(i, j)), i < j} (d(i, j) + 1) \\\ = & \sum_{P(i, j)の全ての要素はnを約数にもつ, i < j} (d(i, j) + 1)
\end{align}
$$
なので、$A_i$が$n$を約数にもつような頂点$i$だけからなるグラフを考えると計算できそうです。パスの長さの総和は、各辺を合計何回通るかで計算できます。
dfs戻りがけを意識して実装すると簡単です。
あとは同様に$F(n)$を倍数メビウス変換をすれば、$f(n)$が求まります。
提出コード
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
| i64 N;
vector<i64> A;
vector<vector<int>> G;
vector<vector<int>> dv(1e5 + 1);
vector<int> vis;
void dfs(int v, int p, int g, vector<i64>& cnt) {
vis[v] = 1;
i64 now = 1;
for(auto u: G[v]) {
if(u == p) continue;
if(A[u] % g != 0) continue;
dfs(u, v, g, cnt);
now += cnt.back();
}
cnt.push_back(now);
}
int main() {
cin >> N;
A.resize(N);
G.resize(N);
rep(i,0,N) cin >> A[i];
rep(i,0,N - 1) {
i64 u, v;
cin >> u >> v;
u--;
v--;
G[u].push_back(v);
G[v].push_back(u);
}
rep(i,0,N) {
for(i64 d = 1; d * d <= A[i]; d++) {
if(A[i] % d == 0) {
dv[d].push_back(i);
if(A[i] / d != d) {
dv[A[i] / d].push_back(i);
}
}
}
}
vector<fp> ans(dv.size());
vis.resize(N, 0);
rep(g,1,dv.size()) {
for(auto v: dv[g]){
if(vis[v] == 0) {
vector<i64> cnt;
dfs(v, -1, g, cnt);
rep(i,0,cnt.size() - 1) {
ans[g] += fp(cnt[i]) * fp(cnt.back() - cnt[i]);
}
ans[g] += fp(cnt.back()) * fp(cnt.back() - 1) / fp(2);
}
}
for(auto v: dv[g]) {
vis[v] = 0;
}
}
inverse_multiple_transform(ans);
fp res;
rep(i,1,dv.size()) {
res += fp(i) * ans[i];
}
cout << res << endl;
}
|