CF1139D Steps to One

发表于 | 阅读数

给一个数列，每次随机选一个 $1$ 到 $m$之间的数加在数列末尾，数列中所有数的 $\gcd=1$ 时停止，求期望长度。

可以先把第一个加的数字确定。

$f[i]$表示从当前$\gcd=i$,变到$1$需要的期望步数。
$f[1]=0$

$f[i]=1+\frac{f[\gcd(i,j)]}{m}$

分析下复杂度$\sum_{i=1}^m\sum_{j=1}^{\frac{m}{i}}D(j)=\sum_{i=1}^m\sum_{d=1}^{\frac{m}{i}}\frac{\frac{m}{i}}{d}=\sum_{i=1}^m\frac{m}{i}\log m=O(m\log^2 m)$

代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97


#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define pii pair<int, int>
const int N = 1e6 + 10;
const int mod = 1e9 + 7;
int inc(int x, int y) { return (x + y >= mod) ? (x + y - mod) : (x + y); }
int del(int x, int y) { return (x - y < 0) ? (x - y + mod) : (x - y); }

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

int pcnt, pr[N], npr[N];
int mu[N], usum[N];

void Prime_init(int X)
{
    npr[1] = 1;
    mu[1] = 1;
    for (int i = 1; i <= X; i++)
    {
        if (!npr[i])
            pr[++pcnt] = i, mu[i] = -1;
        for (int j = 1; j <= pcnt && pr[j] * i <= X; j++)
        {
            npr[pr[j] * i] = 1;

            if (i % pr[j] == 0)
            {
                mu[pr[j] * i] = 0;
                break;
            }
            else
            {
                mu[pr[j] * i] = mu[i] * (-1);
            }
        }
    }
}
int main()
{
    int m;
    cin >> m;
    vector<vector<int>> g(m + 1);
    for (int i = 1; i <= m; i++)
        for (int j = i; j <= m; j += i)
            g[j].push_back(i);
    Prime_init(m);
    auto calc = [&](int a, int b) {
        int res = 0;
        for (int d : g[b])
        {
            res += mu[d] * (a / d);
        }
        //cout << a << " " << b << " " << res << endl;
        return res;
    };
    vector<vector<int>> z(m + 1);
    for (int i = 1, j; i <= m; i = j + 1)
    {
        j = m / (m / i);
        z[m / i].resize(m / i + 1);
        for (int k = 1; k <= m / i; k++)
        {
            z[m / i][k] = calc(m / i, k);
        }
    }
    vector<int> f(m + 1);
    int ans = 0;
    for (int i = 2; i <= m; i++)
    {
        int res = 0;
        for (int d : g[i])
        {
            if (d == i)
                continue;

            res = inc(res, 1ll * f[d] * z[m / d][i / d] % mod);
        }
        f[i] = 1ll * (res + m) * qpow(m - m / i, mod - 2, mod) % mod;
        //cout << f[i] << endl;
        ans = inc(ans, f[i]);
    }
    cout << 1ll * (ans + m) * qpow(m, mod - 2, mod) % mod << endl;
}