Codeforces 1097D Makoto and a Blackboard

题目链接:https://codeforces.com/problemset/problem/1097/D

思路: 显然先将n素因子分解,  n = (a1^b1)*(a2^b2)...(an^bn)   考虑所有结果中只含有a1^1的项, 和只含有a2^2的项，他们除了a1这个素因子其他的都是一样的， 所以ans((a1^b1)*(a2^b2)...(an^bn)) = ans(a1^b1)*ans((a2^b2)...(an^bn)), 就像积性函数. 现在就是怎样处理ans(a1^b1).  ans(a1^b1) = a1*P(a1)+(a1^2)*P(a1^2)+......+(a1^b1)*P(a1^b1).  对于P(a1^1)， P(a1^21)， P(a1^b1)可以用dp推出来.  dp[i][j] 表示 第i次后得到a1^j的概率.

代码:

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod = 1e9+7;
const int maxn = 65;
const int maxx = 1e4+10;

ll ans = 1;
ll n, k;
ll dp[maxx][maxn];
ll inv[maxn];

ll pow_mod(ll p, ll k) {
    ll ans = 1;
    while(k) {
        if(k & 1)  ans = ans*p%mod;
        p = p*p%mod; k >>= 1;
    }
    return ans;
}

void init() {
    for(ll i = 1; i < maxn; i++)
        inv[i] = pow_mod(i, mod-2);
}

void solve(ll p, ll time) {
    memset(dp, 0, sizeof dp);
    dp[0][time] = 1;
    for(int i = 1; i <= k; i++) {
        for(int j = 0; j <= time; j++)
            for(int t = j; t <= time; t++)
                dp[i][j] = (dp[i][j]+dp[i-1][t]*inv[t+1]%mod)%mod;
    }
    ll t = 1, tmp = 0;
    for(int i = 0; i <= time; i++) {
        tmp = (tmp+t*dp[k][i]%mod)%mod;
        t = t*p%mod;
    }
    ans = ans*tmp%mod;
}


int main() {
    init();
    cin >> n >> k;
    for(ll i = 2; i*i <= n; i++) {
        if(n % i == 0) {
            ll cnt = 0;
            while(n % i == 0) {
                n /= i;
                cnt++;
            }
            solve(i, cnt);
        }
    }
    if(n != 1)  solve(n, 1);
    cout << ans << endl;
    return 0;
}