Counting Factorizations
The prime factorization of a positive integer $$$m$$$ is the unique way to write it as $$$\displaystyle m=p_1^{e_1}\cdot p_2^{e_2}\cdot \ldots \cdot p_k^{e_k}$$$, where $$$p_1, p_2, \ldots, p_k$$$ are prime numbers, $$$p_1 < p_2 < \ldots < p_k$$$ and $$$e_1, e_2, \ldots, e_k$$$ are positive integers. For each positive integer $$$m$$$, $$$f(m)$$$ is defined as the multiset of all numbers in its prime factorization, that is $$$f(m)=\{p_1,e_1,p_2,e_2,\ldots,p_k,e_k\}$$$. For example, $$$f(24)=\{2,3,3,1\}$$$, $$$f(5)=\{1,5\}$$$ and $$$f(1)=\{\}$$$.You are given a list consisting of $$$2n$$$ integers $$$a_1, a_2, \ldots, a_{2n}$$$. Count how many positive integers $$$m$$$ satisfy that $$$f(m)=\{a_1, a_2, \ldots, a_{2n}\}$$$. Since this value may be large, print it modulo $$$998\,244\,353$$$.
In the first sample, the two values of $$$m$$$ such that $$$f(m)=\{1,2,3,3\}$$$ are $$$m=24$$$ and $$$m=54$$$. Their prime factorizations are $$$24=2^3\cdot 3^1$$$ and $$$54=2^1\cdot 3^3$$$.In the second sample, the five values of $$$m$$$ such that $$$f(m)=\{2,2,3,5\}$$$ are $$$200, 225, 288, 500$$$ and $$$972$$$.In the third sample, there is no value of $$$m$$$ such that $$$f(m)=\{1,4\}$$$. Neither $$$1^4$$$ nor $$$4^1$$$ are prime factorizations because $$$1$$$ and $$$4$$$ are not primes.
The first line contains one integer $$$n$$$ ($$$1\le n \le 2022$$$).The second line contains $$$2n$$$ integers $$$a_1, a_2, \ldots, a_{2n}$$$ ($$$1\le a_i\le 10^6$$$)  — the given list.
Print one integer, the number of positive integers $$$m$$$ satisfying $$$f(m)=\{a_1, a_2, \ldots, a_{2n}\}$$$ modulo $$$998\,244\,353$$$.
import sys

def gcd(a, b):
    if a == 0:
        return b
    return gcd(b % a, a)

def inv(a):
    g = gcd(a, m)
    if g != 1:
        return -1
    return pow(a, m - 2, m)

def div(a, b):
    a %= m
    pog = inv(b)
    return pog * a % m

def cr(N, K):
    return div(tb[N], tb[K] * tb[N - K])
m = 998244353
mx = int(1000000.0) + 5
sv = [1] * mx
sv[1] = 0
for i in range(2, mx):
    if sv[i]:
        for j in range(i + i, mx, i):
            sv[j] = 0
input = sys.stdin.readline
n = int(input())
a = list(map(int, input().split()))
pms = dict()
cmps = dict()
f = [1] * (n + 2)
for i in range(1, n + 2):
    f[i] = f[i - 1] * i
    f[i] %= m
for i in range(2 * n):
    if sv[a[i]]:
        if a[i] in pms:
            pms[a[i]] += 1
        else:
            pms[a[i]] = 1
    elif a[i] in cmps:
        cmps[a[i]] += 1
    else:
        cmps[a[i]] = 1
if len(pms) < n:
    print(0)
    sys.exit(0)
C = 1
for i in range(1, n + 1):
    C *= i
    C %= m
for i in cmps.keys():
    C = div(C, f[cmps[i]])
for i in pms.keys():
    C = div(C, f[pms[i]])
k = len(pms)
ps = list(pms.values())
t = k - n
dp = [[0] * k for i in range(n)]
for i in range(n):
    if i == 0:
        for j in range(k):
            dp[i][j] = ps[j]
    else:
        c = dp[i - 1][i - 1]
        for j in range(i, k):
            dp[i][j] = ps[j] * c
            c += dp[i - 1][j]
            dp[i][j] %= m
            c %= m
C *= sum(dp[n - 1])
C %= m
print(C)
maxN = 1000005
sieve = list(range(maxN))
for i in range(2, maxN):
    if sieve[i] == i:
        for j in range(2 * i, maxN, i):
            sieve[j] = i

def make_nCr_mod(max_n=2 * 10 ** 5, mod=10 ** 9 + 7):
    max_n = min(max_n, mod - 1)
    (fact, inv_fact) = ([0] * (max_n + 1), [0] * (max_n + 1))
    fact[0] = 1
    for i in range(max_n):
        fact[i + 1] = fact[i] * (i + 1) % mod
    inv_fact[-1] = pow(fact[-1], mod - 2, mod)
    for i in reversed(range(max_n)):
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % mod

    def nCr_mod(n, r):
        res = 1
        while n or r:
            (a, b) = (n % mod, r % mod)
            if a < b:
                return 0
            res = res * fact[a] % mod * inv_fact[b] % mod * inv_fact[a - b] % mod
            n //= mod
            r //= mod
        return res
    return (nCr_mod, fact, inv_fact)
