【DP】【CF1097D】 Makoto and a Blackboard
Description
给定一个数 \(n\)，对它进行 \(k\) 次操作，每次将当前的数改为自己的因数，包括 \(1\) 和自己。写出变成所有因数的概率是相等的。求 \(k\) 次以后 \(n\) 期望会变成多少
Input
一行两个整数 \(n,k\)
Output
一行一个整数代表答案
Hint
\(1~\leq~n~\leq~10^{15}~,~1~\leq~k~\leq~10^4\)
Solution
Hello 2019！
我们考虑整个数字变化的树形图：
以 \(n~=~6~,~k~=~2\) 为例：
（恕我直言，这图是真的丑）
然后我们将这张图改一下，要求每个非叶节点都有 \(4\) 个孩子，如果孩子数不足，就让他们平分这四个。于是新的变化图如下：
（恕我直言，这图更丑了）
我们发现这张图上所有叶节点出现的概率都是等可能的，并且他们包含了所有的情况，所以求出这些数字的和再除以节点数就是期望值。
接着考虑由于一个数字的所有质因数都是独立的，同时每个质因子都可以画出类似的转移图，依据唯一分解定理，我们可以将质因子拆开，对于每个质因子求出他们期望变成多少，然后乘起来即为总期望。这里的质因子包括了质数 \(p\) 和它的指数 \(c\)。
接着考虑我们对质因子DP。一次对 \(p^c\) 的操作相当于将其变成 \(p^0,p^1,p^2,\dots p^c\) 中的任意一个。于是我们可以将一次操作转化成将指数 \(c\) 变成 \([1,c]\) 中任意一个数。我们规定任何一个节点都有 \(d(n)\) 个孩子， \(d(n)\) 为 \(n\) 的因数个数。则对于一个数值为 \(x\) 的节点，她有 \(\frac{1}{x + 1}\) 的概率变成 \(y~(0~\leq~y~\leq~x)\)，那么它会占据 \(x\) 的 \(\frac{d(n)}{x + 1}\) 个孩子。
例如，对于 \(2^2\) 进行 \(2\) 次变化的指数图如下
我们认为粉色的位置每个数占据了 \(\frac{1}{2}\) 个儿子
于是我们考虑当我们进行了 \(i\) 次变化，第 \(w\) 个质因子在树形图上的第 \(i\) 层指数数值为 \(j\) 的节点个数，则显然有
其中 \(c_j\) 为 \(n\) 的唯一分解式中对应项的指数。
于是我们枚举 \(n\) 的因数，求出他们出现的次数与他们值的乘积，除以最下面一层的节点个数即为答案
其中 \(h\) 为 \(d\) 唯一分解式对应项的指数。
我们发现在 DP时，第 \(k\) 层每个位置都被乘了 \(d^k(n)\)，于是可以和式子中的 \(\frac{1}{d^k(n)}\) 约掉，式子变为
看起来舒服多了。
考虑复杂度：我们枚举了 \(k\) 层，每层枚举了 \(n\) 的质因数个数次，这里的质因数不包括指数，但是包括重复的，例如 \(p^2\) 算作两个质因数。由于一个数 \(n\) 的质因数为 \(O(\log n)\)，所以总状态复杂度为 \(O(k \log n)\)。在转移时，我们枚举的上界是 \(O(\log n)\)，所以总复杂度为 \(O(k~\log^2 n)\)。发现转移的位置事实上转移了一个加权后缀和，于是我们再开一个数组维护这个后缀和即可做到 \(O(k \log n)\)。
注意到这样会爆空间，但是发现对于每个DP状态我们只需要第 \(k\) 层的情况，并且第一维互不影响，于是我们可以把第一维省掉，每次DP完一个质因子用另一个数组记录第 \(k\) 层的答案，即可做到空间复杂度 \(O(\log^2 n~+~k~\log n)\)。（其实动态开DP数组可以做到空间复杂度 \(O(k~\log n)\)，但是没啥意义）
最后统计答案时直接爆搜因数即可。
Code
#include <cstdio>
#include <cstring>
#include <vector>
#ifdef ONLINE_JUDGE
#define freopen(a, b, c)
#endif
#define rg register
#define ci const int
#define cl const long long
typedef long long ll;
namespace IPT {
const int L = 1000000;
char buf[L], *front=buf, *end=buf;
char GetChar() {
if (front == end) {
end = buf + fread(front = buf, 1, L, stdin);
if (front == end) return -1;
}
return *(front++);
}
}
template <typename T>
inline void qr(T &x) {
rg char ch = IPT::GetChar(), lst = ' ';
while ((ch > '9') || (ch < '0')) lst = ch, ch=IPT::GetChar();
while ((ch >= '0') && (ch <= '9')) x = (x << 1) + (x << 3) + (ch ^ 48), ch = IPT::GetChar();
if (lst == '-') x = -x;
}
template <typename T>
inline void ReadDb(T &x) {
rg char ch = IPT::GetChar(), lst = ' ';
while ((ch > '9') || (ch < '0')) lst = ch, ch = IPT::GetChar();
while ((ch >= '0') && (ch <= '9')) x = x * 10 + (ch ^ 48), ch = IPT::GetChar();
if (ch == '.') {
ch = IPT::GetChar();
double base = 1;
while ((ch >= '0') && (ch <= '9')) x += (ch ^ 48) * ((base *= 0.1)), ch = IPT::GetChar();
}
if (lst == '-') x = -x;
}
namespace OPT {
char buf[120];
}
template <typename T>
inline void qw(T x, const char aft, const bool pt) {
if (x < 0) {x = -x, putchar('-');}
rg int top=0;
do {OPT::buf[++top] = x % 10 + '0';} while (x /= 10);
while (top) putchar(OPT::buf[top--]);
if (pt) putchar(aft);
}
const int maxt = 60;
const int maxm = 10010;
const int MOD = 1000000007;
ll n, kk, ans;
ll frog[maxm][maxt], inv[maxt], gorf[maxt][maxt], sum[maxm][maxt];
std::vector<ll> d;
int cnt;
int c[maxt];
void Get_Inv(ci, ci);
void dfs(ci, ci);
int mpow(cl x, int);
signed main() {
freopen("1.in", "r", stdin);
qr(n); qr(kk);
Get_Inv(59, MOD);
ll dn = n;d.push_back(0);
for (ll i = 2; (i * i) <= n; ++i) if(!(dn % i)) {
d.push_back(i); ++cnt;
while (!(dn % i)) {dn /= i; ++c[cnt];}
}
if (dn != 1) {++cnt; d.push_back(dn); c[cnt] = 1;}
for (rg int j = 1; j <= cnt; ++j) {
frog[0][c[j]] = sum[0][c[j]] = inv[c[j] + 1];
for (rg int i = 0; i < c[j]; ++i) frog[0][i] = 0, sum[0][i] = sum[0][c[j]];
for (rg int i = 1; i <= kk; ++i) {
rg int di = i - 1;
sum[i][c[j] + 1] = 0;
for (rg int k = c[j]; ~k; --k) {
frog[i][k] = sum[di][k];
sum[i][k] = (sum[i][k + 1] + frog[i][k] * inv[k + 1] % MOD) % MOD;
}
}
for (rg int k = 0; k <= c[j]; ++k) gorf[j][k] = frog[kk][k];
}
dfs(1, 1);
qw((ans + MOD) % MOD, '\n', true);
return 0;
}
void Get_Inv(ci x, ci p) {
inv[1] = 1;
for (rg int i = 2; i <= x; ++i) inv[i] = - p / i * inv[p % i] % MOD;
}
void dfs(ci cur, ci v) {
if (cur > cnt) {ans = (ans + v) % MOD; return;}
for (int i = 0; i <= c[cur]; ++i) dfs(cur + 1, 1ll * gorf[cur][i] * v % MOD * mpow(d[cur], i) % MOD);
}
int mpow(cl x, int y) {
ll _ret = 1, _temp = x % MOD;
while (y) {
if (y & 1) _ret = _ret * _temp % MOD;
_temp = _temp * _temp % MOD;
y >>= 1;
}
return _ret;
}
