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

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

SRE实战 互联网时代守护先锋，助力企业售后服务体系运筹帷幄！一键直达领取阿里云限量特价优惠。

状态的设计和转移如上面博客一样：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+10;
const int MOD=1e9+7;
LL m,mu[N],v[N],dp[N];

LL power(LL x,LL p) {
    LL ret=1;
    for (;p;p>>=1) {
        if (p&1) ret=(ret*x)%MOD;
        x=(x*x)%MOD;
    }
    return ret;
}

vector<int> fac[N];
void prework() {
    for (int i=1;i<=m;i++)
        for (int j=1;j<=m/i;j++)
            fac[i*j].push_back(i);        
    for (int i=1;i<=m;i++) mu[i]=1,v[i]=0;
    for (int i=2;i<=m;i++) {
        if (v[i]) continue;
        mu[i]=-1;
        for (int j=2*i;j<=m;j+=i) {
            v[j]=1;
            if ((j/i)%i==0) mu[j]=0;
            else mu[j]*=-1;
        }
    }
}

LL calc(LL i,LL j) {
    LL ret=0;
    for (int k=0;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=1; dp[1]=0;
    for (int i=2;i<=m;i++) {
        dp[i]=0;
        for (int j=0;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-2)%MOD;
        dp[i]=(dp[i]+1)%MOD;
        dp[i]=(dp[i]*m%MOD*power(m-m/i,MOD-2))%MOD;
        ans=(ans+dp[i]*power(m,MOD-2)%MOD)%MOD;
    }
    cout<<ans<<endl;
    return 0;
}

 

扫码关注我们
微信号：SRE实战
拒绝背锅 运筹帷幄