# 地址

https://codeforces.com/contest/1139/problem/D

# 原文地址

https://www.lucien.ink/archives/409

## 题意

1. 从 $[1, m]$ 中等概率选出一个数字添加到序列里去。
2. 检查这个序列所有元素的 $gcd$ 是否为 $1$，如果为 $1$ 则停止，若否则重复 $1$ 操作直至 $gcd$ 为 $1$ 为止。

## 题解

（废了很大的劲才理解朴素做法，并不会队友 $O(m \cdot log_2m)$ 的写法，Orz 好菜啊）。

$$f_i = 1 + \sum_{j = 1}^{m} \frac{f_{gcd(i, j)}}{m}$$

$$f_i = 1 + \sum_{k \in divisors(i)} \frac{calc(i, k)}{m} \cdot f_k$$

$$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}}$$

$$\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}$$

## 代码

https://pasteme.cn/4936

#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;
}