[2021 CCCC天梯赛] 可怜的简单题 (概率期望 莫比乌斯反演 杜教筛)

题意

每次从 $[1,n]$ 中选择一个数加到一个序列末尾，当 $\gcd (a_1,\cdots ,a_n)=1$ 时停止，求期望长度，对 $p$ 取模

$1\le n\le 10^{11},n<p\le 10^{12}$

分析：

设 $E(x)$ 为长度为 $x$ 的期望，那么根据期望定义

$$E(x)=\sum _{i=1}^{\infty }P(x=i)\times i$$

把 $i$ 改为 $\sum _{j=1}^i$

$$E(x)=\sum _{i=1}^{\infty }P(x=i)\sum _{j=1}^i$$

交换求和次序

$$\sum _{i=1}^{\infty }\sum _{j=i}^{\infty }P(x=j)$$

等价于

$$\sum _{i=1}^{\infty }P(x\ge i)$$

化简一下

$$\sum _{i=1}^{\infty }P(x\ge i)=1+\sum _{i=1}^{\infty }P(x>i)$$

考虑 $P(x>i)$，进行容斥 $1-P(x\le i)$ 就等价于

$$1-P(\gcd (a_1,\cdots ,a_i)=1)$$

枚举 $a_i$ 在 $[1,n]$ 中的取值

$$1-\sum _{a_1=1}^n\cdots \sum _{a_i=1}^n\frac{[\gcd (a_1,\cdots ,a_i)=1]}{n^i}$$

莫比乌斯反演

$$1-\sum _{a_1=1}^n\cdots \sum _{a_i=1}^n\frac{\sum _{d\mid \gcd (a_1,\cdots ,a_i)}\mu (d)}{n^i}$$

交换求和次序

$$1-\frac{\sum _{d=1}^n\mu (d)\lfloor \frac{n}{d}{\rfloor }^i}{n^i}$$

把 $1$ 拿到分子，和第一项抵消了

$$-\frac{\sum _{d=2}^n\mu (d)\lfloor \frac{n}{d}{\rfloor }^i}{n^i}$$

代入到 $1+\sum _{i=1}^{\infty }P(len>i)$ 得

$$1-\sum _{i=1}^{\infty }\frac{\sum _{d=2}^n\mu (d)\lfloor \frac{n}{d}{\rfloor }^i}{n^i}$$

交换求和次序

$$1-\sum _{d=2}^n\mu (d)\sum _{i=1}^{\infty }(\frac{\lfloor \frac{n}{d}\rfloor }{n}{)}^i$$

$\sum _{i=1}^{\infty }(\frac{\lfloor \frac{n}{d}\rfloor }{n}{)}^i$ 这是个等比级数，极限为 $\frac{首项}{1-公比}$

$$1-\sum _{d=2}^n\mu (d)\frac{\lfloor \frac{n}{d}\rfloor }{n-\lfloor \frac{n}{d}\rfloor }$$

就可以用杜教筛了

代码：

#include <bits/stdc++.h>
#define int long long
using namespace std;
signed main() {
    cin.tie(0) -> sync_with_stdio(0);
    int n, mod;
    cin >> n >> mod;
    int cnt = 0, N = 2.2e7 + 5;
    vector<int> primes(N), mobius(N), sum(N);
    vector<bool> st(N);
    auto sieve = [&](int n) {
        mobius[1] = 1;
        for (int i = 2; i <= n; i ++) {
            if (!st[i]) {
                primes[cnt ++] = i;
                mobius[i] = -1;
            }
            for (int j = 0; i * primes[j] <= n; j ++) {
                int t = i * primes[j];
                st[t] = 1;
                if (i % primes[j] == 0) {
                    mobius[t] = 0;
                    break;
                }
                mobius[t] = -mobius[i];
            }
        }
        for (int i = 1; i <= n; i ++) {
        	sum[i] = (sum[i - 1] + mobius[i] + mod) % mod;
        }
    };
    sieve(N - 1);
    auto qmul = [&](int a, int b) {
    	return (__int128)a * b % mod;
    };
    auto qmi = [&](int a, int b) {
    	int res = 1;
    	while (b) {
    		if (b & 1) res = qmul(res, a);
    		a = qmul(a, a);
    		b >>= 1;
    	}
    	return res;
    };
    unordered_map<int, int> mp;
    function<int(int)> Sum = [&](int n) {
    	if (n < N) return sum[n];
    	if (mp[n]) return mp[n];
    	int res = 1;
    	for (int l = 2, r; l <= n; l = r + 1) {
    		r = n / (n / l);
    		res = (res - qmul(r - l + 1, Sum(n / l)) % mod + mod) % mod;
    	}
    	return mp[n] = res;
    };
    int res = 1;
    for (int l = 2, r; l <= n; l = r + 1) {
    	r = n / (n / l);
    	int t = qmul(n / l, qmi(n - n / l, mod - 2));
    	res = (res - qmul(Sum(r) - Sum(l - 1), t) + mod) % mod;
    }
    cout << res << endl;
}