题目描述

可怜按照一种随机过程生成一个整数数列 $A$：
1. 初始时 $A$ 为空。
2. 从整数区间 $[1, n]$ 中等概率随机选择一个整数 $i$，将其加入数列 $A$ 的末尾。
3. 检查当前数列 $A$ 中的所有元素的最大公约数 (GCD)。如果 $\mathrm{GCD}(A) = 1$，则停止过程，返回 $A$。
4. 如果 $\mathrm{GCD}(A) > 1$，则返回步骤 2，继续添加元素。

任务是计算这个过程最终返回的数列 $A$ 的期望长度，结果对一个大质数 $p$ 取模。

样例

输入:

2 998244353

输出:

2

说明：
n=2。
第一次选 1，GCD=1，停止，长度 1。概率 1/2。
第一次选 2，GCD=2 > 1，继续。
第二次选 1，A={2, 1}，GCD=1，停止，长度 2。概率 (1/2)(1/2) = 1/4。
第二次选 2，A={2, 2}，GCD=2 > 1，继续。
第三次选 1，A={2, 2, 1}，GCD=1，停止，长度 3。概率 (1/2)^3 = 1/8。
第三次选 2，A={2, 2, 2}，GCD=2 > 1，继续…
期望长度 E = 1(1/2) + 2(1/4) + 3(1/8) + …
另一种计算方式：$E = 1 + P(\text{length}>1) + P(\text{length}>2) + \dots$
$P(\text{length}>k) = P(\mathrm{GCD}(A_1,\dots,A_k) > 1)$
$P(\text{length}>1) = P(GCD(A_1)>1) = P(A_1=2) = 1/2$.
$P(\text{length}>2) = P(GCD(A_1,A_2)>1)$. A1, A2 只能是 2。概率 (1/2)*(1/2) = 1/4.
$P(\text{length}>3) = P(GCD(A_1,A_2,A_3)>1)$. A1,A2,A3 只能是 2。概率 (1/2)^3 = 1/8.
$E = 1 + 1/2 + 1/4 + 1/8 + \dots = 1 + \frac{1/2}{1 - 1/2} = 1 + 1 = 2$.

输入:

100000000 998244353

输出:

3056898

算法 (期望 + 莫比乌斯反演 + 杜教筛)

$O(L + N^{2/3} + \sqrt{N} \log N)$

思路分析

1. 期望的线性性质:
   设 $L$ 是最终数列 $A$ 的长度。我们要求 $E[L]$。
   考虑随机过程的每一步。设 $A_k = \{i_1, i_2, \dots, i_k\}$ 是加入 $k$ 个元素后的数列。
   过程在第 $k$ 步停止的条件是 $\mathrm{GCD}(A_k) = 1$ 且 $\mathrm{GCD}(A_{k-1}) > 1$ (对于 $k>1$)，或者 $k=1$ 且 $\mathrm{GCD}(A_1)=1$。
   计算每一步停止的概率并求和 $E[L] = \sum_{k=1}^{\infty} k \cdot P(\text{stop at length } k)$ 比较复杂。

   我们使用期望的另一种计算方法。令 $X_k$ 为指示变量，当且仅当过程在长度达到 $k$ 后没有停止时，$X_k=1$，否则 $X_k=0$。即 $X_k=1$ 当且仅当 $\mathrm{GCD}(A_k) > 1$。
   最终的长度 $L$ 可以表示为 $L = 1 + \sum_{k=1}^{\infty} X_k$。（初始长度至少为 1，每当过程继续一次 ($X_k=1$)，长度就加 1）。
   根据期望的线性性质：
   $$E[L] = E[1] + \sum_{k=1}^{\infty} E[X_k] = 1 + \sum_{k=1}^{\infty} P(X_k = 1)$$
   $$E[L] = 1 + \sum_{k=1}^{\infty} P(\mathrm{GCD}(A_k) > 1)$$

