CF1139D Steps to One 题解莫比乌斯反演枚举DP

Posted 2021-02-14 wjyyy

反演套 DP 的好题（不用反演貌似也能做

Description

Vivek initially has an empty array (a) and some integer constant (m).

He performs the following algorithm:

1. Select a random integer (x) uniformly in range from (1) to (m) and append it to the end of (a).
2. Compute the greatest common divisor of integers in (a).
3. In case it equals to (1), break
4. Otherwise, return to step (1).

Find the expected length of (a). It can be shown that it can be represented as (frac PQ) where (P) and (Q) are coprime integers and (Q eq 0pmod{10^9+7}). Print the value of (Pcdot Q^{-1}pmod{10^9+7}).

Input

The first and only line contains a single integer (m)((1le mle 100000)).

Output

Print a single integer — the expected length of the array (a) written as (Pcdot Q^{-1}pmod{10^9+7}).

Examples

input

1

output

1

input

2

output

2

input

4

output

333333338

Note

In the first example, since Vivek can choose only integers from (1) to (1), he will have (a=[1]) after the first append operation, and after that quit the algorithm. Hence the length of (a) is always (1), so its expected value is (1) as well.

In the second example, Vivek each time will append either (1) or (2), so after finishing the algorithm he will end up having some number of (2)‘s (possibly zero), and a single (1) in the end. The expected length of the list is (1?frac 12+2?frac 1{2^2}+3?frac 1{2^3}+dots =2).

题意

每一步在 (1sim m) 中任选一个整数，问期望多少步后选出的数的最大公约数是 (1)。答案对 (1 000 000 007) 取模。

题解

因为每一步不会让已经选了的元素的 (gcd) 和变大，因此认为是一个除自环外的有向无环图。对于自环我们很好处理，所以把它看成是一道期望 DP。

令 (f[i]) 表示当前的 (gcd) 和为 (i)，到 (gcd) 和为 (1) 的状态的期望步数。因此把状态转移方程写出来
[ f[i]=1+frac{sum_{j=1}^m f[gcd(i,j)]}m ]
这样的转移是 (O(m^2)) 的。但是我们发现，对于很多 (j)，(gcd(i,j)) 都是相等的，因此我们把这样的数整合到一起。

令 (F(n)) 表示 (1sim m) 中有多少个数 (i) 满足 (gcd(x,i)=n)，其中视 (x?) 为常数。

则计算 (f[i]) 就转化为了
[ f[i]=1+frac{sum_{d|i}f[d] imes F(d)}n,x=i ]
这样差不多就把枚举优化到了 (log n) 的 (d|i?)。

考虑怎么计算 (F(n))
[ F(n)=sum_{i=1}^m[gcd(x,i)=n] ]
令 (G(n)=sum_{n|d}F(d))，则
[ egin{aligned} G(n)&=sum_{n|d}F(d)&=sum_{n|d}sum_{i=1}^m[gcd(x,i)=d]&=sum_{i=1}^m[n|gcd(x,i)]?&=sum_{i=1}^{leftlfloorfrac mn ight floor}left[1|gcdleft(frac xn,i ight) ight] end{aligned} ]
实际上这样是有问题的，因为（在后面）无法保证 (n|x)，此时 (G(n)) 就一定为 (0) 了。

我们再多化一步：
[ egin{aligned} G(n)&=sum_{i=1}^{leftlfloorfrac mn ight floor}left[1|gcdleft(frac xn,i ight) ight][n|x]&=leftlfloorfrac mn ight floorcdot[n|x] end{aligned} ]
根据 (G(n)=sum_{n|d}F(d))，我们反演到 (F)，得
[ egin{aligned} F(n)&=sum_{n|d}muleft(frac dn ight)G(d)&=sum_{i=1}^{leftlfloorfrac mn ight floor}mu(i)G(ni)&=sum_{i=1}^{leftlfloorfrac mn ight floor}mu(i)leftlfloorfrac{m}{ni} ight floor[(ni)|x] end{aligned} ]
我们发现后面的布尔表达式可以当作条件。原本的条件本来就是 ( o +infty) 的，只不过超过了 (leftlfloorfrac mn ight floor) 没有意义。因此直接把条件换成 ([(ni)|x]) 即可。又因为 (n|x) 在上面的枚举过程中是成立的，同时可以转化为 (left[i|frac xn ight])。
[ F(n)=sum_{i|frac xn}mu(i)leftlfloorfrac{m}{ni} ight floor ]
这样的一次枚举是 (Oleft(dleft(frac xn ight) ight)) 的，由于 (1sim m) 的约数个数和均摊是 (O(log m)) 的，其中最多的有 (128) 个约数，但是这样的数肯定不是很多，并且其中很多被枚举到的数都是质因数，迭代一下并不会造成很大的复杂度。

然后我们需要再把状态转移方程稍微转化一下，把 (f[i]) 移到左边
[ egin{aligned} f[i]&=1+frac{sum_{d|i,d<i}f[d] imes F(d)+f[i] imes F(i)}n,x=i\frac{n-F(i)}{n}cdot f[i]&=1+frac{sum_{d|i,d<i}f[d] imes F(d)}n,x=if[i]&=frac{n+sum_{d|i,d<i}f[d] imes F(d)}{n-F(i)},x=i end{aligned} ]
就得到了真正的转移方程。

时间复杂度 (O(mlog^2 m))。

Code：

#include<cstdio>
#include<cstring>
#include<vector>
#define ll long long
#define p 1000000007
using std::vector;
vector<int> v[100100];//约数用 vector 存一下，每次 √m 枚举不是很稳
ll qpow(ll x,ll y)
{
    ll ans=1;
    while(y)
    {
        if(y&1)
            ans=ans*x%p;
        x=x*x%p;
        y>>=1;
    }
    return ans;
}
bool is[100100];
int pri[100100],mu[100100],cnt=0;
ll f[100100];
int n;
int calc(int x,int y)//1~n 中 gcd(x,i)=y 的数的个数
{
    int g=x/y,ans=0;
    for(int i=0;i<v[g].size();++i)
        ans+=mu[v[g][i]]*(n/v[g][i]/y);
    return ans;
}
int main()
{
    scanf("%d",&n);
    f[1]=1;
    mu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!is[i])
        {
            pri[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=n;++j)
        {
            is[i*pri[j]]=1;
            if(i%pri[j])
                mu[i*pri[j]]=-mu[i];
            else
            {
                mu[i*pri[j]]=0;
                break;
            }
        }
    }
    for(int i=1;i<=n;++i)
        for(int j=i;j<=n;j+=i)
            v[j].push_back(i);
    ll ans=1,inv=qpow(n,p-2);
    for(int i=2;i<=n;++i)
    {
        for(int j=0;j<v[i].size()-1;++j)
            f[i]=(f[i]+calc(i,v[i][j])*f[v[i][j]]%p)%p;
        f[i]=(f[i]*inv+1)%p;
        ll g=n-calc(i,i);
        f[i]=f[i]*n%p*qpow(g,p-2)%p;
        ans=(ans+f[i])%p;
    }
    printf("%lld
",ans*qpow(n,p-2)%p);
    return 0;
}