题目连接:https://codeforces.com/contest/1097/problem/D
首先,对每一类质数分开讨论,现在问题变成了,设\(n\)中质因子\(p\)的次数是\(i\),对于一个数\(x=p^j\),经过\(k\)轮之后它被操作出来的概率。显然这跟\(p\)是多少无关,所以记上述情况为\(F_{ij}\)。\(F_{ij}\)并不好直接计算,所以考虑令\(f_{ilj}\)表示\(n\)中质因子\(p\)的次数是\(i\),对于一个数\(x=p^j\),经过\(l\)轮之后它被操作出来的概率。\(f_{ilj}\)有如下转移
\[f_{ili}=1 (l=0)\]
\[f_{ilj}=\sum_{k=j}^{i}f_{il-1k} \times \frac{1}{k+1} (l>1)\]
第二个转移式子是一段连续区间的和显然可以前缀和优化。注意到\(l\)是由\(l-1\)转移过来的,故空间上也可以用滚动数组优化。接下来深搜枚举\(n\)的约数分解后的各个质数的次数并计算答案即可。
时间复杂度\(O(k\log_{2}^{2}n+\sqrt n)\)。
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <iostream>
#include <algorithm>
#include <stack>
#include <queue>
#include <map>
#include <set>
#include <iomanip>
#include <assert.h>
#include <fstream>
using namespace std;
typedef long long ll;
const int MAXP = 55;
const int MAXK = 10005;
const ll MOD = 1000000007;
int k,tot;
int cnt[MAXP];
ll n,ans;
ll p[MAXP];
ll inv[MAXP];
ll f[MAXP][MAXP];
ll g[MAXP][MAXP];
ll h[MAXP][MAXP];
ll power(ll a,ll b)
{
ll res = 1;
while (b)
{
if (b & 1)
res = res * a % MOD;
a = a * a % MOD;
b >>= 1;
}
return res;
}
void init()
{
for (int i = 1;i <= 50;i++)
inv[i] = power(i,MOD - 2);
for (int i = 0;i <= 50;i++)
f[i][i] = 1;
for (int i = 0;i <= 50;i++)
{
h[i][0] = f[i][0] * inv[1] % MOD;
for (int j = 1;j <= i;j++)
h[i][j] = (h[i][j - 1] + f[i][j] * inv[j + 1]) % MOD;
}
for (int i = 1;i <= k;i++)
{
for (int j = 0;j <= 50;j++)
for (int k = 0;k <= j;k++)
g[j][k] = (h[j][j] - (k > 0 ? h[j][k - 1] : 0)) % MOD;
memcpy(f,g,sizeof(f));
memset(g,0,sizeof(g));
for (int j = 0;j <= 50;j++)
{
h[j][0] = f[j][0] * inv[1] % MOD;
for (int k = 1;k <= j;k++)
h[j][k] = (h[j][k - 1] + f[j][k] * inv[k + 1]) % MOD;
}
}
}
void dfs(int now,ll N,ll P)
{
if (now > tot)
{
(ans += N % MOD * P) %= MOD;
return;
}
for (int i = 0;i <= cnt[now];i++)
{
dfs(now + 1,N,P * f[cnt[now]][i] % MOD);
N *= p[now];
}
}
int main()
{
cin >> n >> k;
init();
ll N = n;
for (ll i = 2;i * i <= N;i++)
if (N % i == 0)
{
cnt[++tot] = 0;
p[tot] = i;
while (N % i == 0)
{
cnt[tot]++;
N /= i;
}
}
if (N != 1)
{
cnt[++tot] = 1;
p[tot] = N;
}
dfs(1,1,1);
cout << (ans + MOD) % MOD << endl;
return 0;
}