题解 CF1139D 【Steps to One】

xyz32768 · 2019-03-29 11:06:52 · 题解

题意

• 给定正整数 m ，你有一个空集合 a

• 每次你可以等概率随机一个 [1,m] 内的正整数加入集合 a

• 加入之后，如果集合 a 内所有数的 \gcd 为 1 则结束过程，否则继续随机

• 求过程结束后集合 a 的期望大小

• 1\le m\le 100000

算法：莫比乌斯反演 + 期望 DP

• 状态 f[i] 表示当前集合 a 的 \gcd 为 i 的情况下，期望加入多少个数之后这个 \gcd 变成 1

• 显然有

• f[1]=0
• f[i]=1+\frac 1m\sum_{j=1}^mf[\gcd(i,j)]

• 发现 \gcd(i,j) 的取值个数只有 i 的约数个数，所以把枚举 j 改成枚举 \gcd(i,j)

• f[i]+1+\frac 1m\sum_{d|i}f[d]\times cnt(d,i)
• cnt(d,i)$ 表示 $[1,m]$ 内有多少个数与 $i$ 的 $\gcd$ 等于 $d

• 另设 g(i,j) 表示 [1,i] 内有多少个数与 j 互质

• 显然有 cnt(d,i)=g(\lfloor\frac md\rfloor,\frac id)

• 注意到 \sum_{i=1}^m\frac mi=O(m\log m)

• 所以考虑求出所有的 ij\le m ，对应 g(\lfloor\frac mi\rfloor,j) 的值

• 那么对于给定的 a\ge b ，如何求有多少个 1\le i\le a 满足 \gcd(i,b)=1 呢？

• 如果不是满足 \gcd(i,b)=1 而是满足 x|\gcd(i,b) ，则这显然为 \lfloor\frac ax\rfloor

• 所以，考虑反演，得

• g(a,b)=\sum_{x|b}\mu(x)\lfloor\frac ax\rfloor

• 设 D(x) 为 x 的约数个数

• 那么线性筛 \mu 之后就可以在 O(D(b)) 的时间内求出 g(a,b)

• 对每对 ij\le m 进行计算的总复杂度为

• \sum_{i=1}^m\sum_{j=1}^{\lfloor\frac mi\rfloor}D(j)=\sum_{i=1}^m\sum_{j=1}^{\lfloor\frac mi\rfloor}\lfloor\frac{\lfloor\frac mi\rfloor}j\rfloor=\sum_{i=1}^m\lfloor\frac mi\rfloor\log m=O(m\log^2m)

• 回到这个转移

• f[i]+1+\frac 1m\sum_{d|i}f[d]\times cnt(d,i)

• 发现 d 可以等于 i ，得不到转移的目的，所以做一些移项

• (1-\frac{cnt(i,i)}m)f[i]=1+\frac1m\sum_{d|i,d\ne i}f[d]\times cnt(d,i)
• f[i]=\frac1{1-\frac{cnt(i,i)}m}(1+\frac1m\sum_{d|i,d\ne i}f[d]\times cnt(d,i))

• 于是我们就可以转移了

• 考虑如何算答案

• 由于第一个数可以随便取，所以答案为

• \frac1m\sum_{i=1}^m(1+f[i])

• 看到某人的提交记录里有吊打标算的非 DP 做法，复杂度仿佛是 O(m\log m) ，所以有更优做法欢迎提出

代码

• 该代码中的 f[i] 表示的是题解中的 f[i] 加一，答案为 \frac1m\sum_{i=1}^mf[i] ，边界为 f[1]=1 ，见谅

#include <bits/stdc++.h>

// 20030830

inline int read()
{
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}

const int N = 1e5 + 5, ZZQ = 1e9 + 7;

int qpow(int a, int b)
{
    int res = 1;
    while (b)
    {
        if (b & 1) res = 1ll * res * a % ZZQ;
        a = 1ll * a * a % ZZQ;
        b >>= 1;
    }
    return res;
}

int n, f[N], miu[N], tot, pri[N], inv[N], wt[N], num[N], res[N], ans;
bool mark[N];
std::vector<int> divi[N], cont[N];

int solve(int n, int m)
{
    int tt = divi[m].size(), res = 0;
    for (int i = 0; i < tt; i++)
    {
        int x = divi[m][i];
        res += miu[x] * (n / x);
    }
    return res;
}

void sieve()
{
    mark[0] = mark[miu[1] = 1] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!mark[i]) miu[pri[++tot] = i] = -1;
        for (int j = 1; j <= tot; j++)
        {
            if (1ll * i * pri[j] > n) break;
            mark[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
            else miu[i * pri[j]] = -miu[i];
        }
    }
    for (int i = 1; i <= n; i++)
        for (int j = i; j <= n; j += i)
            divi[j].push_back(i);
    for (int i = 1; i <= n;)
    {
        int nxt = n / (n / i);
        for (int j = 1; j <= n / i; j++) res[j] = solve(n / i, j);
        for (int j = i; j <= nxt; j++)
            for (int k = j; k <= n; k += j)
                cont[k].push_back(res[k / j]);
        i = nxt + 1;
    }
}

int main()
{
    n = read();
    inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = 1ll * (ZZQ - ZZQ / i) * inv[ZZQ % i] % ZZQ;
    sieve();
    f[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        int tt = divi[i].size();
        for (int j = 0; j < tt; j++)
            wt[j + 1] = 1ll * inv[n] * cont[i][j] % ZZQ, num[j + 1] = divi[i][j];
        int qt = qpow((1 - wt[tt] + ZZQ) % ZZQ, ZZQ - 2);
        for (int j = 1; j < tt; j++) wt[j] = 1ll * wt[j] * qt % ZZQ;
        f[i] = qt;
        for (int j = 1; j < tt; j++)
            f[i] = (1ll * wt[j] * f[num[j]] + f[i]) % ZZQ;
    }
    for (int i = 1; i <= n; i++)
        ans = (1ll * inv[n] * f[i] + ans) % ZZQ;
    std::cout << ans << std::endl;
    return 0;
}