MOD = 998244353
(_, fact, inv_fact) = make_nCr_mod(max_n=3000, mod=MOD)

def is_prime(x):
    if x == 1:
        return False
    return sieve[x] == x

def main():
    n = int(input())
    a = readIntArr()
    unique_p_cnts = dict()
    unique_other_cnts = dict()
    for x in a:
        if is_prime(x):
            if x not in unique_p_cnts.keys():
                unique_p_cnts[x] = 0
            unique_p_cnts[x] += 1
        else:
            if x not in unique_other_cnts.keys():
                unique_other_cnts[x] = 0
            unique_other_cnts[x] += 1
    if len(unique_p_cnts.keys()) < n:
        print(0)
        return
    base_multiplier = fact[n]
    for c in unique_other_cnts.values():
        base_multiplier = base_multiplier * inv_fact[c] % MOD
    arr = list(unique_p_cnts.values())
    m = len(arr)
    arr = list(unique_p_cnts.values())
    m = len(arr)
    dp = [0] * (n + 1)
    dp[0] = 1
    for i in range(m):
        dp2 = [0] * (n + 1)
        for j in range(n + 1):
            if j - 1 >= 0:
                dp2[j] += dp[j - 1] * inv_fact[arr[i] - 1]
            dp2[j] += dp[j] * inv_fact[arr[i]]
            dp2[j] %= MOD
        dp = dp2
    ans = base_multiplier * dp[n] % MOD
    print(ans)
    return
import sys
input = sys.stdin.buffer.readline

def oneLineArrayPrint(arr):
    print(' '.join([str(x) for x in arr]))

def multiLineArrayPrint(arr):
    print('\n'.join([str(x) for x in arr]))

def multiLineArrayOfArraysPrint(arr):
    print('\n'.join([' '.join([str(x) for x in y]) for y in arr]))

def readIntArr():
    return [int(x) for x in input().split()]

def makeArr(defaultValFactory, dimensionArr):
    dv = defaultValFactory
    da = dimensionArr
    if len(da) == 1:
        return [dv() for _ in range(da[0])]
    else:
        return [makeArr(dv, da[1:]) for _ in range(da[0])]

def queryInteractive(a, b):
    print('? {} {}'.format(a, b))
    sys.stdout.flush()
    return int(input())

def answerInteractive(ans):
    print('! {}'.format(ans))
    sys.stdout.flush()
inf = float('inf')
from math import gcd, floor, ceil
import math
for _abc in range(1):
    main()
def redirect_io():
    import pathlib, sys
    fname = pathlib.Path(__file__).parent / 'input.txt'
    sys.stdin = open(fname, 'r')

def ina():
    return list(map(int, input().split()))
from math import sqrt
''

def is_prime(n: int) -> bool:
    if n == 1:
        return False
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    for i in range(3, int(sqrt(n)) + 1, 2):
        if n % i == 0:
            return False
    return True
from collections import defaultdict
mod = 998244353
n = int(input())
arr = ina()
primes_freq = defaultdict(int)
non_prime_freq = defaultdict(int)
for i in arr:
    if is_prime(i):
        primes_freq[i] += 1
    else:
        non_prime_freq[i] += 1
if len(primes_freq.keys()) < n:
    print(0)
    exit()

def factorials(MAXN, mod):
    """"""
    fac = [1]
    for i in range(1, MAXN + 1):
        fac.append(fac[-1] * i % mod)
    invf = [pow(fac[-1], mod - 2, mod)]
    for i in range(MAXN, 0, -1):
        invf.append(i * invf[-1] % mod)
    invf = invf[::-1]
    return (fac, invf)
(fac, invf) = factorials(4045, mod)
sorted_primes = sorted(list(primes_freq.keys()))
''
dp = [[0] * len(sorted_primes) for _ in range(2)]
for prime_idx in range(len(sorted_primes)):
    freq = primes_freq[sorted_primes[prime_idx]]
    if prime_idx == 0:
        dp[0][prime_idx] = invf[freq]
    else:
        dp[0][prime_idx] = invf[freq] * dp[0][prime_idx - 1]
    dp[0][prime_idx] %= mod
remains = primes_freq[sorted_primes[0]] - 1
dp[1][0] = invf[remains]
for num_primes_left in range(1, n + 1):
    for idx in range(num_primes_left - 1, len(sorted_primes)):
        if (num_primes_left, idx) == (1, 0):
            continue
        remains_if_used = primes_freq[sorted_primes[idx]] - 1
        ans = 0
        ans = invf[remains_if_used] * dp[0][idx - 1]
        if num_primes_left - 1 != idx:
            ans += invf[remains_if_used + 1] * dp[1][idx - 1]
        ans %= mod
        dp[1][idx] = ans
    dp[0] = dp[1]
    dp[1] = [0] * len(sorted_primes)
ans = fac[n] * dp[0][-1]
ans %= mod
for j in non_prime_freq.values():
    ans *= invf[j]
    ans %= mod
print(ans)