CF 1097D Makoto and a Blackboard

算是记一下昨天晚上都想了些什么

官方题解   点我

简单题意

给定两个正整数$n$和$k$，定义一步操作为把当前的数字$n$等概率地变成$n$的任何一个约数，求$k$步操作后的期望数字，模$1e9 + 7$。

$$n leq 10^{15}, k leq 10^4$$

我的思路

设$f(n, k)$表示$n$在$k$步操作之后的期望数字，假设$n$的约数有$m$个，分别为$d_1, d_2, dots, d_m$，有递推式

$$f(n, k) = frac{1}{m}sum_{i = 1}^{m}f(d_i, k - 1)$$

边界条件显然是$f(n, 0) = n$。

直接暴算的话一共有$n * k$个状态，无法承受。

接下来证明：$f(a, k) * f(b, k) = f(ab, k)$，（其中$a, b$互质）。

数学归纳法来了（逃）

1、$k = 0$的时候显然成立。

2、假设在$k - 1$的时候成立。

我们设$a$有$n$个约数，$b$有$m$个约数，因为约数个数$sigma$是一个积性函数，所以$ab$的约数个数有$nm$个。

那么根据递推式，有

$$f(ab, k) = frac{1}{nm}sum_{d | ab}f(d, k - 1)$$

$$f(a, k) * f(b, k) = frac{1}{n}sum_{i | a}f(i, k - 1)frac{1}{m}sum_{j | b}f(j, k - 1) = frac{1}{nm}sum_{i | a}sum_{j | b}f(i, k - 1) * f(j, k - 1)$$

要证$f(ab, k) = f(a, k) * f(b, k)$，

即证$sum_{d | ab}f(d, k - 1) = sum_{i | a}sum_{j | b}f(i, k - 1) * f(j, k - 1)$，

右边的式子变形一下

$$sum_{i | ab}sum_{j | i & j | a}f(j, k - 1) * f(frac{i}{j}, k - 1)[gcd(j, frac{i}{j} == 1)]$$

把$ab$分解质因数变成$prod_{i = 1}^{m}p_i^{c_i}$的形式。

注意到$a$、$b$互质，那么$gcd(j, frac{i}{j}) == 1$的时候其实只有一种，那就是$j$恰好取完了某个或某几个$p_i^{c_i}$的时候，这时候有$f(d, k - 1) = f(i, k - 1) * f(frac{d}{i}, k - 1)[gcd(d, a) == i]$。

代进去之后发现两式相等了。

这样子的话我们就得到了$f$函数一个类似于积性的性质，于是我们可以直接把$n$分解成$prod_{i = 1}^{m}p_i^{c_i}$的形式，然后分别计算每一个$p_i^{c_i}$的答案最后乘起来。

发现这样子状态数十分有限，只有$klogn$个，所以直接暴力算就可以了。

时间复杂度应当是$O(sqrt{n} + klogn)$。

代码非常乱。

Code：

#include <cstdio>
#include <cstring>
#include <map>
#include <vector>
#include <algorithm>
#define rep(i, a, b) for (int i = (a); i <= (b); i++)
#define per(i, a, b) for (int i = (a); i >= (b); i--)
using namespace std;
typedef long long ll;
typedef pair <ll, int> pin;

const ll P = 1e9 + 7;

ll ans = 1LL, g[80][10005], inv[80];
bool vis[80][10005];

map <pin, ll> mp;

template <typename T>
inline void inc(T &x, T y) {
    x += y;
    if (x >= P) x -= P; 
}

inline ll fpow(ll x, ll y) {
    ll res = 1LL;
    for (; y > 0; y >>= 1) {
        if (y & 1) res = res * x % P;
        x = x * x % P;
    }
    return res;
}

ll f(ll p, int m, int k) {
    if (p == 1) return 1LL;
    if (k == 0) return fpow(p, m);
    if (vis[m][k]) return g[m][k];
    
    vis[m][k] = 1;
    
    ll res = 0;
    rep(i, 0, m) inc(res, f(p, i, k - 1) * inv[m + 1] % P);
    
    return g[m][k] = res;
}

inline void solve(ll p, int m, int k) {
    rep(i, 0, m) rep(j, 0, k) g[i][j] = 0, vis[i][j] = 0; 
    ans = ans * f(p, m, k) % P;
}

int main() {
//    freopen("Sample.txt", "r", stdin);
//    freopen("out.txt", "w", stdout);

    rep(i, 0, 79) inv[i] = fpow(i, P - 2);
    
    ll n; int k;
    scanf("%I64d%d", &n, &k);
    
    ll tmp = n;
    for (ll i = 2; i * i <= n; i++) {
        if (tmp % i == 0) {
            int m = 0;
            for (; tmp % i == 0; tmp /= i, ++m);
            solve(i, m, k);
        }
    }    
    if (tmp > 1) solve(tmp, 1, k);
    
    printf("%I64d
", ans);
    return 0;
}