CF1139 D. Steps to One(期望dp，数论)

最新推荐文章于 2021-04-14 17:28:10 发布

原创最新推荐文章于 2021-04-14 17:28:10 发布 · 539 阅读

2 ·

CC 4.0 BY-SA版权

dp 同时被 2 个专栏收录

9 篇文章

订阅专栏

数论

4 篇文章

订阅专栏

该博客讨论了CF1139 D题，题目涉及在给定范围内的随机整数操作，计算序列的最小公倍数，直到达到某个条件。博主分析了问题并提出用期望dp的方法来解决，通过建立数学模型计算到达特定状态的期望循环次数，最终给出O(n^2)复杂度的解决方案。

题目链接：D. Steps to One

题意：

给你一个整数m(1≤m≤100000)，刚开始序列a为空，循环进行以下操作：

随机在[1,m]选择一个整数x，将x加入序列a尾部
计算集合a的最小公倍数
如果 $GCD(a)=i$ ，结束。
否则，跳到1.

问循环结束时序列a的长度期望，对1e9+7取模。

题目分析：

设 $F(i)$ 表示 $GCD(a)=i$ 时，到达 $GCD(a)=1$ 时所需的期望循环次数

那么所求答案可以表示为： $Ans=1+\sum_{i=1}^{m} \frac{F(i)}{m}$

1代表第一次取出一个数字，求和是选到第一个数字之后的期望次数*概率

假设当前 $GCD(a)=n$ ，下一个数字为 $i$ ，那么 $GCD(a)=gcd(n,i)$ ，我们可以得到：

$F(n)=1+\sum_{i=1}^{m}\frac{F(gcd(n,i))}{m}$

将式子左右同乘 $m$ ，设 $d = gcd(n,i)$ ，则有：

$mF(n)=m+\sum_{d|n}\left \[ F(d)*\sum_{i=1}^{m}[gcd(n,i)=d]\right \]$

对内部的和式进行化简：

$\large \sum_{i=1}^{m}[gcd(n,i)=d]\\=\sum_{i=1}^{[m/d]}[gcd(\frac{n}{d},i)=1]\\=\sum_{i=1}^{[m/d]}\sum_{t|\frac{n}{d}\&t|i}^{ }\mu(t)\\=[\frac{[m/d]}{t}]\sum_{t|\frac{n}{d}}^{ }\mu(t)$

代回原式：

$mF(n) \\=m+\sum_{d|n}\left \[ F(d)*\sum_{t|(n/d)}\mu(t)*[[m/d]/t]\right \] \\=m+\sum_{d|n\&d\neq n}\left \[ F(d)*\sum_{t|(n/d)}\mu(t)*[[m/d]/t]\right \]+F(n)*[\frac{m}{n}]$

将F(n)移到同一边得到：

$F(n) = \left ( m+\sum_{d|n\&d\neq n}\left \[ F(d)*\sum_{t|(n/d)}\mu(t)*[[m/d]/t]\right \] \right )/\left ( m-[\frac{m}{n}] \right )$

计算这个式子需要枚举 $n$ 的因子 $d$ ，对于每个 $n/d$ 再次枚举因子 $t$ ，复杂度玄学，O(能过)

代码：

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll INF = 0x3f3f3f3f3f3f3f3f;
const ll mod = 1e9+7;
const int MAXN = 100005;
int m;
ll F[MAXN], inv[MAXN];

int mu[MAXN], p[MAXN], tot;
bool nprime[MAXN];

struct node{
    int fac, next;
}e[MAXN*100];
int head[MAXN], cnt;
inline void inser(int k, int d){//d|k
    e[cnt] = (node){d,head[k]}, head[k] = cnt++;
}

inline void add(int &a, ll b){
    a += b; if(a >= mod) a -= mod;
}
void init(){
    tot = cnt = 0;
    //mu函数
    mu[1] = 1;  nprime[0] = nprime[1] = 1;
    for(int i = 2;i<MAXN;i++){
        if(!nprime[i]) mu[i] = -1, p[tot++] = i;
        for(int j = 0;j<tot && i*p[j] < MAXN;j++){
            nprime[i*p[j]] = 1;
            if(i%p[j]==0){
                mu[i*p[j]] = 0;
                break;
            }
            mu[i*p[j]] = -mu[i];
        }
    }
    //处理因子d|n && d != n
    memset(head, -1, sizeof head);
    for(int i = 2;i<(MAXN>>1);i++)
        for(int j = (i<<1);j<MAXN;j+=i)
            inser(j,i);
    //逆元
	inv[1] = 1;
	for (ll i = 2; i < MAXN; ++i)
        inv[i] = 1ll*(mod-mod/i)*inv[mod%i]%mod;
    //F[i]表示当前gcd为i，到达gcd为1所需的期望长度
    F[1] = 0;
    for(int i = 2;i <= m;i++){
        int tmp = 0;
        for(int j = head[i];~j;j = e[j].next){
            int d = e[j].fac;
            int x = m/d;
            int tmpp = 1ll*mu[i/d]*(x/(i/d))%mod;
            for(int k = head[i/d];~k;k = e[k].next)
                add(tmpp, 1ll*mu[e[k].fac]*(x/e[k].fac)%mod);
            add(tmpp, x);
            add(tmp, 1ll*tmpp*F[d]%mod);
        }
        add(tmp, m);
        F[i] = 1ll*tmp*inv[m-(m/i)]%mod;
    }
}
int main() {
    scanf("%d",&m);
    init();
    int ans = 0;
    for(int i = 1;i<=m;i++)add(ans,F[i]);
    ans = 1ll*ans*inv[m]%mod;
    add(ans, 1);
    printf("%d\n",ans);
    return 0;
}