CF1139D step to one(数论+概率dp)

原题地址
题意：给一个数m和一个空数组a，每次可以选1~m中任意一个数加入到a中，求最后a数组中所有数的gcd为1时的期望长度；
令dp[i] 表示当前数组的gcd为i时，还期望添加dp[i]个数使得所有数的gcd为1；那么初始条件dp[1] = 0;
状态转移公式：
解释一下，1表示添一个数；求和表示从1~m中选一个数加入数组，期望是dp[gcd(i,j)]，选到的概率是1/m。直接做复杂度O(m^2)，过不了，所以要优化，观察这个式子,gcd(i,j)肯定会出现重复项，比如，当i=4时，j=6和j=8都会让gcd(i,j)=2。从这里入手，假设gcd(i,j)的结果是a1,a2…ax共x项，每一项出现的次数是c1,c2…cx；那么原式就是:

并且a[1~x]全部是i的因子。由于dp[i]未知，所以存在a等于i，这个式子还要化解；

此时a1~ax不存在值为i的数；整理得：

求取val：对于枚举i的因子，对于当前枚举到的因子j(注意跳过j==1)
求它的c值，即求[1~m]中有多少个数k满足gcd(k,i)=j，再化简就是[1 ~ m/j]中有多少个数k和i/j互质。这样枚举i的因子，求取的复杂度变成了O(logn)，总共O(n*logn)。

#include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<map>
#include<bitset>
#include<queue>
#define lson (rt << 1)
#define rson (rt << 1 | 1)
#define lowbit(x) ((x)&(-x))
const int maxn = 1e5+10;
const int maxm = 200 + 10;
const int inf_max = 0x3f3f3f;
const int mod = 1000000007;
const double eps = 1e-4;
using namespace std;
typedef long long ll;
inline int read() {
    int x = 0,f = 1;
    char ch = getchar();
    while(ch<'0' || ch>'9') {if(ch == '-') f = -1;ch = getchar();}
    while(ch>='0' && ch<='9') x = x*10 + ch -'0',ch = getchar();
    return x*f;
}
int m,cnt[maxn],ct,sum,dp[maxn];
vector<int>fac[maxn];
int fast_power(int base,int power) {
    int ret = 1;
    while(power) {
        if(power&1) ret = 1ll*ret*base % mod;
        base = 1ll*base*base % mod;
        power >>= 1;
    }
    return ret;
}
void get_fac(int x) {
    int tt = x;
    for(int i = 2;i*i <= x; ++i) {
        if(x%i == 0) {
            fac[tt].push_back(i);
            while(x%i == 0) x /= i;
        }
    }
    if(x > 1) fac[tt].push_back(x);
}
void dfs(int pos,int use,int cur,int up,int now) {
    if(pos == fac[now].size()) {
        if(use&1) sum -= up/cur;
        else sum += up/cur;
        return ;
    }
    dfs(pos+1,use,cur,up,now);
    dfs(pos+1,use+1,cur*fac[now][pos],up,now);
}
int main()
{
    for(int i = 1;i <= 100000; ++i) get_fac(i);
    m = read();
    dp[1] = 0;
    for(int i = 2; i <= m; ++i) {
        int tot = m/i,al = 0;
        for(int j = 1;j*j <= i ;++j) {
            if(i%j) continue;
            sum = 0;
            dfs(0,0,1,m/j,i/j);
            //printf("1:%d\n",sum);
            al = (al+1ll*sum*dp[j]%mod)%mod;  
            if(i/j != j && j != 1) {
                sum = 0;
                dfs(0,0,1,m/(i/j),j);
                al = (al+1ll*sum*dp[i/j]%mod)%mod;
               // printf("2:%d\n",sum);
            }
        }
        dp[i] = 1ll*(al+m)%mod*fast_power(m-tot,mod-2)%mod;
    }
    int ans = 0;
    for(int i = 1;i <= m; ++i) ans = (ans + dp[i])%mod;
    ans = (1ll*ans*fast_power(m,mod-2)%mod+1)%mod;
    cout<<ans<<endl;
    return 0;
}