# Codeforces Round 548 (Div. 2) (D)
> https://codeforces.com/contest/1139/problem/D

給一個整數$m$，現在有個空序列，每次操作為隨機從$[1,m]$中取一數加進序列中，如果整個序列之$gcd$等於1則停下，求序列期望長度模$mod$。$1 \le m \le 10^5$。

這題狀態還蠻好想的，定義$dp(i)$為目前序列之$gcd$為$i$時執行操作到結束的期望長度，假如$i=6$則我們可以寫出：
\begin{array}{ll}
dp(6) = 1 + \frac{1}{m} \left[
F(6, 1) \times dp(1) +
F(6, 2) \times dp(2) + \\
F(6, 3) \times dp(3) +
F(6, 6) \times dp(6)
\right],
\end{array}
其中$F(i, j)$為1到m中有幾個數使得其和$i$的$gcd$等於$j$，$j$必須是$i$的因數，於是有：
\begin{array}{ll}
F(6, 6) = \left \lfloor \frac{m}{6} \right \rfloor, \\
F(6, 3) = \left \lfloor \frac{m}{3} \right \rfloor - \left \lfloor \frac{m}{6} \right \rfloor, \\
F(6, 2) = \left \lfloor \frac{m}{2} \right \rfloor - \left \lfloor \frac{m}{6} \right \rfloor, \\
F(6, 1) = m - \left \lfloor \frac{m}{2} \right \rfloor - \left \lfloor \frac{m}{3} \right \rfloor + \left \lfloor \frac{m}{6} \right \rfloor,
\end{array}
其實就是對$i/j$的質因數做容斥，如果要算整個$dp(i)$的話，我們可以考慮枚舉每個$i$的因數$j$，對$dp(i)$的貢獻有$\left \lfloor \frac{m}{j} \right \rfloor dp(j/j的質因數乘積)$這東西，如果是奇數個質因數就扣掉，否則就加，這邊複雜度應為$\sigma_0(i)2^6 \times 6$，$\sigma_0(i)$為$i$之因數個數，然後$\le 10^5$最多質因數有6個。

假如目前為空序列，則我們操作一次，假如說我們得到x則表示目前序列之$gcd$等於x，於是我們可以得到所求應為：
\begin{array}{ll}
\sum_{k=1}^{m} \frac{dp(k)}{m},
\end{array}
其中$dp(1) \equiv 1$。

如何求$dp(i)$呢?列一下式子：
\begin{array}{ll}
&dp(i) = 1 + \frac{1}{m} \sum_{d|i}F(i, d)dp(d) \\
\Rightarrow &dp(i) (1 - \frac{F(i, i)}{m}) = 1 + \frac{1}{m} \sum_{d|i, d<i} F(i, d) dp(d) \\
\Rightarrow &dp(i)(m - F(i, i)) = m + \sum_{d|i, d<i} F(i, d) dp(d),
\end{array}
預處理每個數的因數和質因數，於是1到m中因數個數和有$mlog(m)$個複雜度應該過得去。

#include <bits/stdc++.h>
#pragma GCC optimize(3)
#define ll long long
#define pii pair<int, int>
#define pll pair<long long, long long>
#define F first
#define S second
#define endl '\n'
using namespace std;
const int inf = 0x3f3f3f3f;
const ll Inf = 1e18;
const int mod = 1e9 + 7;
const int N = 1e5 + 5;

vector<vector<int>> d(N), pd(N);

const int maxn = 1e5+5;
vector<bool> p(maxn, 1);
vector<int> prime;

ll inv[maxn];

void init(){
    p[0] = p[1] = 0;
    for(int i=2; i<maxn; i++){
        if(p[i]){
            prime.push_back(i);
            for(ll j=1LL*i*i; j<1LL*maxn; j+=i){
                p[j] = 0;
            }
        }
    }

    inv[0] = inv[1] = 1;
    for(int i=2; i<maxn; i++){
        inv[i] = mod - mod/i * inv[mod%i] % mod;
    }

    return;
}

void solve(){
    int m; cin >> m;
    for(int i=1; i<=m; i++){
        for(int j=i; j<=m; j+=i){
            d[j].push_back(i);
        }
    }

    init();
    for(auto i : prime){
        for(int j=i; j<=m; j+=i){
            pd[j].push_back(i);
        }
    }

    ll dp[m+1];
    dp[1] = 1;
    
    for(int i=2; i<=m; i++){
        ll tmp = m;
        for(auto j : d[i]){
            int l = pd[j].size();
            for(int k=0; k<(1<<l); k++){
                int cur = 1;
                for(int kk=0; kk<l; kk++){
                    if(~(k>>kk) & 1) continue;
                    cur *= pd[j][kk];
                }
                if(j == i && cur == 1) continue;
                tmp = (tmp + (__builtin_popcount(k) & 1 ? -1LL : 1LL) * (m/j) * dp[j/cur]) % mod;
            }
        }
        if(tmp < 0) tmp += mod;
        dp[i] = inv[m - m/i] * tmp % mod;
    }

    ll ans = 0;
    for(int i=1; i<=m; i++) ans = (ans + dp[i]) % mod;
    ans = ans * inv[m] % mod;
    cout << ans << endl;

} 
 
 
int main(){
    ios::sync_with_stdio(0); cin.tie(0);
    //freopen("input.txt", "r", stdin);
    //freopen("output.txt", "w", stdout);

    int t = 1;
    // cin >> t;
    while(t--)
    solve();
    return 0;
}

這題還有莫比烏斯函數的解法，可以參考：
> https://codeforces.com/blog/entry/66101?#comment-501072
> https://codeforces.com/blog/entry/66101?#comment-500999

### 時間複雜度 ：$O(mlog(m)) \times 2^6 \times 6$
###### tags: `number theory` `math` `dp`