题意
初始有一个空数组\(a\)，接下来每次操作会这么做：
- 在\([1,n]\)中选择一个数，将其拼接在数组\(a\)后。
计算数组\(a\)的\(\gcd\)。
如果结果是\(1\)，退出。
否则，回到步骤1.
试问数组\(a\)长度的期望，答案对\(1e^9+7\)取模。
题解
莫反解法
考虑dp。
设\(dp[x]\)为当前数组的\(\gcd\)为\(x\)时的期望长度，那么答案便是
\[ans=\sum_{i=1}^n\frac{1+dp[i]}{n}=1+\sum_{i=1}^n\frac{dp[i]}{n}
\]
由\(dp[x]\)的意义，考虑对其因子进行转移，即
\[dp[x]=1+\frac{1}{n}\sum_{i=1}^ndp[\gcd(i,x)]
\]
然而\(O(n^2)\)的时间复杂度必定是超时的，于是考虑优化。
对于\(\gcd(i,x)\)，它必定是\(x\)的因子，于是可以优化成
\[dp[x]=1+\frac{1}{n}\sum_{d|x}dp[d]*h(x,d)
\]
其中\(h(x,d)\)表示\([1,n]\)中与\(x\)的\(\gcd\)为\(d\)的个数，即
\[h(x)=\sum_{i=1}^n[\gcd(x,i)==d]
\]
对于优化后的\(dp[x]\)的表达式，由于左右两边可能都含有\(dp[x]\)，于是我们把\(dp[x]\)单独分离出来
\[n*dp[x]=n+\sum_{d|x,d\neq x}dp[d]*h(x,d)+dp[x]*h(x,x)
\]
其中
\[h(x,x)=\lfloor \frac{n}{x} \rfloor
\]
通过移项，变形得
\[(n-\lfloor\frac{n}{x}\rfloor)dp[x]=n+\sum_{d|x,d\neq x}dp[d]*h(x,d)
\]
接下来再分析\(h(x,d)\)：
\[h(x,d)=\sum_{i=1}^n [\gcd(x,i)==d]=\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}[\gcd(\frac{x}{d}, i)==1]=\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{t|\frac{x}{d}, t|i}\mu(t)
\]
变换求和顺序得
\[h(x,d)=\sum_{t|\frac{x}{d}}\sum_{i=1}^{{\lfloor \frac{n}{d}\rfloor}}[t|i]\mu(t)=\sum_{t|\frac{x}{d}}\lfloor\frac{n}{dt}\rfloor \mu(t)
\]
带入原式
\[dp[x]=\frac{n+\sum_{d|x,d\neq x}(dp[d]*\sum_{t|\frac{x}{d}}\lfloor\frac{n}{dt}\rfloor \mu(t))}{n-\lfloor\frac{n}{x}\rfloor}
\]
AC代码
#include <bits/stdc++.h>
#define IO ios::sync_with_stdio(NULL)
#define sc(z) scanf("%d", &(z))
#define _ff(i, a, b) for (ll i = a; i <= b; ++i)
#define _rr(i, a, b) for (ll i = b; i >= a; --i)
#define _f(i, a, b) for (ll i = a; i < b; ++i)
#define _r(i, a, b) for (ll i = b - 1; i >= a; --i)
#define pii pair<int, int>
#define mkp make_pair
#define endl "\n"
using namespace std;
typedef long long ll;
const int N = 1e6 + 5;
const ll mod = 1e9 + 7;
ll dp[N];
vector<ll> mind[N];
bool prime[N];
ll primes[N], cnt, mu[N], inv[N];
void precompute() {
mu[1] = 1, cnt = 0;
_ff(i, 2, N) {
if (!prime[i]) {
primes[cnt++] = i;
mu[i] = -1;
}
_f(j, 0, cnt) {
if (primes[j] * i > N) break;
prime[i * primes[j]] = 1;
if (i % primes[j]) {
mu[i * primes[j]] = -mu[i];
} else {
mu[i * primes[j]] = 0;
break;
}
}
}
_ff(i, 1, N) {
mu[i] = (mu[i] + mod) % mod;
for (int j = i + i; j < N; j += i) mind[j].push_back(i);
}
inv[1] = 1;
_ff(i, 2, N) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
}
int main() {
ll n; cin >> n;
precompute();
_ff(x, 2, n) {
dp[x] = n;
for (ll d : mind[x]) {
//由于之前筛因子的时候并没有把自身筛进去，于是这里要单独计算t=x/d的情况
ll s = n / x * mu[x / d] % mod;
for (ll t : mind[x / d]) {
s = (s + n / d / t * mu[t]) % mod;
}
dp[x] = (dp[x] + s * dp[d]) % mod;
}
dp[x] = dp[x] * inv[n - n / x] % mod;
}
ll ans = 0;
_ff(i, 1, n) ans = (ans + dp[i]) % mod;
ans = ans * inv[n] % mod;
cout << ans + 1 % mod;
return 0;
}
无穷级数解法
根据期望的定义，我们可以得到
\[ans=\sum_{i=1}^niP(i)
\]
如何计算\(P(i)\)呢？我们设\(P(E)\)为事件\(E\)发生的概率，设当前长度为\(len\),那么我们可以计算\(len\)的期望\(P(len>x)\):
\[E(len)=1+\sum_{x=1}^{\infty}P(len>x)
\]
注意到\(len>x\)就意味着所有前\(x\)个数是不互质的，那么\(P(len>x)\)就可以表示为
\[P(len>x)=\sum_{i=2}^nP(\gcd(a_1,a_2,...,a_x)=i)
\]
我们令\(ans_i=\sum_{x=1}^{\infty}P(\gcd(a_1,a_2,...,a_x)=i)\)，那么最终答案就为
\[ans=1+\sum_{i=2}^nans_i
\]
那么如何计算\(ans_i\)呢？我们可以根据容斥原理，先计算前\(x\)个数都是x的倍数的概率，再减去前\(x\)个数的\(\gcd\)是\(i\)的倍数的概率，即
\[ans_i=\sum_{x=2}^{\infty}P(gcd(a_1,a_2,...,a_x)=i)=\sum_{x=2}^{\infty}P(前x个数都能整除i)-\sum_{j=2}^{\lfloor\frac{n}{i}\rfloor}ans_{j*i}
\]
我们很容易知道，在\([1,n]\)中取一个\(i\)的倍数的概率为\(\frac{\lfloor \frac{n}{i}\rfloor}{n}\)，那么原式就为
\[ans_i=\sum_{x=1}^{\infty}(\frac{\lfloor \frac{n}{i}\rfloor}{n})^x-\sum_{j=2}^{\lfloor\frac{n}{i}\rfloor}ans_{j*i}
\]
很容易知道，级数\(\sum_{x=1}^{\infty}(\frac{\lfloor \frac{n}{i}\rfloor}{n})^x\)是收敛的，收敛为\(\frac{1}{1-\frac{\lfloor \frac{n}{i}\rfloor}{n}}-1\)。
AC代码
#include <bits/stdc++.h>
#define IO ios::sync_with_stdio(NULL)
#define sc(z) scanf("%d", &(z))
#define _ff(i, a, b) for (int i = a; i <= b; ++i)
#define _rr(i, a, b) for (int i = b; i >= a; --i)
#define _f(i, a, b) for (int i = a; i < b; ++i)
#define _r(i, a, b) for (int i = b - 1; i >= a; --i)
#define pii pair<int, int>
#define mp make_pair
using namespace std;
typedef long long ll;
const int N = 1e5 + 5;
const ll mod = 1e9 + 7;
ll ans[N];
ll qpow(ll a, ll b) {
ll res = 1;
while (b) {
if (b & 1) res = (res * a) % mod;
a = (a * a) % mod;
b >>= 1;
}
return res;
}
ll inv(ll x) { return qpow(x, mod - 2); }
int main() {
ll n; cin >> n;
ll dans = 1;
_rr(i, 2, n) {
ans[i] = (inv((1 - n / i * inv(n) + mod) % mod) - 1 + mod) % mod;
for (ll j = i + i; j <= n; j += i) {
ans[i] = (ans[i] - ans[j] + mod) % mod;
}
dans = (dans + ans[i]) % mod;
}
cout << dans;
return 0;
}
