Codeforces - 1139D - Steps to One (概率DP+莫比乌斯反演)

蒟蒻数学渣呀，根本不会做。

解法是参考 https://blog.csdn.net/xs18952904/article/details/88785210 这位大佬的。

状态的设计和转移如上面博客一样：dp[i]代表当前序列的gcd为i的期望长度。

那么可以写出状态转移方程：dp[i]=(1+(x/m)∑(j|i,j≠i)dp[j]) / (1-(m/i)/m) (写得有点乱，其实和上面大佬的一样的)

这里要说一下的是 x=∑(t=1,t<=m) [ gcd(t,i)==j ] 就是怎么求1<=t<=m 中gcd(t,i)=j的t的个数。

这里考虑莫比乌斯反演：

x=∑(t=1,t<=m)[gcd(t,i)=j]

把j提出来 x=∑(t=1,t<=m/j) [gcd(t,i/j)=1]

代入莫比乌斯性质：x=∑(t=1,t<=m) ∑(d|gcd(t,i/j)) μ(d)

套路，改为枚举d : x=(m/jd)(d|(i/j)) μ(d)

这样就可以求出x了，这道题就可以解决了。

时间复杂度为O(m*log(m)*因子个数)

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=1e5+;

const int MOD=1e9+;

LL m,mu[N],v[N],dp[N];

LL power(LL x,LL p) {

    LL ret=;

    for (;p;p>>=) {

        if (p&) ret=(ret*x)%MOD;

        x=(x*x)%MOD;

    }

    return ret;

}

vector<int> fac[N];

void prework() {

    for (int i=;i<=m;i++)

        for (int j=;j<=m/i;j++)

            fac[i*j].push_back(i);

    for (int i=;i<=m;i++) mu[i]=,v[i]=;

    for (int i=;i<=m;i++) {

        if (v[i]) continue;

        mu[i]=-;

        for (int j=*i;j<=m;j+=i) {

            v[j]=;

            if ((j/i)%i==) mu[j]=;

            else mu[j]*=-;

        }

    }

}

LL calc(LL i,LL j) {

    LL ret=;

    for (int k=;k<fac[i/j].size();k++) {

        int d=fac[i/j][k];

        LL tmp=(mu[d]+MOD)%MOD*(m/j/d)%MOD;

        ret=(ret+tmp)%MOD;

    }

    return ret;

}

int main()

{

    cin>>m;

    prework();

    LL ans=; dp[]=;

    for (int i=;i<=m;i++) {

        dp[i]=;

        for (int j=;j<fac[i].size();j++) {

            if (fac[i][j]==i) continue;

            LL x=calc(i,fac[i][j]);

            dp[i]=(dp[i]+x*dp[fac[i][j]]%MOD)%MOD;

        }

        dp[i]=dp[i]*power(m,MOD-)%MOD;

        dp[i]=(dp[i]+)%MOD;

        dp[i]=(dp[i]*m%MOD*power(m-m/i,MOD-))%MOD;

        ans=(ans+dp[i]*power(m,MOD-)%MOD)%MOD;

    }

    cout<<ans<<endl;

    return ;

}

巴特西

Codeforces - 1139D - Steps to One (概率DP+莫比乌斯反演)

最新文章

热门文章