算是記一下昨天晚上都想了些什麼

官方題解   點我

簡單題意

給定兩個正整數$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}f(i, 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\n", ans);
    return 0;
}