整除分块, 莫比乌斯函数, 莫比乌斯反演
整除分块
对于式子
\[\sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor\]其可以在 $O(n)$ 的时间复杂度内完成计算，但这个复杂度不够优秀，通过整除分块可以将复杂度降低至 $O(\sqrt n)$。
通过直接观察，我们可以发现，对于每一个值相同的分块，其最后一个数是 $\frac{n}{n / i}$。
整除分块的实现：
for(int l = 1; l <= n; l = r + 1){
r = n / (n / i);
ans += (r) * (n / l);
}
莫比乌斯函数
\[\mu(n)=\begin{cases}1, &n = 1 \\0 ,& n\text{含有平方因子} \\(-1)^k,& k\text{为}n\text{的本质不同质因子个数}\end{cases}\]通俗的来讲，莫比乌斯函数可以被如下的语句描述：
- 莫比乌斯函数 $\mu(n)$ 的定义域是 $N$；
- $\mu(1)=1$；
- 当 $n$ 存在平方因子时，$\mu(n)=0$；
- 当 $n$ 是素数或奇数个不同素数之积时，$\mu(n)=-1$；
- 当 $n$ 是偶数个不同素数之积时，$\mu(n)=1$。
莫比乌斯函数的部分重要性质：
性质一：
\[\sum_{d \mid n} \mu(d) = \begin{cases} 1, &n=1 \\ 0, &n\neq1 \end{cases}\]证明：
若 $n=1$，显然成立。
若 $n \neq 1$，令 $n=\prod_{i=1}^{k} p_i^{c_i}, n^{\prime}=\prod_{i=1}^{k}p_i$，有
$\sum_{d \mid n} \mu(d) = \sum_{d \mid n^{\prime}} \mu(d)=\sum_{i=0}^k C_k^i \cdot (-1)^i=(1 + (-1))^k=0$。
性质二：
\[\sum_{d \mid n} \frac{\mu(d)}{d} = \frac{\varphi(n)}{n}\]使用迪利克雷卷积后两边同时除以 $n$ 即可得到。
性质三：
莫比乌斯函数是不完全积性函数：
\[\mu(nm) = \mu(m) \cdot \mu(n), \text{当且仅当}gcd(n,m) = 1\]给出线性筛法求莫比乌斯函数的方法：
void get_mu(int n){
mu[1] = 1;
for(int i = 2; i <= n; i++){
if(!st[i]){prime[++cnt] = i; mu[i] = -1;}
for(int j = 1; j <= cnt && prime[j] * i <= n; j++)
{
st[prime[j] * i] = 1;
if(i % prim[j] == 0) break;
else mu[i * prime[j]] = -mu[i];
}
}
}
迪利克雷卷积的三个重要公式：
\[\mu \times I = \sum_{d \mid n} \mu(d) I(\frac{n}{d}) = \sum_{d \mid n} \mu(d) = [n = 1] = \varepsilon.\]如果 $n \neq 1$，有
\[\sum_{d \mid n} = \mu(d) = \sum_{i = 0}^{k} C_{k}^{i} (-1)^{i} = 0.\]并且有
\[\varphi \times I = \sum_{d \mid n} \varphi(d) I(\frac{n}{d}) = \sum_{d \mid n} \varphi(d) = \text{id}(n).\]上述公式中 $\sum_{d \mid n} \varphi(d) = \text{id}(n)$ 的证明： 对集合：${\frac{1}{n}, \frac{2}{n}, \dots, \frac{n}{n}}$， 可以化为最简形式，对于最简形式 $\frac{a}{b}$，有 $a \leq b, \gcd(a, b) = 1, b \mid n$。 那么，以 $b$ 为分母的集合中的分数有 $\varphi(b)$ 个，一共有 $n$ 个元素。 由此，$\sum_{d \mid n} \varphi(d) = n$。
综上，$\mu \times \text{id} = \varphi$，$\varphi \times I = \text{id}$，$\varphi \times I \times \mu = \text{id} \times \mu = \varphi$。
莫比乌斯反演
定义：对于定义在正整数集合上的两个函数 $g(n)$，$f(n)$，若他们满足条件：
\[g(n)=\sum_{d \mid n} f(d)\]那么有结论：
\[f(n) = \sum_{d \mid n} \mu(d)g(\frac{n}{d})\]证明：
由
\[g(n)=\sum_{d \mid n} f(d)\]有
\[g(\frac{n}{d})=\sum_{k \mid \frac{n}{d}}f(k)\]那么
\[\sum_{d \mid n} \mu(d)g(\frac{n}{d}) = \sum_{d \mid n} \mu(d) \sum_{k \mid \frac{n}{d}}f(k)\]变换求和次序：
\[\sum_{d \mid n} \mu(d)g(\frac{n}{d}) = \sum_{d \mid n} \sum_{k \mid \frac{n}{d}} \mu(d) f(k)=\sum_{k \mid n} \sum_{d \mid \frac{n}{k}} \mu(d) f(k) =\sum_{k \mid n}f(k)\sum_{d \mid \frac{n}{k}} \mu(d)\]由性质一，当 $\frac{n}{k} = 1$ 即 $n = k$ 时才有 $\sum_{d \mid \frac{n}{k}} \mu(d) = 1$，此时即得 $\sum_{d \mid n} \mu(d)g(\frac{n}{d}) = f(n)$。
莫比乌斯反演的倍数形式：
\[g(n)=\sum_{n \mid d}f(d) \iff f(n) = \sum_{n \mid d} g(d) \mu(\frac{d}{n})\]反演练习题
Crash的数字表格
题目描述：对 $n, m(1\leq n, m \leq 10^7)$，求 $\sum_{i = 1}^{n} \sum_{j = 1}^m lcm(i,j)$。
思路：
先将求最小公倍数变为求最大公约数，再改写枚举最大公因数的形式：
\[\begin{aligned} \sum_{i = 1}^{n} \sum_{j = 1}^m lcm(i,j) &= \sum_{i = 1}^{n} \sum_{j = 1}^m \frac{ij}{gcd(i,j)} \\ &= \sum_{i = 1}^{n} \sum_{j = 1}^m \sum_{d \mid i,d \mid j, gcd(\frac{i}{d}, \frac{j}{d})=1}\frac{ij}{d} \end{aligned}\]改变求和次序：
\[=\sum_{d=1}^nd \sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[gcd(i,j)=1]\cdot ij\]可以发现，$\sum_{d=1}^nd$ 很好求，现在需要对后面的式子进行一定的化简。
令
\[sum(n,m) = \sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=1]\cdot ij\]枚举约数：
\[sum(n,m) =\sum_{d=1}^n \sum_{d \mid i}^n\sum_{d\mid j}^m \mu(d) \cdot ij\]令 $i = d i^{\prime}, j=dj^{\prime}$，有：
\[\begin{aligned} sum(n,m) &=\sum_{d=1}^n \mu(d)\cdot d^2 \sum_{i^{\prime}=1}^{\frac{n}{d}}\sum_{j^{\prime}=1}^{\frac{m}{d}}i^{\prime}j^{\prime} \\ &=\sum_{d=1}^n \mu(d)\cdot d^2 \cdot \frac{(1 + \frac{n}{d})(\frac{n}{d})}{2}\cdot \frac{(1 + \frac{m}{d})(\frac{m}{d})}{2} \end{aligned}\]对于 $\sum_{d=1}^{n}\mu(d) \cdot d^2$，可以对其进行 $O(n)$ 的预处理前缀和，后面的部分可以 $O(1)$ 计算。
将 $sum(n,m)$ 代回原式，得到：
\[ans = \sum_{d=1}^nd \cdot sum(\lfloor \frac{n}{d} \rfloor, \lfloor \frac{m}{d} \rfloor)\]#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e7 + 1;
int n, m;
int mod = 20101009;
int primes[N], cnt, st[N], mu[N], smu[N], sum[N];
void init(){ // 线性求莫比乌斯函数
mu[1] = 1;
int k = min(m, n);
for(int i = 2; i <= k; ++i){
if(!st[i]) primes[cnt++] = i, mu[i] = -1;
for(int j = 0; primes[j] * i < N; ++j){
st[primes[j] * i] = 1;
if(i % primes[j] == 0) {
mu[i * primes[j]] = 0;
break;
}
mu[primes[j] * i] = -mu[i];
}
}
for(int i = 1; i <= k; ++i){
smu[i] = smu[i - 1] + mu[i];
}
for(int i = 1; i <= k; ++i){ // O(k) 预处理 mu[d] * d^2
sum[i] = (sum[i - 1] + (LL) i * i % mod * (mu[i] + mod) % mod) % mod;
}
}
int g(int x, int i){ // 整除分块函数
return x / (x / i);
}
LL calc(int a, int b){ // O(1) 求 sum 的后半部分
return (LL)(1 + a) * a / 2 % mod * ((LL)(1 + b) * b / 2 % mod) % mod;
}
int func(int a, int b){ // 即公式中的 sum(n, m) 函数
int res = 0, r;
for(int l = 1; l <= min(a, b); l = r + 1){
r = min(g(a, l), g(b, l));
res = (res + (LL)(sum[r] - sum[l - 1]) * calc(a / l, b / l) % mod) % mod;
}
return res;
}
int work(){
int res = 0;
for(int l = 1, r; l <= min(n, m); l = r + 1){
r = min(g(n, l), g(m, l));
res = (res + (LL)(r - l + 1) * (l + r) / 2 % mod * func(n / l, m / l) % mod) % mod;
}
return (res + mod) % mod;
}
int main(){
cin >> n >> m;
init();
cout << work() << endl;
return 0;
}
简单题
题目描述：对 $n, k(1\leq n \leq 5*10^6, 1 \leq k \leq 10^{18})$，求：
\[\left(\sum_{i=1}^{n}\sum_{j=1}^n(i+j)^k\cdot \mu^2(\gcd(i,j))\cdot \gcd(i,j)\right)\bmod 998244353\]思路：
枚举 $\gcd$，有：
\[\sum_{d=1}^{n}d \mu^2(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} (di+dj)^k[\gcd(i,j)=1]\]代换 $[\gcd(i,j)=1]$ 并交换求和次序：
\[\sum_{d=1}^{n}d \mu^2(d) \sum_{x=1}^{\lfloor \frac{n}{d} \rfloor} \mu(x) \sum_{i=1}^{\lfloor \frac{n}{dx} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{dx} \rfloor} (dxi+dxj)^k\]提取因子：
\[\sum_{d=1}^{n} d^{k+1} \mu^2(d) \sum_{x=1}^{\lfloor \frac{n}{d} \rfloor} \mu(x)\cdot x^k \sum_{i=1}^{\lfloor \frac{n}{dx} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{dx} \rfloor} (i+j)^k\]令 $dx=T$，有：
\[\sum_{T=1}^{n} T^k \sum_{i=1}^{\lfloor \frac{n}{T} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{T} \rfloor} (i+j)^k \sum_{d \mid T} \mu^2(d) \cdot d \cdot \mu(\lfloor \frac{T}{d} \rfloor)\]其中，$\sum_{T=1}^nT^k$ 是好求的，下面考虑后面的部分。
令
\[\begin{aligned} f(n)&= \sum_{d \mid n} \mu^2(d) \cdot d \cdot \mu(\lfloor \frac{n}{d} \rfloor) =(\mu^2 \cdot id) \ast\mu \\ p(n)&= \sum_{i=1}^{n} \sum_{j=1}^{n} (i+j)^k \end{aligned}\]首先考虑 $f(n)$，它显然是一个积性函数。
对 $n=p_1^{\alpha_1} \cdot p_2^{\alpha_2} \cdot \dots \cdot p_c^{\alpha_c}$，考虑 $c=1$，即 $n=p_1^{\alpha_1}$。
- 若 $\alpha_1=1$，$f(n)=p_1 - 1$。
- 若 $\alpha_1=2$，$f(n)=\mu^2(1) \cdot 1 \cdot \mu(p_1^2) + \mu^2(p_1) \cdot p_1 \cdot \mu(p_1) + \mu^2(p_1^2) \cdot p_1^2 \cdot \mu(1) = -p_1$。
- 若 $\alpha_1>2$，观察即发现 $f(n)$ 中必定有一个因子为 $0$，因此 $f(n)=0$。
据此，可以得到：
\[f(n)=\prod_{i=1}^c ([\alpha_i = 1](p_i - 1) + [\alpha_2=2](-p_i))\]此时，$f(n)$ 可以通过线性筛得到。
考虑 $p(n)$，通过画图可以发现其求和面积可等价为：
\[p(n)= \sum_{i=1}^{2n} \sum_{j=1}^{i} j^k - 2\sum_{i=1}^n \sum_{j=1}^{i} j^k\]此时，$p(n)$ 也可以通过自然数的 $k$ 次幂配合两次前缀和得到，需要注意的是，前缀和需要计算到 $2n$。
需要注意的是，在代码实现中将 $\sum_{T=1}^nT^k$ 和 $f(n)$ 放在一起计算了前缀和。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5 * 1e6 + 10, mod = 998244353;
ll n, k;
int f[N * 2 + 10], primes[N * 2 + 10], cnt, st[N * 2 + 10], tpowk[N * 2 + 10];
int qmi(ll a, ll b){
int res = 1;
while(b){
if(b & 1) res = (ll) res * a % mod;
a = (ll) a * a % mod;
b >>= 1;
}
return res % mod;
}
void init(){
f[1] = 1, tpowk[1] = 1;
for(int i = 2; i <= N << 1; ++i){
if(!st[i]) primes[cnt++] = i, f[i] = i - 1, tpowk[i] = qmi(i, k);
for(int j = 0; primes[j] * i <= N << 1; ++j){
st[primes[j] * i] = 1;
tpowk[primes[j] * i] = (ll) tpowk[primes[j]] * tpowk[i] % mod;
if(i % primes[j] == 0){
int q = i / primes[j];
if(q % primes[j]){
f[primes[j] * i] = (ll) (mod - primes[j]) * f[q] % mod;
}
break;
}
f[primes[j] * i] = (ll) f[i] * (primes[j] - 1) % mod;
}
}
for(int i = 1; i <= N << 1; ++i){
f[i] = ((((ll) f[i - 1] + (ll) f[i] * tpowk[i] % mod) + mod) % mod)% mod;
tpowk[i] = ((((ll) tpowk[i - 1] + tpowk[i]) % mod + mod) % mod) % mod;
}
for(int i = 1; i <= N << 1; ++i){
tpowk[i] = (((ll) tpowk[i - 1] + tpowk[i]) % mod + mod) % mod;
}
}
int main(){
cin >> n >> k;
init();
ll res = 0;
for(int l = 1, r; l <= n; l = r + 1){
r = n / (n / l);
res = (res + (ll) ((((ll) f[r] - f[l - 1]) % mod + mod) % mod) * ((((ll) tpowk[2 * (n / l)] - 2 * tpowk[(n / l)]) % mod + mod) % mod) % mod) % mod;
}
cout << res << endl;
return 0;
}
亚线性筛
杜教筛
杜教筛是一种在亚线性时间复杂度下求出积性函数前缀和的一种方法，其主要思想如下。
$f$ 为积性函数，要求 $S(n) = \sum_{i = 1}^{n} f(i)$。 我们可以尝试构造一个积性函数 $g$，考虑 $f \times g$ 的前缀和：$\sum_{i = 1}^{n} (f \times g)(i)$。 对上式进行变形：
\[\begin{aligned} \sum_{i = 1}^{n} (f \times g)(i) &= \sum_{i = 1}^{n} \sum_{d \mid i} f(d) g(\frac{i}{d}) \\ &= \sum_{d = 1}^{n} g(d) \sum_{i = 1}^{\frac{n}{d}} f(i) \\ &= \sum_{d = 1}^{n} g(d) S(\frac{n}{d}). \end{aligned}\]考虑 $g(1)S(n)$：
\[g(1)S(n) = \sum_{i = 1}^{n} g(i) S(\frac{n}{i}) - \sum_{i = 2}^{n} g(i) S(\frac{n}{i})\]又 $g(1) = 1$，有 $S(n) = \sum_{i = 1}^{n} g(i) S(\frac{n}{i}) - \sum_{i = 2} g(i) S(\frac{n}{i})$。 据此，时间复杂度为 $O(n^{\frac{3}{4}})$，可以线性筛预处理出前 $n^{\frac{2}{3}}$ 个答案优化复杂度到 $O(n^{\frac{2}{3}})$。
杜教筛模板题
题目描述：求 $\sum_{i=1}^n \varphi(i), \sum_{i=1}^n \mu(i), (1 \leq n \lt 2^{31})$。
思路：利用 $\mu \ast I = \epsilon, \varphi \ast I = id$ 进行杜教筛即可。
注意数论分块时 $l, r$ 存在越界的问题，需要使用 $long\ long$。此外，要通过记忆化加快筛法的速度。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1000000 + 10;
int primes[N], cnt;
bool st[N];
int mu[N];
ll phi[N];
unordered_map<int, int> sum_mu;
unordered_map<int, ll> sum_phi;
void init() {
mu[1] = 1, phi[1] = 1;
for(int i = 2; i < N; ++i) {
if(!st[i]) {
primes[cnt++] = i;
mu[i] = -1;
phi[i] = i - 1;
}
for(int j = 0; primes[j] * i < N; ++j) {
st[primes[j] * i] = 1;
if(i % primes[j] == 0) {
phi[i * primes[j]] = phi[i] * primes[j];
break;
}
phi[i * primes[j]] = phi[i] * (primes[j] - 1);
mu[i * primes[j]] = -mu[i];
}
}
for(int i = 1; i < N; ++i) {
mu[i] += mu[i - 1];
phi[i] += phi[i - 1];
}
}
ll get_phi(int x) {
if(x < N) return phi[x];
if(sum_phi[x]) return sum_phi[x];
ll res = (1ll + x) * x / 2;;
for(ll l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
res -= (r - l + 1) * get_phi(x / l);
}
return sum_phi[x] = res;
}
int get_mu(int x) {
if(x < N) return mu[x];
if(sum_mu[x]) return sum_mu[x];
int res = 1;
for(ll l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
res -= (r - l + 1) * get_mu(x / l);
}
return sum_mu[x] = res;
}
void work() {
int n;
cin >> n;
cout << get_phi(n) << " " << get_mu(n) << endl;
}
int main() {
cin.tie(0), cout.tie(0);
ios::sync_with_stdio(false);
int T;
init();
cin >> T;
while(T--) work();
return 0;
}
PN 筛
2022_杭电多校训练 round5 1002
思路：推公式题，之前没有见过 PN筛，通过这道题了解一下该算法并加深对迪利克雷卷积的理解。
公式推导过程如下。 有
\[\left| Good_i \right| = \prod C_k, (i = p_1^{C_1} \cdot p_2^{C_2} \cdot \dots \cdot p_k^{C_k})\]那么，答案可表示为 $ans = \frac{1}{m} \sum_{i = 1}^{m} \frac{i}{\prod C_k}$。 令 $f(x) = \frac{x}{\prod C_k}$，$g(x) = x$， 有 $f(x) = \sum_{d \mid x} h(d) g(\frac{x}{d})$。 那么，有
\[\begin{aligned} ans &= \frac{1}{m} \sum_{n = 1}^{m} f(n) \\ &= \frac{1}{m} \sum_{n = 1}^{m} \sum_{d \mid n} h(d) g(\frac{n}{d}) \\ &= \frac{1}{m} \sum_{d = 1}^{m} h(d) \sum_{n = 1}^{m} g(\frac{n}{d}) [d \mid n]. \end{aligned}\]上式中，
\[\sum_{n = 1}^{m} g(\frac{n}{d}) [d \mid n] = \sum_{n = 1}^{\frac{m}{d}} g(n) = \frac{(1 + \frac{m}{d}) \cdot \frac{m}{d}}{2}\]又 $f(p) = h(1)g(p) + h(p)g(1) = ph(1) + h(p) = p$，且 $h(p) = 0$， 利用结论 $\left| PN \right| \sim O(\sqrt{N})$，有
\[ans = \frac{1}{m} \sum_{d \in PN} h(d) \cdot \frac{(1 + \frac{m}{d}) \cdot \frac{m}{d}}{2}\]对 $\frac{p^k}{k}$ 有：
\[\begin{aligned} \frac{p^k}{k} &= f(n) \\ &= \sum_{d \mid n} h(d) g(\frac{n}{d}) \\ &= \sum_{i = 0}{k} h(p^i) g(p^{k - i}) \\ &= \sum_{i = 0}^{k} h(p^i) \cdot p^{k - i} \\ &= p \sum_{i = 0}^{k - 1} (h(p^i) \cdot p^{k - i - 1}) + h(p^k) \\ &= p f(p^{k - 1}) + h(p^k) \\ &= p \frac{p^{k - 1}}{k - 1} + h(p^k). \end{aligned}\]由上式可以得到 $h(p^k) = \frac{p^k}{k} - \frac{p^k}{k - 1} = -\frac{p^k}{k(k - 1)}$。
注意 $work$ 函数的搜索方式。
#pragma optimize(2)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N = 1e12 + 10, M = 1e7 + 10, mod = 4179340454199820289;
ll m;
ll primes[M], cnt;
ll inv[65];
bool st[M];
ll mul(ll a, ll b) {
return (a * b - (unsigned long long)((long double) a / mod * b) * mod + mod) % mod;
}
ll qmi(ll a, ll b) {
ll res = 1;
while(b) {
if(b & 1) res = mul(res, a);
a = mul(a, a);
b >>= 1;
}
return res;
}
void init() {
for(int k = 1; k <= 64; ++k) {
inv[k] = qmi(k, mod - 2);
}
ll mx = sqrt(N);
for(ll i = 2; i <= mx + 5; ++i) {
if(!st[i]) primes[cnt++] = i;
for(ll j = 0; i * primes[j] <= mx + 5; ++j) {
st[primes[j] * i] = 1;
if(i % primes[j] == 0) break;
}
}
}
ll calch(ll p, ll k, ll pk) {
return mod - mul(pk, mul(inv[k], inv[k - 1]));
}
ll work(ll n, ll idx, ll hd) {
ll res = mul(mul(n, mul(n + 1, inv[2])), hd);
for(ll i = idx; i < cnt; ++i) {
ll prime = primes[i], k = 1;
ll cur = n / prime, pk = prime;
if(cur < prime) break;
while(cur >= prime) {
cur /= prime;
pk *= prime;
++k;
res = (res + work(cur, i + 1, mul(hd, calch(prime, k, pk)))) % mod;
res = (res + mod) % mod;
}
}
return res;
}
int main() {
cin.tie(0), cout.tie(0);
ios::sync_with_stdio(false);
int T;
scanf("%d", &T);
init();
while(T--) {
scanf("%lld", &m);
printf("%lld\n", mul(work(m, 0, 1), qmi(m, mod - 2)));
}
return 0;
}
反演与其他
反演与期望
CF1139D Steps to One
题目描述
一个数列，初始时为空，每次随机选一个 $1$ 到 $m$ 之间的数加入数列末尾，数列中所有数的 $\gcd = 1$ 时停止操作。求数列的期望长度。
思路
设 $E(len)$ 为长度为 $len$ 的期望，那么显然有 $E(len) = \sum_{i = 1}^{\infty} p(len = i) * i$，其中，$p(len = x)$ 为数列长度为 $x$ 的概率。
对 $E(len)$ 进行处理：
\[\begin{aligned} E(len) &= \sum_{i= 1}^{\infty} p(len = i) * i \\ &= \sum_{i = 1}^{\infty}p(len = i) \sum_{j = 1}^{i}1 \\ &=\sum_{j = 1}^{\infty} \sum_{i = j}^{\infty}p(len = i) \\ &=\sum_{j = 1}^{\infty} p(len \geq j) \\ &= 1 + \sum_{i = 1}^{\infty} p(len \gt i) \end{aligned}\]考虑 $p(len>i)$：
\[\begin{aligned} p(len \gt i) &= p((\gcd_{j = 1}^{i} a_{j}) > 1) \\ &= 1 - p((\gcd_{j = 1}^{i}a_{j}) = 1) \\ &= 1 - \frac{\sum_{a_{1} = 1}^{m} \sum_{a_{2} = 1}^{m} \dots \sum_{a_{i} = 1}^{m} [(\gcd_{j = 1}^{i}a_{j}) = 1]} {m^{i}} \\ &= 1 - \frac{\sum_{a_{1} = 1}^{m} \sum_{a_{2} = 1}^{m} \dots \sum_{a_{i} = 1}^{m} \sum_{d \mid \gcd} \mu(d)} {m^{i}} \\ &= 1 - \frac{\sum_{d = 1}^{m} \mu(d) \cdot (\lfloor \frac{m}{d} \rfloor)^{i}} {m^{i}} \end{aligned}\]将 $p(len \gt i)$ 带回 $E(len)$：
\[\begin{aligned} E(len) &= 1 + \sum_{i = 1}^{\infty} (1 - \frac{\sum_{d = 1}^{m} \mu(d) (\lfloor \frac{m}{d} \rfloor)^{i}}{m^{i}}) \\ &= 1 - \sum_{i = 1}^{\infty} (\frac{\sum_{d = 2}^{m} \mu(d) (\lfloor \frac{m}{d} \rfloor)^{i}}{m^{i}}) \\ &= 1 - \sum_{i = 1}^{\infty} \frac{1}{m^{i}} \sum_{d = 2}^{m} \mu(d)(\lfloor \frac{m}{d} \rfloor)^{i} \\ &= 1 - \sum_{d=2}^{m}\mu(d) \sum_{i=1}^{\infty}(\frac{ \lfloor \frac{m}{d} \rfloor}{m})^{i} \\ &= 1 - \sum_{d=2}^{m}\mu(d) \frac{\lfloor \frac{m}{d} \rfloor}{m - \lfloor \frac{m}{d} \rfloor} \end{aligned}\]由此，可以预处理出莫比乌斯函数后进行递推，$O(m)$。
#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 10, mod = 1e9 + 7;
int m;
int primes[N], cnt;
bool st[N];
int mu[N];
int inv[N];
int qmi (int a, int b) {
int res = 1;
while (b) {
if (b & 1) {
res = (long long) res * a % mod;
}
a = (long long) a * a % mod;
b >>= 1;
}
return res;
}
void init () {
mu[1] = 1;
for (int i = 2; i <= m; ++i) {
if (!st[i]) {
primes[cnt++] = i;
mu[i] = -1;
}
for (int j = 0; primes[j] * i <= m; ++j) {
st[primes[j] * i] = 1;
if (i % primes[j] == 0) {
break;
}
mu[primes[j] * i] = -mu[i];
}
}
}
int main () {
cin >> m;
init();
int res = 0;
for (int i = 1; i <= m; ++i) {
inv[i] = qmi(i, mod - 2);
}
for (int d = 2; d <= m; ++d) {
res = ((res - (long long) mu[d] * (m / d) * inv[m - (m / d)]) % mod + mod) % mod;
}
res = (res + 1) % mod;
cout << res << endl;
return 0;
}