Codeforces.1139D.Steps to One(DP 莫比乌斯反演)

题目链接

啊啊啊我在干什么啊。怎么这么颓一道题做这么久。。
又记错莫比乌斯反演式子了(╯‵□′)╯︵┻━┻


\(Description\)

给定\(n\)。有一个初始为空的集合\(S\)。令\(g\)表示S中所有数的\(\gcd\)。每次随机选择一个\([1,n]\)中的数加到集合\(S\)中去,直到\(g=1\)。求集合\(S\)的期望大小。(原题目描述为数列长度,\(n\)是指\(m\)我自己都看混了=-=
\(n\leq10^5\)

\(Solution\)

首先不要想\(f[i][j]\)这种奇奇怪怪的二维DP状态。。
\(f[x]\)表示\(g=x\)时集合的期望大小。那么有\[\begin{aligned}f[x]&=1+\frac1n\sum_{i=1}^nf[\gcd(i,x)]\\&=1+\frac1n\sum_{d\mid x}f[d]\sum_{i=1}^n[(i,x)=d]\\&=1+\frac1n\sum_{d\mid x}f[d]\sum_{i=1}^{\lfloor\frac nd\rfloor}[(i,\lfloor\frac xd\rfloor)=1]\end{aligned}\]

问题在于求\(\sum_{i=1}^{\lfloor\frac nd\rfloor}[(i,\lfloor\frac xd\rfloor)=1]\)。不妨设\(g(n,k)=\sum_{i=1}^n[(i,k)=1]\)
若限制不是\([(i,k)=1]\)而是\([x\mid (i,k)]\),显然答案就是\([x\mid k]\lfloor\frac nx\rfloor\)
所以令\(f(d)=\sum_{i=1}^n[(i,k)=d]\)\(F(d)=\sum_{i=1}^n[d\mid(i,k)]=[d\mid k]\lfloor\frac nd\rfloor\),有\[\begin{aligned}g(n,k)=f(1)&=\sum_{i=1}^n\mu(i)F(i)\\&=\sum_{i=1}^n\mu(i)[i\mid k]\lfloor\frac ni\rfloor\\&=\sum_{d\mid k}\mu(d)\lfloor\frac nd\rfloor\end{aligned}\]

那么直接这样计算一次\(g(n,k)\)的复杂度是\(O(约数个数)\)的。\(1\sim n\)的约数个数均摊应该是\(O(\log n)\)的?

发现转移的时候两边都有\(f[x]\)。移下项,\[(1-\frac1n\sum_{i=1}^{\lfloor\frac nx\rfloor}[(i,1)=1])f[x]=1+\frac1n\sum_{d\mid x,d\neq x}f[d]\sum_{i=1}^{\lfloor\frac nd\rfloor}[(i,\lfloor\frac xd\rfloor)=1]\\f[x]=\frac{1+\frac1n\sum_{d\mid x,d\neq x}f[d]\sum_{i=1}^{\lfloor\frac nd\rfloor}[(i,\lfloor\frac xd\rfloor)=1]}{1-\frac1n\lfloor\frac nx\rfloor}=\frac{n+\sum_{d\mid x,d\neq x}f[d]\sum_{i=1}^{\lfloor\frac nd\rfloor}[(i,\lfloor\frac xd\rfloor)=1]}{n-\lfloor\frac nx\rfloor}\]

这样就可以啦。答案是\(\sum_{i=1}^n\frac1nf[i]\)
复杂度大概是\(O(n\log^2n)\)的?
还有其它更优的做法,不想看惹(甚至有线性的,瑟瑟发抖)。


//140ms 10900KB
#include <cstdio>
#include <vector>
#include <algorithm>
#define mod 1000000007
typedef long long LL;
const int N=1e5+5;

int P[N>>2],mu[N],f[N],inv[N];
bool notP[N];
std::vector<int> fac[N];

inline int FP(int x,int k)
{
    int t=1;
    for(; k; k>>=1,x=1ll*x*x%mod) k&1&&(t=1ll*t*x%mod);
    return t;
}
inline int G(int n,int k)
{
    LL ans=0;
    for(int i=0,l=fac[k].size(); i<l; ++i) ans+=mu[fac[k][i]]*(n/fac[k][i]);
    return ans%mod;
}
void Init(const int n)
{
    mu[1]=1;
    for(int i=2,cnt=0; i<=n; ++i)
    {
        if(!notP[i]) P[++cnt]=i, mu[i]=-1;
        for(int j=1,v; j<=cnt&&(v=i*P[j])<=n; ++j)
        {
            notP[v]=1;
            if(i%P[j]) mu[v]=-mu[i];
            else break;
        }
    }
    for(int i=1; i<=n; ++i)
        for(int j=i; j<=n; j+=i) fac[j].push_back(i);
}

int main()
{
    int n; scanf("%d",&n);
    Init(n);
    f[1]=1, inv[1]=1;
    for(int i=2; i<=n; ++i) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    for(int x=2; x<=n; ++x)
    {
        LL tmp=n;
        const std::vector<int> &v=fac[x];
        for(int i=0,l=v.size(),d; i+1<l; ++i) d=v[i], tmp+=1ll*f[d]*G(n/d,x/d);
        f[x]=tmp%mod*inv[n-n/x]%mod;
    }
    LL ans=0;
    for(int i=1; i<=n; ++i) ans+=f[i];
    printf("%d\n",int(ans%mod*inv[n]%mod));

    return 0;
}

转载于:https://www.cnblogs.com/SovietPower/p/10692976.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值