Codeforces1097D. Makoto and a Blackboard（数论+dp+概率期望）

题目链接：传送门

题目大意：

　　给出一个整数n写在黑板上，每次操作会将黑板上的数（初始值为n）等概率随机替换成它的因子。

　　问k次操作之后，留在黑板上的数的期望。

　　要求结果对10⁹+7取模，若结果不是整数，则用分数表示，并对10⁹+7取逆元。

　　（1 ≤ n ≤ 10¹⁵， 1 ≤ k ≤ 10⁴）

思路：

　　首先我们要知道，在模10⁹+7的范围内，可以任意进行模10⁹+7的加减乘除运算，因为一个给定的数值，它在模10⁹+7条件下的值是唯一确定的。

　　然后我们只要正常地计算，在每次运算之后对10⁹+7取模就好了。

　　假设一个数p^j的出现概率为x，其中p为质数，j为任意非负整数（假设初始的数为p^cc的话，那么0 ≤ j ≤ cc）。那么一次操作留下的数只能是p^t（0 ≤ t ≤ j），因为是等概率分布，所以每个留下的数的概率均为x/j = x * inv(j)。（inv表示在模10⁹+7下取逆元）

　　标记状态：dp[i][j]为第i次操作后p的j次幂出现的概率；

　　那么状态转移方程为：dp[i+1][t] += dp[i][j] * inv(j)；（0 ≤ t ≤ j）

　　而题目给出的数可能不是一个质数的幂，这怎么办呢？我们知道任意一个正整数可以表示为p₁^cc₁*p₂^cc₂*…*p_k^cc_k，对于所有的p_i^cc_i，他们留下的数p_i^j_i的乘积就是最后留下的数，所以他们对答案的贡献的乘积就是最后的答案。

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int md = 1e9 + 7;

inline void add(int& a, int b) {
    a += b;
    if (a > md)
        a -= md;
}

inline void sub(int& a, int b) {
    a -= b;
    if (a < 0)
        a += md;
}

inline int mul(int a, int b) {
    return (int) (1LL * a * b % md);
}

int fpow(int a, int p) {
    int res = 1;
    for (; p; p >>= 1) {
        if (p & 1)
            res = mul(res, a);
        a = mul(a, a);
    }
    return res;
}

int inv(int x) {
    return fpow(x, md-2);
}

int main()
{
    std::ios::sync_with_stdio(false);
    cin.tie(nullptr);
    ll n;
    int k;
    cin >> n >> k;
    vector <pair<ll, int> > f;
    for (ll i = 2; i <= n/i; i++) {
        if (n % i == 0) {
            int cc = 0;
            while (n % i == 0) {
                n /= i;
                cc++;
            }
            f.emplace_back(i, cc);
        }
    }
    if (n > 1)
        f.emplace_back(n, 1);

    int ans = 1;
    for (auto& p : f) {
        int cc = p.second;
        vector <vector<int> > dp(k+1, vector<int>(cc+1, 0));
        dp[0][cc] = 1;
        for (int i = 0; i < k; i++) {
            for (int j = 0; j <= cc; j++) {
                int tmp = mul(dp[i][j], inv(j+1));
                for (int t = 0; t <= j; t++)
                    add(dp[i+1][t], tmp);
            }
        }

        int x = 1, res = 0;
        for (int i = 0; i <= cc; i++) {
            add(res, mul(x, dp[k][i]));
            x = mul(x, (int)(p.first % md));
        }
        ans = mul(ans, res);
    }
    cout << ans << endl;
    return 0;
}