2. 计算概率 $P(\mathrm{GCD}(A_k) > 1)$:
   直接计算 $\mathrm{GCD}(A_k) > 1$ 的概率比较困难，我们转而计算 $P(\mathrm{GCD}(A_k) = 1)$。
   $P(\mathrm{GCD}(A_k) > 1) = 1 - P(\mathrm{GCD}(A_k) = 1)$。
   $P(\mathrm{GCD}(A_k) = 1)$ 表示随机从 $[1, n]$ 中选择 $k$ 个数（可重复，独立同分布），它们的 GCD 为 1 的概率。
   利用莫比乌斯反演来计算这个概率。
   设 $f(d) = P(\mathrm{GCD}(A_k) = d)$。
   设 $g(d) = P(\mathrm{GCD}(A_k) \text{ 是 } d \text{ 的倍数})$。
   则 $g(d) = \sum_{m=1}^{\lfloor n/d \rfloor} f(m \cdot d)$。
   $g(d)$ 发生的条件是 $A_1, \dots, A_k$ 都是 $d$ 的倍数。
   在 $[1, n]$ 中，$d$ 的倍数有 $\lfloor n/d \rfloor$ 个。
   随机选一个数是 $d$ 的倍数的概率是 $\frac{\lfloor n/d \rfloor}{n}$。
   随机选 $k$ 个数都是 $d$ 的倍数的概率是 $\left( \frac{\lfloor n/d \rfloor}{n} \right)^k$。
   所以 $g(d) = \left( \frac{\lfloor n/d \rfloor}{n} \right)^k$。
   根据莫比乌斯反演公式：
   $$f(1) = \sum_{d=1}^{n} \mu(d) g(d)$$
   $$P(\mathrm{GCD}(A_k) = 1) = \sum_{d=1}^{n} \mu(d) \left( \frac{\lfloor n/d \rfloor}{n} \right)^k$$
   其中 $\mu(d)$ 是莫比乌斯函数。

3. 代入期望公式:
   $$E[L] = 1 + \sum_{k=1}^{\infty} \left( 1 - \sum_{d=1}^{n} \mu(d) \left( \frac{\lfloor n/d \rfloor}{n} \right)^k \right)$$
   注意到 $\sum_{d=1}^{n} \mu(d) (\dots)^k = \mu(1)(\lfloor n/1 \rfloor / n)^k + \sum_{d=2}^{n} \mu(d) (\dots)^k = 1 \cdot (n/n)^k + \sum_{d=2}^{n} \dots = 1 + \sum_{d=2}^{n} \mu(d) (\dots)^k$。
   所以 $1 - \sum_{d=1}^{n} \mu(d) (\dots)^k = -\sum_{d=2}^{n} \mu(d) \left( \frac{\lfloor n/d \rfloor}{n} \right)^k$。
   $$E[L] = 1 + \sum_{k=1}^{\infty} \left( -\sum_{d=2}^{n} \mu(d) \left( \frac{\lfloor n/d \rfloor}{n} \right)^k \right)$$
   交换求和顺序：
   $$E[L] = 1 - \sum_{d=2}^{n} \mu(d) \sum_{k=1}^{\infty} \left( \frac{\lfloor n/d \rfloor}{n} \right)^k$$

4. 计算几何级数:
   内部的和是一个无穷等比级数 $\sum_{k=1}^{\infty} r^k$，其中 $r = \frac{\lfloor n/d \rfloor}{n}$。
   因为 $d \ge 2$，所以 $\lfloor n/d \rfloor \le n/2 < n$，故 $0 \le r < 1$。
   级数和为 $\frac{r}{1-r} = \frac{\lfloor n/d \rfloor / n}{1 - \lfloor n/d \rfloor / n} = \frac{\lfloor n/d \rfloor}{n - \lfloor n/d \rfloor}$。
   $$E[L] = 1 - \sum_{d=2}^{n} \mu(d) \frac{\lfloor n/d \rfloor}{n - \lfloor n/d \rfloor}$$

5. 计算求和 (模 p):
   我们需要计算 $S = \sum_{d=2}^{n} \mu(d) \frac{\lfloor n/d \rfloor}{n - \lfloor n/d \rfloor} \pmod p$。
   由于 $n$ 非常大 ($10^{11}$)，直接计算是不可行的。
   我们注意到 $\lfloor n/d \rfloor$ 的取值只有 $O(\sqrt{n})$ 种。这提示我们使用整除分块 (数论分块)。
   设 $v = \lfloor n/l \rfloor$。对于 $d \in [l, r]$，其中 $r = \lfloor n / \lfloor n/l \rfloor \rfloor = \lfloor n/v \rfloor$，$\lfloor n/d \rfloor$ 的值都等于 $v$。
   因此，可以将求和式按块计算：
   $$S = \sum_{\text{块 } [l, r], l \ge 2} \left( \sum_{d=l}^{r} \mu(d) \right) \frac{v}{n - v} \pmod p$$
   $$S = \sum_{\text{块 } [l, r], l \ge 2} (S_\mu(r) - S_\mu(l-1)) \cdot v \cdot (n - v)^{-1} \pmod p$$
   其中 $S_\mu(x) = \sum_{i=1}^x \mu(i)$ 是莫比乌斯函数的前缀和。 $(n - v)^{-1}$ 是 $n-v$ 在模 $p$ 意义下的逆元。

