题目地址：洛谷CF1097D

首先考虑 \(n=p^c\) （ \(p\) 为质数）的情况，显然DP：

令 \(f_{i,j}\) 为第 \(i\) 次替换后出现 \(p^j\) 的概率

边界：

\[f_{0,c}=1\]

状态转移方程：

\[f_{i,j}=\sum_{t=j}^{c} \frac{f_{i-1,t}}{t+1}\]

目标：

\[\sum_{j=0}^{c} (f_{k,j}p^j)\]

考虑一般情况，将 \(n\) 分解质因数：

\[n=\prod_{i=1}^{m} {p_i}^{c_i}\]

按照上述方法DP每个 \({p_i}^{c_i}\)

由于期望是积性函数，直接将所有答案乘起来即可 （我就是卡在这一步上，难受QWQ）

代码：

#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int K = 10006, C = 56, P = 1000000007;
ll n, f[K][C], inv[C];
int k;
vector<pair<ll, int> > d;

void divide(ll n) {
    for (ll i = 2; i <= sqrt(n); i++)
        if (n % i == 0) {
            int c = 0;
            while (n % i == 0) {
                n /= i;
                ++c;
            }
            d.push_back(make_pair(i, c));
        }
    if (n > 1ll) d.push_back(make_pair(n, 1));
}

ll work(ll p, int c) {
    for (int i = 0; i <= k; i++)
        for (int j = 0; j <= c; j++)
            f[i][j] = 0;
    f[0][c] = 1;
    for (int i = 1; i <= k; i++)
        for (int j = c; j >= 0; j--)
            for (int t = j; t <= c; t++)
                f[i][j] = (f[i][j] + f[i-1][t] * inv[t+1] % P) % P;
    ll ans = 0, now = 1;
    for (int j = 0; j <= c; j++) {
        ans = (ans + f[k][j] * now % P) % P;
        now = now * p % P;
    }
    return ans;
}

int main() {
    inv[1] = 1;
    for (int i = 2; i < C; i++)
        inv[i] = -(P / i) * inv[P%i] % P;
    cin >> n >> k;
    divide(n);
    ll ans = 1;
    for (unsigned int i = 0; i < d.size(); i++)
        ans = ans * work(d[i].first, d[i].second) % P;
    cout << (ans + P) % P << endl;
    return 0;
}