地址
https://codeforces.com/contest/1139/problem/D
原文地址
https://www.lucien.ink/archives/409
题意
模拟一个循环,一开始有一个空序列,之后每次循环:
- 从 $[1, m]$ 中等概率选出一个数字添加到序列里去。
- 检查这个序列所有元素的 $gcd$ 是否为 $1$,如果为 $1$ 则停止,若否则重复 $1$ 操作直至 $gcd$ 为 $1$ 为止。
求这个序列的长度期望。
题解
(废了很大的劲才理解朴素做法,并不会队友 $O(m \cdot log_2m)$ 的写法,Orz
好菜啊)。
设 $f_i$ 为当前序列的 $gcd$ 为 $i$ 时,将这个序列的 $gcd$ 变为 $1$ 的期望次数。
显然有 $f_1 = 0$ 。
接下来讨论 $i > 1$ 的情况,根据意义我们可以得到:
$$f_i = 1 + \sum_{j = 1}^{m} \frac{f_{gcd(i, j)}}{m}$$
这是一个优秀的 $O(m^2log_2{10^5})$ 的超时算法,考虑优化。
观察上式子中的 $gcd(i, j)$ ,我们发现可以只枚举 $i$ 的因子,设 $calc(x, y)$ 为计算集合 $\{i | i \in [1, m], gcd(x, i) = y\}$ 的元素数量 $|set|$ ,那么上式就演化为:
$$f_i = 1 + \sum_{k \in divisors(i)} \frac{calc(i, k)}{m} \cdot f_k$$
发现等式右边也会含有 $f_i$ ,为了消除影响,将其挪至等式左侧并化简:
$$f_i = 1 + \frac{\lfloor \frac{m}{i} \rfloor}{m} \cdot f_i + \sum_{k \in divisors(i) \cap k \neq i} \frac{calc(i, k)}{m} \cdot f_k$$
$$f_i - \frac{\lfloor \frac{m}{i} \rfloor}{m} \cdot f_i = 1 + \sum_{k \in divisors(i) \cap k \neq i} \frac{calc(i, k)}{m} \cdot f_k$$
$$f_i \cdot (1 - \frac{\lfloor \frac{m}{i} \rfloor}{m}) = 1 + \sum_{k \in divisors(i) \cap k \neq i} \frac{calc(i, k)}{m} \cdot f_k$$
$$f_i = \frac{1 + \sum_{k \in divisors(i) \cap k \neq i} \frac{calc(i, k)}{m} \cdot f_k}{1 - \frac{\lfloor \frac{m}{i} \rfloor}{m}}$$
忽略 $calc$ 的复杂度,那么外围复杂度为 $O(mlog_2m)$ ,下面考虑 $calc$ 函数:
对于集合 $set = \{i | i \in [1, m], gcd(x, i) = y\}$ ,令:
$$ \begin{cases} x = k_x \cdot y& (1 \leq k_x, k_x \in Z)\\ i = k_i \cdot y& (1 \leq k_i, k_i \in Z) \end{cases} $$
那么:
$$ \begin{cases} x = k_x \cdot y \\ i = k_i \cdot y \\ 1 \leq i \leq m \\ gcd(x, i) = y \end{cases} \Rightarrow \begin{cases} k_x = \frac{x}{y} \\ k_i = \frac{i}{y} \\ 1 \leq k_i \leq \frac{m}{y} \\ gcd(k_x, k_i) = 1 \end{cases} $$
显然,当且仅当 $m' = \frac{m}{y}, x' = \frac{x}{y}$ 时,集合 $set' = \{i | i \in [1, m'], gcd(x', i) = 1\}$ 的元素数量 $|set'| = |set|$ 。
考虑计算 $|set'|$ ,那么问题就变成了,给定 $m'$ 和 $x'$ ,求 $[1, m']$ 中与 $x'$ 互质的元素的数量,设 $n$ 为 $x'$ 的质因子的数量,这个子问题的复杂度为 $O(2 ^ n \cdot n)$。
最后的总复杂度为 $O(m \cdot log_2m \cdot 2 ^ n \cdot n)$ 。
代码
#include <bits/stdc++.h>
typedef long long ll;
const int maxn = int(1e5) + 7, mod = int(1e9) + 7;
bool not_prime[maxn];
std::vector<int> factor[maxn], prime_factor[maxn];
int f[maxn], m;
ll qpow(ll p, ll q) {
ll ret = 1;
while (q) {
if (q & 1) ret = ret * p % mod;
p = p * p % mod;
q >>= 1;
}
return ret;
}
void init() {
for (int i = 2; i < maxn; i++) {
for (int j = 2; i * j < maxn; j++) factor[i * j].push_back(i);
if (!not_prime[i]) {
for (int j = 1; i * j < maxn; j++) {
int num = i * j;
not_prime[num] = true;
prime_factor[num].push_back(i);
}
}
}
}
int __calc(int a, int b) {
int ret = b;
std::vector<int> buf;
for (int &it : prime_factor[a]) buf.push_back(it);
int len = int(buf.size());
for (int i = 1; i < 1 << len; i++) {
int tmp = 1;
for (int j = 0; j < len; j++) {
if ((i >> j) & 1) {
tmp = -tmp * buf[j];
}
}
ret += b / tmp;
}
return ret;
}
int calc(int x, int y) {
return __calc(x / y, m / y);
}
int main() {
init();
std::cin >> m;
int inv_m = int(qpow(m, mod - 2)), ans = 1;
for (int i = 2; i <= m; i++) {
f[i] = 1;
for (auto &k : factor[i]) {
f[i] = int((f[i] + 1ll * calc(i, k) * f[k] % mod * inv_m % mod) % mod);
}
f[i] = int(1ll * f[i] * qpow((1 - 1ll * (m / i) * inv_m % mod + mod) % mod, mod - 2) % mod);
ans = int((ans + 1ll * f[i] * inv_m % mod) % mod);
}
std::cout << ans << std::endl;
return 0;
}
2 条评论
所以 __calc() 计算互质个数是用容斥吧
是的