Codeforces #548 (Div2)

Problem   Codeforces #548 (Div2) - D.Steps to One

Time Limit: 2000 mSec

Problem Description

Input

 The first and only line contains a single integer mm (1≤m≤100000，1≤m≤100000).

Output

Print a single integer — the expected length of the array aa written as P⋅Q^−1(mod10^9+7)

Sample Input

 4

Sample Output

333333338

题解：概率dp做的太少了，不会做。。。

　　首先是状态定义，dp[x]表示当前序列的gcd为x的情况下，还能添加数字个数的期望值，看了题解之后感觉这样很自然，但是自己想就想不到，有了这个定义之后状态转移就比较简单了，枚举下一个数字，根据期望的线性性质，有：

　　

这样一来得到了一个O(n^2)的算法，显然不行，不过到这里的转化就很套路了，枚举gcd，式子化为：

　　

f(y, x)指的是和x的gcd是y的数有多少个（从1到m），预处理出1~m的约数，复杂度mlogm，这样不用根号m枚举约数，枚举y这一步题解中给出的均摊复杂度是logm，（不会证qwq，不过很多地方都是这么分析的），这样一来就只剩下f的求解了，设x = y * a，则和x的gcd为b的数必定可表达为 p = y * b且gcd(a, b) = 1，这样一来就相当于计算从1到m/y中和a互质的数的个数，也就是和a没有相同的质因子，打打表（看题解）可以发现在题目的范围内b的质因子数至多为6，容斥一下即可计数，最终复杂度O(mlogm*2^6*6)，是可以接受的。代码是看了题解的代码之后写的，主要学习这里质因数分解的操作。

#include <bits/stdc++.h>

using namespace std;

#define REP(i, n) for (int i = 1; i <= (n); i++)
#define sqr(x) ((x) * (x))

const int maxn = 100000 + 100;
const int maxm = 200000 + 100;
const int maxs = 10000 + 10;

typedef long long LL;
typedef pair<int, int> pii;
typedef pair<double, double> pdd;

const LL unit = 1LL;
const int INF = 0x3f3f3f3f;
const LL mod = 1000000007;
const double eps = 1e-14;
const double inf = 1e15;
const double pi = acos(-1.0);

LL pow_mod(LL x, LL n, LL mod)
{
    LL base = x % mod;
    LL ans = 1;
    while (n)
    {
        if (n & 1)
        {
            ans = ans * base % mod;
        }
        base = base * base % mod;
        n >>= 1;
    }
    return ans % mod;
}

LL m, invm;
unordered_map<int, int> prime[maxn];
vector<LL> fact[maxn];
bool is_prime[maxn];
LL dp[maxn];

void premanagement()
{
    memset(is_prime, true, sizeof(is_prime));
    is_prime[0] = is_prime[1] = false;
    for (LL i = 2; i < maxn; i++)
    {
        if (is_prime[i])
        {
            for (LL j = i; j < maxn; j += i)
            {
                LL cnt = 0;
                LL tmp = j;
                while (tmp % i == 0)
                {
                    cnt++;
                    tmp /= i;
                }
                prime[j][i] = cnt;
                is_prime[j] = false;
            }
            is_prime[i] = true;
        }
    }

    for (LL i = 1; i < maxn; i++)
    {
        for (LL j = 2 * i; j < maxn; j += i)
        {
            fact[j].push_back(i);
        }
    }
}

LL cal(LL x, LL n)
{
    vector<LL> a;
    LL curcnt = m / x;
    for (auto &item : prime[n])
    {
        if (!prime[x].count(item.first))
        {
            a.push_back(item.first);
            continue;
        }
        if (prime[x][item.first] == item.second)
        {
            continue;
        }
        else
        {
            a.push_back(item.first);
        }
    }

    LL sz = a.size();
    LL lim = m / x;
    for (LL sit = 1; sit < (unit << sz); sit++)
    {
        LL tag = 1;
        LL val = 1;
        for (LL j = 0; j < sz; j++)
        {
            if (sit & (unit << j))
            {
                val *= a[j];
                tag *= -1;
            }
        }
        curcnt += tag * (lim / val);
    }
    return curcnt;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    //freopen("input.txt", "r", stdin);
    //freopen("output.txt", "w", stdout);
    cin >> m;
    invm = pow_mod(m, mod - 2, mod);
    premanagement();
    dp[1] = 0;
    for (LL i = 2; i <= m; i++)
    {
        LL &ans = dp[i];
        ans = 1;
        for (LL item : fact[i])
        {
            ans += cal(item, i) * dp[item] % mod * invm % mod;
            ans %= mod;
        }
        LL cnt = m / i;
        ans = ans * m % mod * pow_mod(m - cnt, mod - 2, mod) % mod;
    }
    LL ans = 1;
    for (LL i = 1; i <= m; i++)
    {
        ans = (ans + dp[i] * invm) % mod;
    }
    cout << ans << endl;
    return 0;
}