随机数列期望长度计算与杜教筛优化实践

1. 问题背景与数学建模

这道题目描述了一个随机数列生成过程，并需要我们计算该数列的期望长度。让我们先拆解题目描述：

1. 初始时数列A为空
2. 每次从[1,n]中等概率随机选取一个整数加入A
3. 检查是否存在大于1的整数w，使得A中所有元素都是w的倍数
   • 如果存在这样的w，则继续步骤2
   • 否则，返回当前A作为结果

这个过程的终止条件是数列A中的元素不全部是任何大于1的整数的倍数。换句话说，终止条件是数列中所有元素的最大公约数(GCD)为1。

1.1 概率模型分析

设E(n)为数列A的期望长度。考虑第一次选择的数i：

• 如果i=1，那么GCD立即变为1，过程终止，此时长度为1
• 如果i>1，我们需要考虑后续选择的数必须与i互质（否则GCD会大于1）

这引导我们想到使用容斥原理和莫比乌斯函数。实际上，经过推导可以发现：

E(n) = 1 + Σ (μ(d) / (n - ⌊n/d⌋)) ，其中d从2到n，μ是莫比乌斯函数

这个公式的推导基于以下观察：

1. 数列长度可以看作是一个几何分布的变种
2. 每次失败（GCD>1）相当于所有元素都是某个d>1的倍数
3. 莫比乌斯函数帮助我们处理容斥关系

2. 算法设计与优化

2.1 杜教筛基础

题目要求我们高效计算莫比乌斯函数的前缀和S(n)=Σμ(i)。对于大n（1e11），我们需要使用杜教筛。

杜教筛的核心思想是：
S(n) = 1 - Σ S(⌊n/i⌋) ，i从2到n

直接计算这个递归式的时间复杂度是O(n^(3/4))，通过预处理前n^(2/3)项可以将复杂度优化到O(n^(2/3))。

2.2 分段数论分块优化

传统数论分块在计算ΣS(n/i)时，对于i较小的情况效率不高。我们采用两阶段优化：

1. 直接处理i ≤ √n的情况，因为此时n/i变化剧烈
2. 对于i > √n的情况，通过枚举v=n/i的值来反推i的范围

这种优化避免了大量64位除法操作，显著提升了性能。

2.3 批量逆元计算

题目需要对多个(n - v_i)计算模逆元。传统方法是逐个计算，时间复杂度为O(k log p)。我们可以使用以下优化：

1. 计算前缀积数组pre[i] = Π(n - v_j) for j ≤ i
2. 计算总逆元inv_total = pre[k]^(-1)
3. 逆向遍历计算每个逆元：inv[i] = inv_total * pre[i-1]

这样将k次逆元计算优化为O(k + log p)时间。

3. 代码实现详解

3.1 预处理阶段

java复制static void precompute() {
    sumMu = new int[MAX + 1];
    primes = new int[MAX / 10 + 100];
    int[] isNotPrime = new int[(MAX >> 5) + 1];

    int cnt = 0;
    sumMu[1] = 1;
    for (int i = 2; i <= MAX; i++) {
        if ((isNotPrime[i >> 5] & (1 << (i & 31))) == 0) {
            primes[cnt++] = i;
            sumMu[i] = -1; 
        }
        for (int j = 0; j < cnt; j++) {
            int p = primes[j];
            int next = i * p;
            if (next > MAX) break;
            isNotPrime[next >> 5] |= (1 << (next & 31)); 
            if (i % p == 0) {
                sumMu[next] = 0;
                break;
            }
            sumMu[next] = -sumMu[i];
        }
    }
    for (int i = 1; i <= MAX; i++) sumMu[i] += sumMu[i - 1];
}

这段代码实现了线性筛法预处理莫比乌斯函数：

1. 使用位数组isNotPrime压缩存储质数标记
2. 同时计算莫比乌斯函数μ(i)及其前缀和
3. 时间复杂度O(MAX)，空间复杂度O(MAX)

3.2 杜教筛实现

java复制static long getSumMu(long n) {
    if (n <= MAX) return sumMu[(int)n] % P;
    
    int idx = (int)(N / n);
    if (visitedLarge[idx]) return memoLarge[idx];

    long res = 1;
    long sqrtn = (long) Math.sqrt(n);
    
    // 直接处理小i
    for (long l = 2; l <= sqrtn; l++) {
        res -= getSumMu(n / l);
    }
    res %= P;
    
    // 通过v=n/i的值处理大i
    long lastR = sqrtn;
    long vStart = n / (sqrtn + 1);
    for (long v = vStart; v >= 1; v--) {
        long r = n / v;
        long len = r - lastR;
        if (len > 0) {
            long nextS = sumMu[(int)v] % P; 
            if (len == 1) res -= nextS;
            else res -= mul(len % P, nextS);
            lastR = r;
        }
    }
    
    res %= P;
    if (res < 0) res += P;
    visitedLarge[idx] = true;
    memoLarge[idx] = res;
    return res;
}

