Codeforces Round 856 (Div2）

Counting Factorizations

任何一个正整数 $m$ 都可以被唯一的分解为 $p_1^{e_1}\cdot p_2^{e_2}\dots p_k^{e_k}$ 的形式。将正整数 $m$ 的唯一质数分解转化为一个长度为 $2k$ 的 可重集合 记为 $f(m)$。

$$f(m)=\{p_1,e_1,p_2,e_2,p_3,e_3,\dots ,p_k,e_k\}$$

给定正整数 $n$ 和一个长度为 $2n$ 的可重集 $A$。求出满足 $f(m)=A$ 的正整数 $m$ 的个数。答案对 $998244353$ 取模。

$1\le n\le 2022$

$1\le a_i\le {10}^6$

题解：欧拉筛 + 组合计数 + $DP$

只有质数才能作为底数，也就是说我们要在$2n$个数中选择$n$个的不同的质数作为底数，多余$n$个的质数和非质数作为指数，设有$s_1$个非质数，每个非质数数量为$b_i$,有$s_2$个质数，每个质数的数量为$c_i$ ，我们在$s_2$个质数中选择n个质数作为底数也就是使对应的数量减$1$，$c_i^{\prime }:=c_i-1$,根据组合计数的原理我们得到答案为

$$\sum \frac{n!}{b_1!b_2!...b_{s_1!}{c}_1^{\prime }{c}_2^{\prime }...c_{s_2}^{\prime }}$$

实际上左半部分贡献是固定的$\frac{n!}{b_1!b_2!...b_{s_1!}}$,我们只需求出$\sum \frac{1}{c_1^{\prime }{c}_2^{\prime }...c_{s_2}^{\prime }}$即可，那么对于这部分贡献我们可以通过$DP$求出

状态表示：$f[i][j]$：在前$i$个质数中选择$j$个作为底数的方案数 $O(n^2)$

状态属性：数量

状态计算：

• 不选择第$i$个质数作为底数

$$f[i][j]=f[i][j]+f[i-1][j]\ast c[i]!,j<=i$$

• 选择第$i$个质数作为底数

$$f[i][j]=f[i][j]+f[i-1][j-1]\ast (c[i]-1)!,j<=i$$

状态初始：$f[0][0]=1$

答案呈现：从$s_2$个质数中选$n$个质数作为底数的方案数，$f[s_2][n]$

我们应该提前利用快速幂来预处理阶乘的逆元$O(nlogn)$，同时利用欧拉筛提前预处理出$1e6$范围内的质数

#include <bits/stdc++.h>
#define Zeoy std::ios::sync_with_stdio(false), std::cin.tie(0), std::cout.tie(0)
#define debug(x) cerr << #x << '=' << x << endl
#define all(x) (x).begin(), (x).end()
#define rson id << 1 | 1
#define lson id << 1
#define int long long
#define mpk make_pair
#define endl '\n'
using namespace std;
typedef unsigned long long ULL;
typedef long long ll;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
const int inf = 0x3f3f3f3f;
const ll INF = 0x3f3f3f3f3f3f3f3f;
const int mod = 998244353;
const double eps = 1e-9;
const int N = 1e6 + 10, M = 5e3 + 10;

int n;
int p[N], vis[N], idx;
int fact[M], inv[M], c[M];
int f[M][M];
unordered_map<int, int> uprime, prime;

int qpow(int a, int b, int p)
{
    int res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % p;
        b >>= 1;
        a = a * a % p;
    }
    return res % p;
}

void get_primes(int n)
{
    vis[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!vis[i])
            p[++idx] = i;
        for (int j = 1; j <= idx && p[j] * i <= n; ++j)
        {
            vis[p[j] * i] = 1;
            if (i % p[j] == 0)
                break;
        }
    }
}

void solve()
{
    cin >> n;
    get_primes(1e6);
    int tot = 0;
    for (int i = 1; i <= 2 * n; ++i)
    {
        int x;
        cin >> x;
        if (!vis[x])
            prime[x]++;
        else
            uprime[x]++;
    }
    fact[0] = inv[0] = 1;
    for (int i = 1; i <= n; ++i)
    {
        fact[i] = fact[i - 1] * i % mod;
        inv[i] = inv[i - 1] * qpow(i, mod - 2, mod) % mod;
    }
    int ans = fact[n];
    for (auto [x, y] : uprime)
        ans = ans * inv[y] % mod;
    int cnt = 0;
    for (auto [x, y] : prime)
        c[++cnt] = y;
    f[0][0] = 1;
    for (int i = 1; i <= cnt; ++i)
    {
        for (int j = 0; j <= i; ++j)
            f[i][j] = (f[i][j] % mod + f[i - 1][j] * inv[c[i]] % mod) % mod;
        for (int j = 1; j <= i; ++j)
            f[i][j] = (f[i][j] % mod + f[i - 1][j - 1] * inv[c[i] - 1] % mod) % mod;
    }
    cout << (ans % mod * f[cnt][n] % mod) % mod << endl;
}
signed main(void)
{
    Zeoy;
    int T = 1;
    // cin >> T;
    while (T--)
    {
        solve();
    }
    return 0;
}