CF1097D Makoto and a Blackboard 积性函数、概率期望、DP

传送门

---

比赛秒写完ABC结果不会D……最后C还fst了qwq

首先可以想到一个约数个数(^2)乘上(K)的暴力DP，但是显然会被卡

在(10^{15})范围内因数最多的数是(978217616376000=2^6 imes 3^4 imes 5^3 imes 7^2 imes 11 imes 13 imes 17 imes 19 imes 23 imes 29)，它有(26880)个因数

但是不难发现：在我们的答案中参与计算的只有约数个数函数和约数和函数。它们都是传统的积性函数。这给我们一些启示：可以考虑分解质因数然后DP。

考虑对一个数(x=p^k)进行(DP)，设(f_{i,j})表示初始数字为(p^j)、做(i)轮操作的期望值，转移为(f_{i,j} = frac{sumlimits_{k=0} ^ j f_{i-1,k}}{j+1})，使用前缀和优化转移。最后将所有质因数得到的答案乘起来就是最后的答案。

#include<bits/stdc++.h>
#define int long long
//This code is written by Itst
using namespace std;

inline int read(){
    int a = 0;
    char c = getchar();
    bool f = 0;
    while(!isdigit(c) && c != EOF){
        if(c == '-')
            f = 1;
        c = getchar();
    }
    if(c == EOF)
        exit(0);
    while(isdigit(c)){
        a = a * 10 + c - 48;
        c = getchar();
    }
    return f ? -a : a;
}

const int MOD = 1e9 + 7;
int dp[50][10010] , inv[52];

inline int poww(int a , int b){
    int times = 1;
    while(b){
        if(b & 1)
            times = times * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return times;
}

signed main(){
    for(int i = 1 ; i <= 51 ; ++i)
        inv[i] = poww(i , MOD - 2);
    int N = read() , K = read() , ans = 1;
    for(int i = 2 ; i * i <= N ; ++i)
        if(N % i == 0){
            int cnt = 0;
            while(N % i == 0){
                ++cnt;
                N /= i;
            }
            dp[0][0] = 1;
            int tms = i;
            for(int j = 1 ; j <= cnt ; ++j , tms = tms * i % MOD)
                dp[j][0] = (dp[j - 1][0] + tms) % MOD;
            for(int j = 1 ; j <= K ; ++j){
                dp[0][j] = 1;
                for(int k = 1 ; k <= cnt ; ++k)
                    dp[k][j] = (dp[k][j - 1] * inv[k + 1] + dp[k - 1][j]) % MOD;
            }
            ans = ans * (dp[cnt][K] - dp[cnt - 1][K] + MOD) % MOD;
        }
    if(N != 1)
        ans = ans * ((poww(poww(2 , K) , MOD - 2) * (N % MOD) + MOD + 1 - (poww(poww(2 , K) , MOD - 2))) % MOD) % MOD;
    cout << ans;
    return 0;
}