1. 程式人生 > >Codeforces Hello 2019 D

Codeforces Hello 2019 D

題目連線:https://codeforces.com/contest/1097/problem/D
首先,對每一類質數分開討論,現在問題變成了,設\(n\)中質因子\(p\)的次數是\(i\),對於一個數\(x=p^j\),經過\(k\)輪之後它被操作出來的概率。顯然這跟\(p\)是多少無關,所以記上述情況為\(F_{ij}\)\(F_{ij}\)並不好直接計算,所以考慮令\(f_{ilj}\)表示\(n\)中質因子\(p\)的次數是\(i\),對於一個數\(x=p^j\),經過\(l\)輪之後它被操作出來的概率。\(f_{ilj}\)有如下轉移
\[f_{ili}=1 (l=0)\]
\[f_{ilj}=\sum_{k=j}^{i}f_{il-1k} \times \frac{1}{k+1} (l>1)\]


第二個轉移式子是一段連續區間的和顯然可以字首和優化。注意到\(l\)是由\(l-1\)轉移過來的,故空間上也可以用滾動陣列優化。接下來深搜列舉\(n\)的約數分解後的各個質數的次數並計算答案即可。
時間複雜度\(O(k\log_{2}^{2}n+\sqrt n)\)

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <iostream>
#include <algorithm>
#include <stack>
#include <queue>
#include <map>
#include <set>
#include <iomanip>
#include <assert.h>
#include <fstream>

using namespace std;

typedef long long ll;

const int MAXP = 55;
const int MAXK = 10005;
const ll MOD = 1000000007;

int k,tot;
int cnt[MAXP];

ll n,ans;
ll p[MAXP];
ll inv[MAXP];
ll f[MAXP][MAXP];
ll g[MAXP][MAXP];
ll h[MAXP][MAXP];

ll power(ll a,ll b)
{
    ll res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return res;
}

void init()
{
    for (int i = 1;i <= 50;i++)
        inv[i] = power(i,MOD - 2);
    for (int i = 0;i <= 50;i++)
        f[i][i] = 1;
    for (int i = 0;i <= 50;i++)
    {
        h[i][0] = f[i][0] * inv[1] % MOD;
        for (int j = 1;j <= i;j++)
            h[i][j] = (h[i][j - 1] + f[i][j] * inv[j + 1]) % MOD;
    }
    for (int i = 1;i <= k;i++)
    {
        for (int j = 0;j <= 50;j++)
            for (int k = 0;k <= j;k++)
                g[j][k] = (h[j][j] - (k > 0 ? h[j][k - 1] : 0)) % MOD;
        memcpy(f,g,sizeof(f));
        memset(g,0,sizeof(g));
        for (int j = 0;j <= 50;j++)
        {
            h[j][0] = f[j][0] * inv[1] % MOD;
            for (int k = 1;k <= j;k++)
                h[j][k] = (h[j][k - 1] + f[j][k] * inv[k + 1]) % MOD;
        }
    }
}

void dfs(int now,ll N,ll P)
{
    if (now > tot)
    {
        (ans += N % MOD * P) %= MOD;
        return;
    }
    for (int i = 0;i <= cnt[now];i++)
    {
        dfs(now + 1,N,P * f[cnt[now]][i] % MOD);
        N *= p[now];
    }
}

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