6. 计算莫比乌斯函数前缀和 $S_\mu(x)$:
   当 $x$ 较小（例如 $x \le L = 5 \times 10^6$）时，可以通过线性筛预处理出 $\mu(i)$ 和 $S_\mu(i)$。
   当 $x$ 较大时，需要使用更高级的算法，如 杜教筛 (Motor-Roller Sieve)。
   杜教筛计算 $\sum_{i=1}^n f(i)$ 的原理基于狄利克雷卷积和整除分块。
   对于莫比乌斯函数 $\mu$，我们知道 $\mu * \mathbf{1} = \epsilon$，其中 $\mathbf{1}(n)=1$ 是常数函数，$\epsilon(n) = [n=1]$ 是单位元。
   设 $S_\mu(n) = \sum_{i=1}^n \mu(i)$。
   根据杜教筛的公式（由 $\sum_{i=1}^n (f*g)(i) = \sum_{i=1}^n g(i) F(\lfloor n/i \rfloor)$ 推导）：
   令 $f=\mu, g=\mathbf{1}$。 $f*g=\epsilon$。 $G(n) = \sum_{i=1}^n \mathbf{1}(i) = n$。 $\sum_{i=1}^n \epsilon(i) = 1$ (for $n \ge 1$)。
   $1 = \sum_{i=1}^n \mathbf{1}(i) S_\mu(\lfloor n/i \rfloor)$
   $1 = \mathbf{1}(1) S_\mu(\lfloor n/1 \rfloor) + \sum_{i=2}^n \mathbf{1}(i) S_\mu(\lfloor n/i \rfloor)$
   $1 = S_\mu(n) + \sum_{i=2}^n S_\mu(\lfloor n/i \rfloor)$
   $$S_\mu(n) = 1 - \sum_{i=2}^n S_\mu(\lfloor n/i \rfloor)$$
   这个递归式可以使用整除分块和记忆化来计算 $S_\mu(n)$。
   $S_\mu(n) = 1 - \sum_{\text{块 } [l, r], l \ge 2} (r - l + 1) S_\mu(\lfloor n/l \rfloor)$。
   代码中的 getS 函数正是实现了这个递归计算，结合了线性筛预处理小范围的值和 unordered_map 进行记忆化。其复杂度约为 $O(L + n^{2/3})$。

7. 实现细节:

   • 模运算: 所有计算都在模 $p$ 下进行。注意负数的处理，(x % p + p) % p。
   • 模逆元: 使用快速幂计算 $a^{p-2} \pmod p$ (费马小定理)，因为 $p$ 是质数。代码中提供了 power 和 modInverse 函数。
   • 大数乘法: 中间乘法可能超过 long long 范围（例如 $p \approx 10^{12}$ 时 $p \times p$）。需要使用 __int128 进行中间计算再取模。代码中已使用。
   • 整除分块: 循环 for (i64 l = 2, r; l <= n; l = r + 1) 实现。计算 v = n / l，r = n / v。计算块内 $\mu$ 的和 sum_mu_block = getS(r) - getS(l - 1)。计算项的值并累加到 ans。
   • 最终结果: $E[L] = (1 - S) \pmod p$。代码计算的是 $S’ = \sum_{d=2}^{n} (-\mu(d)) \frac{\lfloor n/d \rfloor}{n - \lfloor n/d \rfloor} \pmod p$ (对应代码中的 ans)，然后最终结果是 $(1 + ans) \pmod p$。这两种形式是等价的。

时间复杂度

• 线性筛预处理: $O(L)$。
• 杜教筛计算 $S_\mu(x)$: 每次调用（非记忆化）需要 $O(\sqrt{x})$ 时间。总复杂度为 $O(N^{2/3})$。
• 整除分块计算最终和: $O(\sqrt{N})$ 次迭代，每次迭代调用两次 getS (可能命中缓存或递归)。总复杂度主要受 getS 影响，约为 $O(N^{2/3})$。
• 总时间复杂度: $O(L + N^{2/3})$。对于 $N=10^{11}, L=5 \times 10^6$，这是可行的。

C++ 代码

#include <bits/stdc++.h> // 引入所有标准库

using namespace std;

// 定义类型别名
using i64 = long long;
using u64 = unsigned long long; // 未使用
using u32 = unsigned; // 未使用