这段代码实现了优化的杜教筛：

1. 对小规模n直接返回预处理结果
2. 使用记忆化存储避免重复计算
3. 采用两阶段分块优化计算效率

3.3 主逻辑流程

java复制public static void main(String[] args) throws Exception {
    // 输入处理
    BufferedReader br = new BufferedReader(new InputStreamReader(System.in));
    StringTokenizer st = new StringTokenizer(br.readLine());
    N = Long.parseLong(st.nextToken());
    P = Long.parseLong(st.nextToken());

    // 预处理
    if (N < MAX) MAX = (int) N;
    precompute();
    
    // 杜教筛初始化
    int threshold = (int)(N / MAX);
    memoLarge = new long[threshold + 1];
    visitedLarge = new boolean[threshold + 1];

    // 数论分块收集项
    long sqrtN = (long) Math.sqrt(N);
    int maxBlocks = (int)(sqrtN * 2 + 10);
    long[] vals = new long[maxBlocks];    
    long[] weights = new long[maxBlocks]; 
    int blockCnt = 0;

    // 阶段1：直接处理小i
    for (int l = 2; l <= sqrtN; l++) {
        long val = N / l;
        long muL = sumMu[l] - sumMu[l - 1]; 
        if (muL != 0) {
            vals[blockCnt] = val;
            weights[blockCnt] = muL;
            blockCnt++;
        }
    }

    // 阶段2：通过v处理大i
    long lastR = sqrtN;
    long vStart = N / (sqrtN + 1);
    long prevSum = getSumMu(lastR); 
    for (long v = vStart; v >= 1; v--) {
        long r = N / v; 
        if (r > lastR) {
            long currSum = getSumMu(r);
            long sMu = currSum - prevSum; 
            if (sMu != 0) {
                vals[blockCnt] = v;
                weights[blockCnt] = sMu;
                blockCnt++;
            }
            prevSum = currSum;
            lastR = r;
        }
    }

    // 批量逆元计算
    long[] invs = new long[blockCnt];
    long[] pre = new long[blockCnt];
    pre[0] = (N - vals[0]) % P;
    for (int i = 1; i < blockCnt; i++) {
        pre[i] = mul(pre[i - 1], (N - vals[i]) % P);
    }

    long totalInv = modInverse(pre[blockCnt - 1]);
    for (int i = blockCnt - 1; i > 0; i--) {
        invs[i] = mul(totalInv, pre[i - 1]);
        totalInv = mul(totalInv, (N - vals[i]) % P);
    }
    invs[0] = totalInv;

    // 计算最终答案
    long ans = 1;
    for (int i = 0; i < blockCnt; i++) {
        long term = mul(vals[i] % P, invs[i]);
        ans -= mul(weights[i], term);
    }
    ans %= P;
    if (ans < 0) ans += P;
    System.out.println(ans);
}

4. 关键优化技巧

4.1 快速模乘实现

java复制static long mul(long a, long b) {
    long res = a * b - (long)((double)a * b / P) * P;
    while (res < 0) res += P;
    while (res >= P) res -= P;
    return res;
}

这个实现避免了直接计算a*b可能导致的溢出，通过浮点估算商值来高效计算模乘。

4.2 记忆化存储设计

对于杜教筛的大n结果，我们使用以下方式存储：

java复制static long[] memoLarge;
static boolean[] visitedLarge;

// 索引计算
int idx = (int)(N / n);
memoLarge[idx] = res;
visitedLarge[idx] = true;

这种设计将n映射到紧凑的数组索引，节省内存空间。

4.3 分块策略选择

代码中根据i的大小采用不同处理策略：

1. i ≤ √n：直接遍历，因为n/i变化剧烈
2. i > √n：通过v=n/i的值连续处理，避免除法

这种混合策略在实践中比纯分块更高效。

5. 复杂度分析

1. 预处理阶段：O(MAX)时间和空间
2. 杜教筛阶段：O(n^(2/3))时间
3. 数论分块：O(√n)时间
4. 批量逆元：O(k + log p)时间

总体时间复杂度主要由杜教筛决定，为O(n^(2/3))，可以处理n=1e11的情况。

6. 实际测试与调优

在实现过程中，有几个关键点需要注意：

1. 预处理大小MAX的选择：通常取n^(2/3)左右，在Java中要考虑内存限制
2. 分块阈值的选择：√n是一个经验值，可以根据实际测试调整
3. 模运算的优化：避免频繁取模，但要注意中间结果溢出

经过测试，这个实现可以在PTA平台上在时间限制内完成最大规模的计算。