HELLO 2019 - Makoto and a Blackboard （积性函数+DP期望）

Makoto and a Blackboard

题目链接： D - Makoto and a Blackboard

题意

给你一个N，定义一个操作：将N替换为他的一个因子（包括1和N）

现在重复K次以上操作，问最后期望的值是多少？

数据范围： $N<=10^{15},K<=10^4$

---

思路

我的思路

求期望的有两种做法，

1. DP存的是概率，直接是最后所有可能的结果乘以他们的概率，加起来
2. DP存的是期望，利用公式 $总期望=所有可以转移的子期望\ast 他们发生的概率$

但其实这两种期望并不关键，最关键的发现积性函数这个性质

官方思路

How to solve the problem if $n$ is a prime power $p^{\alpha }$? We can simply use dynamic programming. Let $D{P}_{i,j}$ denote the probability that after ii steps the current number on the blackboard is $p^j$. At the beginning, $D{P}_{0,\alpha }=1$. The transitions in DP are rather straightforward. The expected value of the number on the blackboard is then
$$\sum _{j=0}^{\alpha }D{P}_{k,j}{p}^j$$

What if $n$ is not a prime power? We can observe that the result is multiplicative on $n$. Therefore, we factorize $n=p_1^{\alpha 1}{p}_2^{\alpha 2}\dots p_{\ell }^{\alpha \ell }$, apply the solution above for each prime power $p_1^{\alpha 1}$ through $p_{\ell }^{\alpha \ell }$, and eventually multiply the results.

Another problem is that we need to write a single integer modulo $10^9+7$ even though the answer is a rational number. To cope with this, we notice that $10^9+7$ is prime. Therefore, thanks to the Fermat Little Theorem, dividing by some number xx is equivalent to multiplying by $x^{10^9+5}$, modulo $10^9+7$ . We use the quick exponentiation to compute the results.

The overall complexity is $ O(\sqrt{n}+k \cdot \log n)$

这里非常重要的就是积性函数，既 $f(6)=f(2)\ast f(3),f(30)=f(5)\ast f(2)\ast f(3)$

关于积性函数链接：https://www.cnblogs.com/zhoushuyu/p/8275530.html

---

代码

#include<bits/stdc++.h>
using namespace std;
#define rep(i,j,k) for(ll i = (ll)j;i <= (ll)k;i ++)
#define debug(x) cerr<<#x<<":"<<x<<endl
#define pb push_back

typedef long long ll;
const ll MAXN = 1e6+7;
const ll MOD = 1e9+7;

ll dp[10007][107]; //dp[i][j] : the possibility of after i times the j power now value

ll inv[107],n,k;

ll fpow(ll a,ll x) {
    ll res = 1;
    while (x) {
        if (x&1) res = res*a%MOD;
        a = a*a%MOD;
        x >>= 1;
    }
    return res;
}

ll solve (ll x,ll ci) {
    rep(i,0,k) rep(j,0,ci) dp[i][j] = 0;
    dp[0][ci] = 1;
    rep(i,1,k) {
        rep(j,0,ci) {
            rep(m,j,ci) {
                dp[i][j] = (dp[i][j] + dp[i-1][m]*(inv[m+1])%MOD)%MOD;
            }
        }
    }

    ll res = 0;
    rep(i,0,ci) {
        res = (res + fpow(x,i)*dp[k][i]%MOD)%MOD;
    }

    return res;
}

int main()
{
    ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);

    rep(i,1,100) inv[i] = fpow(i*1LL,MOD-2);

    cin >> n >> k;

    ll ans = 1;
    for (ll i = 2;i*i <= n;i ++) {
        if (n%i==0) {
            ll cnt = 0;
            while (n%i==0) n/=i,cnt++;
            ans = ans*solve(i,cnt)%MOD;
        }
    }
    if (n!=1) ans = ans*solve(n,1)%MOD;

    cout << ans << endl;
}