Codeforce 1139D Step to one（概率dp+容斥原理）

#include<iostream>
#include<cmath>
#include<ctime>
#include<cstdlib>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N=1e5+10;
const ll mod=1e9+7;
ll dp[N];
int fac[20];
int p=1;
ll f(int n,int x)
{
    p=1;
    ll ans=n;
    int tp=x;
    for(int i=2;i*i<=x;i++)
    {
        if(tp%i==0)
        {
            fac[p]=i;
            p++;
            while(tp%i==0)
            {
                tp=tp/i;
            }
        }
        if(tp==1)break;
    }
    if(tp!=1)
    {
        fac[p]=tp;
        p++;
    }
    for(int i=1;i<(1<<(p-1));i++)
    {
        int cnt=0,pow=1;
        for(int j=1;j<p;j++)
        {
            if((i&(1<<(j-1))))
            {
                cnt++;
                pow=pow*fac[j];
            }
        }
        if(cnt%2==0)
        {
            ans+=n/pow;
        }
        else
        {
            ans-=n/pow;
        }
    }
    return ans;
}
ll powx(ll a,ll b)
{
    if(b==1)return a;
    else
    {
        ll tp=powx(a,b/2);
        if(b%2==0)
        {
            return (tp*tp)%mod;
        }
        else return ((tp*tp)%mod*a)%mod;
    }
}
int gcd(int a,int b)
{
    if(!b)return a;
    else return gcd(b,a%b);
}
int main ()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    int m;
    cin>>m;
    ll ans=0;
    for(int i=2;i<=m;i++)
    {
        dp[i]=m;
        int cnt=m-m/i;
        for(int j=2;j*j<=i;j++)
        {
            if(i%j==0)
            {
                dp[i]=(dp[i]+dp[j]*f(m/j,i/j))%mod;
                if(j*j!=i)
                dp[i]=(dp[i]+dp[i/j]*f(m/(i/j),j))%mod;
            }
        }
        //cout<<dp[i]<<endl;
        ll x=powx(cnt,mod-2);
        dp[i]=(dp[i]*x)%mod;
        ans+=dp[i];
    }
    ans=(((m+ans)%mod)*powx(int(m),int(mod-2)))%mod;
    cout<<ans<<endl;
    return 0;
}