CodeForces 1139D Steps to One（概率dp 容斥/莫比乌斯反演）

题目链接https://codeforces.com/contest/1139/problem/D

题意：给定一个m，每次在1-m中随机取一个数放到容器中，当容器的gcd为1时停止，求期望步数，用分数形式的逆元输出

 

题解思路：

表示只会容斥的思路啊，反演的题解那步概率的转移没看懂= =(p/(1-p))哪个不知道怎么来的

下面主要提一下容斥，以后有能力了补莫比乌斯反演的方法。

 

设dp[x]表示当前gcd为x，转到gcd为1的期望步数

 

这里的cnt表示在中取数p满足的数量

这个地方求就需要容斥一下了

因为gcd的变化只能往小于等于的方向变，如果取到,满足，这时新的gcd还是x，需要在容斥完之后解一下方程

容斥的过程假设

在不考虑去重的情况下，且 

我们对x和y打质因数分解，假设当前x有的一个因子，y有的一个因子，一定有

假如那么p在y的基础上就不能取这个

那么就不能再有，遍历的时候这个地方奇减偶加进行容斥

 

最后的，加一是第一步取什么数就转到哪个状态

 

#include <bits/stdc++.h>
#define quickio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
#define ll long long
const int mod=1e9+7;
const int N=1e5+5;
const int maxn=1e5;
ll dp[maxn+5];
bool is_prime[N+5];
vector<int>factor[maxn+5];
map<int,int>ma[N+5];


ll quick_pow_mod(ll a,ll b,ll c)
{
    ll res=1;
    while(b)
    {
        if(b & 1)
        {
            res=(res*a)%c;
        }
        a=(a*a)%c;
        b=b>>1;
    }
    return res;
}

void find_fac()
{
    memset(is_prime,true,sizeof(is_prime));
    for(int i=2;i<=maxn;i++)
    {
        if(is_prime[i]==1)
        {
            for(int j=i;j<=maxn;j+=i)
            {
                int cur=j;
                int cnt=0;
                while(cur%i==0)
                {
                    cur/=i;
                    cnt++;
                }
                ma[j][i]=cnt;
                is_prime[j]=false;
            }
        }
    }

    for (int i = 1; i < N; i++)
    {
        for (int j = 2 * i; j < N; j += i)
            factor[j].push_back(i);
    }
}
int m;
ll find(ll x,ll n)//gcd(n,p)=x的p的数量
{
    vector<int>a;//p里不能有的质因子
    for(auto &it:ma[n])
    {
        if(!ma[x].count(it.first))//如果x中没有n的因子
        {
            a.push_back(it.first);
            continue;
        }
        if(it.second==ma[x][it.first])continue;
        else
            a.push_back((it.first));
    }
    ll sz =(int)a.size();
    ll lim=m/x;
    ll ans=m/x;//初始不考虑重复的情况下
    for(ll i=1;i<(1<<sz);i++)
    {
        ll now=1;
        int cnt=0;
        for(int j=0;j<sz;j++)
        {
            if(i&(1<<j))
            {
                now*=a[j];
                cnt++;
            }
        }
        if(cnt%2==1)ans-=lim/now;
        else ans+=lim/now;
    }
    return ans;
}

int main()
{
    quickio
    find_fac();//打每个m对应的质因数分解表
    dp[1]=0;
    cin>>m;
    ll invm=quick_pow_mod(m,mod-2,mod);
    for(int i=2;i<=m;i++)
    {
        ll a=1;
        for(int x:factor[i])
        {
            ll cnt=find(x,i);//满足new_gcd=x的数量
            a+=(dp[x]*cnt)%mod*invm%mod;
            if(a>mod)a-=mod;
        }
        ll cnt=m/i;
        dp[i]=a*m%mod*quick_pow_mod(m-cnt,mod-2,mod)%mod;
        //最后这一步的目的是解方程，gcd为i如果取到p,i|p,新的gcd仍为i
    }
    ll res=1;//第一步
    for(int i=1;i<=m;i++)
    {
        res+=dp[i]*invm%mod;
        if(res>=mod)res-=mod;
    }
    cout<<res<<endl;

    
    return 0;
}

 

2019.3.31更新莫比乌斯反演的方法

现在定义一个E(i)，表示有的所有序列的期望长度 

 

P(随机取一个数满足)=

对于gcd为i的序列，长度为1的序列的期望长度为

对于单独求gcd为i的时候的长度为x的序列的期望长度，就是 ，但是这样计算会存在重复，对于长度为1的序列，他有可能在下次随机取数的操作中仍旧满足，所以对于每一个x，我们只需要累加到无穷，即可求得

而这个满足等比数列的求和公式，当项数趋于无穷的时候，=

下面是莫比乌斯反演的内容：

得到答案需要计算的是，也就是我们这里定义的

答案需要的E和我们上面设的E的关系就是一个容斥的关系，在不考虑重复的情况下，=

但是对于4来说，前面2已经出现了一次，再加上自身本身，对答案的贡献的系数=-1+1=0

由莫比乌斯反演的知识可以得到，在计算的时候，对每个需要乘上一个 ，得到的就是结果了

最后再+1，是由变到的一步

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int maxn=1e5;
int pri[maxn+5];
int mu[maxn+5];
bool check[maxn+5];
int m;
const int mod=1e9+7;

ll quick_pow_mod(ll a,ll b,ll c)
{
    ll res=1;
    while(b)
    {
        if(b & 1)
        {
            res=(res*a)%c;
        }
        a=(a*a)%c;
        b=b>>1;
    }
    return res;
}
void init()
{
    memset(check,0,sizeof(check));
    mu[1]=1;
    int tot=0;
    for(int i=2;i<=maxn;i++)
    {
        if(!check[i])
        {
            pri[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot;j++)
        {
            if(i*pri[j]>maxn)break;
            check[i*pri[j]]=1;
            if(i%pri[j]==0)
            {
                mu[i*pri[j]]=0;
                break;
            }
            else
            {
                mu[i*pri[j]]=-mu[i];
            }
        }
    }
}
int main()
{
    cin>>m;
    //ll invm=quick_pow_mod(m,mod-2,mod);
    init();

    ll res=1;
    for(int i=2;i<=m;i++)
    {
        res=res-mu[i]*m/i*quick_pow_mod(m-m/i,mod-2,mod);
        res=(res+mod)%mod;
    }
    cout<<res<<endl;

    return 0;
}