int main() {
    // 加速 IO
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    i64 n, p; // 输入 n 和模数 p
    std::cin >> n >> p;

    // 模快速幂函数 (使用 __int128 防止中间溢出)
    auto power = [&](i64 base, i64 exp) {
        i64 res = 1;
        base %= p;
        while (exp > 0) {
            if (exp % 2 == 1) res = (__int128)res * base % p;
            base = (__int128)base * base % p;
            exp /= 2;
        }
        return res;
    };

    // 模逆元函数 (基于费马小定理)
    auto modInverse = [&](i64 num) {
        // i64 M = p; // M 未使用
        return power(num, p - 2); // 计算 num^(p-2) mod p
    };

    const int L = 5000000; // 线性筛预处理的上界
    vector<int> mu(L + 1);          // 莫比乌斯函数 mu[i]
    vector<int> prime; prime.reserve(L / 5); // 存储素数
    vector<bool> is_prime(L + 1, true); // 标记是否为素数
    vector<i64> Smu(L + 1);        // 莫比乌斯函数前缀和 Smu[i] = sum_{j=1..i} mu[j]
    unordered_map<i64, i64> memo;   // 杜教筛记忆化存储

    // 线性筛函数：预计算 mu 和 Smu
    auto sieve = [&]() {
        is_prime[0] = is_prime[1] = false;
        mu[1] = 1; // 定义 mu[1] = 1
        for (int i = 2; i <= L; ++i) {
            if (is_prime[i]) { // i 是素数
                prime.push_back(i);
                mu[i] = -1; // 素数的 mu 值为 -1
            }
            // 筛合数
            for (int pr : prime) {
                 if ((i64)i * pr > L) break; // 超出范围
                 is_prime[i * pr] = false; // 标记为合数
                 if (i % pr == 0) { // pr 是 i 的因子
                    mu[i * pr] = 0; // 如果含平方因子，mu 值为 0
                    break;          // 每个合数只被其最小质因子筛一次
                 } else { // pr 不是 i 的因子
                    mu[i * pr] = -mu[i]; // mu 是积性函数
                 }
            }
        }
        Smu[0] = 0; // 前缀和 Smu[0] = 0
        // 计算前缀和 (模 p)
        for(int i = 1; i <= L; ++i) {
            // Smu[i] = (Smu[i-1] + mu[i]); // 不需要取模，因为要存原始和
             Smu[i] = (Smu[i-1] + mu[i]); // 原始代码没有取模，但在下面取模。这里存原始值是对的。
        }
    };

    sieve(); // 执行线性筛

    // 杜教筛计算莫比乌斯函数前缀和 S_mu(x)
    // 使用 std::function 方便递归调用
    function<i64(i64)> getS =
        [&](i64 x) -> i64 {
        if (x <= L) { // 如果 x 在预处理范围内
            return Smu[x]; // 直接返回预计算的值
        }
        if (memo.count(x)) { // 如果已经计算过 (记忆化)
            return memo[x]; // 返回缓存的值
        }

        // 递归计算 S_mu(x) = 1 - sum_{i=2..x} S_mu(floor(x/i))
        i64 res = 1; // 对应公式中的 1
        // 使用整除分块计算 sum_{i=2..x} S_mu(floor(x/i))
        for (i64 l = 2, r; l <= x; l = r + 1) {
            i64 v = x / l; // 当前块 floor(x/i) 的值
            r = x / v;     // 当前块的右端点
            // 累减 (块长度 * 块对应的 S_mu 值)
            res -= (r - l + 1) * getS(v); // 递归调用 getS(v)
        }
        // 存储结果并返回
        return memo[x] = res;
    };

    // 计算 E[L] = 1 - sum_{d=2..n} mu(d) * floor(n/d) / (n - floor(n/d)) mod p
    // 代码计算 ans = sum_{d=2..n} (-mu(d)) * floor(n/d) / (n - floor(n/d)) mod p
    i64 ans = 0;
    // 使用整除分块计算求和部分
    for (i64 l = 2, r; l <= n; l = r + 1) {
        i64 v = n / l; // 当前块的 floor(n/d) 值
        if (v == 0) break; // 如果 floor(n/d) 为 0，后续块也为 0，停止
        r = n / v; // 当前块的右端点

        // 计算块内 mu 的和: S_mu(r) - S_mu(l-1)
        i64 sum_mu_block = (getS(r) - getS(l - 1));

        // 计算项中的 floor(n/d) / (n - floor(n/d)) 部分 (模 p)
        i64 term_val = v % p; // floor(n/d) mod p

        i64 denom = (n - v); // 分母 (n - floor(n/d))
        denom = (denom % p + p) % p; // 确保分母在 [0, p-1] 范围内
        i64 term_inv = modInverse(denom); // 计算分母的模逆元

        // 计算项中的 (-mu(d)) 部分 (模 p)
        // sum_mu_block 是块内 mu 的和 (原始值，可能负)
        i64 term_mu_raw = sum_mu_block;
        // term_mu 代表块内 -mu 的和 (模 p)
        i64 term_mu = (p - (term_mu_raw % p) + p) % p; // 计算 (-sum_mu_block) mod p

        // 计算当前块的贡献 (模 p，使用 __int128 防止溢出)
        i64 term = (__int128)term_mu * term_val % p; // (-mu_sum) * v
        term = (__int128)term * term_inv % p;        // * (n - v)^(-1)

        // 累加到 ans
        ans = (ans + term) % p;
    }

    // 最终结果是 1 + ans (模 p)
    i64 final_ans = (1 + ans) % p;
    std::cout << final_ans << "\n";

    return 0; // 程序正常结束
}