Problem:


Score : 600 points

Problem Statement

There are 2N people numbered 1 through 2N.
The height of Person i is h_i.

How many ways are there to make N pairs of people such that the following conditions are satisfied?
Compute the answer modulo 998,244,353.

Each person is contained in exactly one pair.

For each pair, the heights of the two people in the pair are different.

Two ways are considered different if for some p and q, Person p and Person q are paired in one way and not in the other.

Constraints

1 \leq N \leq 50,000

1 \leq h_i \leq 100,000

All values in input are integers.

Input

Input is given from Standard Input in the following format:

N
h_1
:
h_{2N}

Output

Print the answer.

Sample Input 1

2
1
1
2
3

Sample Output 1

2

There are two ways:

Form the pair (Person 1, Person 3) and the pair (Person 2, Person 4).

Form the pair (Person 1, Person 4) and the pair (Person 2, Person 3).

Sample Input 2

5
30
10
20
40
20
10
10
30
50
60

Sample Output 2

516



Sample input:
2
1
1
2
3


Sample output:
2

Submission status: Accepted

Solution code:
import os
import sys
from heapq import heappop, heappush

import numpy as np


# https://github.com/atcoder/ac-library/blob/master/atcoder/convolution.hpp

def solve(inp):
    def ceil_pow2(n):
        x = 0
        while 1 << x < n:
            x += 1
        return x

    def bit_scan_forward(n):
        x = 0
        while n & 1 == 0:
            n >>= 1
            x += 1
        return x

    def pow_mod(x, n, m):
        r = 1
        y = x % m
        while n:
            if n & 1:
                r = (r * y) % m
            y = y * y % m
            n >>= 1
        return r

    def get_primitive_root(m):
        if m == 2: return 1
        if m == 167772161: return 3
        if m == 469762049: return 3
        if m == 754974721: return 11
        if m == 998244353: return 3
        divs = [2]
        x = (m - 1) // 2
        while x & 1 == 0:
            x >>= 1
        i = 3
        while i * i <= x:
            if x % i == 0:
                divs.append(i)
                x //= i
                while x % i == 0:
                    x //= i
            i += 2
        if x > 1:
            divs.append(x)
        g = 2
        while True:
            ok = True
            for d in divs:
                if pow_mod(g, (m - 1) // d, m) == 1:
                    ok = False
                    break
            if ok:
                return g

    def butterfly_prepare(mod, primitive_root):
        sum_e = np.zeros(30, np.int64)
        sum_ie = np.zeros(30, np.int64)
        es = np.zeros(30, np.int64)
        ies = np.zeros(30, np.int64)
        cnt2 = bit_scan_forward(mod - 1)
        e = pow_mod(primitive_root, (mod - 1) >> cnt2, mod)
        ie = pow_mod(e, mod - 2, mod)
        for i in range(cnt2, 1, -1):
            es[i - 2] = e
            ies[i - 2] = ie
            e = e * e % mod
            ie = ie * ie % mod
        now_e = 1
        now_ie = 1
        for i in range(cnt2 - 1):
            sum_e[i] = es[i] * now_e % mod
            sum_ie[i] = ies[i] * now_ie % mod
            now_e = now_e * ies[i] % mod
            now_ie = now_ie * es[i] % mod
        return sum_e, sum_ie

    def butterfly(aaa, mod, sum_e):
        n = aaa.size
        h = ceil_pow2(n)

        for ph in range(1, h + 1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            now = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = aaa[i + offset]
                    r = aaa[i + offset + p] * now % mod
                    aaa[i + offset] = (l + r) % mod
                    aaa[i + offset + p] = (l - r) % mod
                now = now * sum_e[bit_scan_forward(~s)] % mod

    def butterfly_inv(aaa, mod, sum_ie):
        n = aaa.size
        h = ceil_pow2(n)

        for ph in range(h, 0, -1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            inow = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = aaa[i + offset]
                    r = aaa[i + offset + p]
                    aaa[i + offset] = (l + r) % mod
                    aaa[i + offset + p] = ((l - r) * inow) % mod
                inow = inow * sum_ie[bit_scan_forward(~s)] % mod

    MOD = 998244353
    primitive_root = get_primitive_root(MOD)
    sum_e, sum_ie = butterfly_prepare(MOD, primitive_root)

    def convolution(aaa, bbb, MOD):
        n = aaa.size
        m = bbb.size
        z = 1 << ceil_pow2(n + m - 1)
        raaa = np.zeros(z, np.int64)
        rbbb = np.zeros(z, np.int64)
        raaa[:n] = aaa
        rbbb[:m] = bbb
        butterfly(raaa, MOD, sum_e)
        butterfly(rbbb, MOD, sum_e)
        ccc = raaa * rbbb % MOD
        butterfly_inv(ccc, MOD, sum_ie)
        iz = pow_mod(z, MOD - 2, MOD)
        result = ccc[:n + m - 1] * iz % MOD
        return result

    def mod_pow(x, a, MOD):
        ret = 1
        cur = x
        while a:
            if a & 1:
                ret = ret * cur % MOD
            cur = cur * cur % MOD
            a >>= 1
        return ret

    def precompute_factorials(n, MOD):
        factorials = np.ones(n + 1, dtype=np.int64)
        for m in range(2, n + 1):
            factorials[m] = factorials[m - 1] * m % MOD
        inversions = np.ones(n + 1, dtype=np.int64)
        inversions[n] = mod_pow(factorials[n], MOD - 2, MOD)
        for m in range(n, 2, -1):
            inversions[m - 1] = inversions[m] * m % MOD
        return factorials, inversions

    def counter(arr):
        cnt = {0: 0}
        for a in arr:
            if a in cnt:
                cnt[a] += 1
            else:
                cnt[a] = 1
        if cnt[0] == 0:
            del cnt[0]
        return cnt

    n = inp[0]
    hhh = inp[1:]
    hhh_cnt = counter(hhh)

    same_heights = []
    max_pair_cnt = 0
    max_same_height = 0
    for c in hhh_cnt.values():
        if c <= 1:
            continue
        max_pair_cnt += c // 2
        max_same_height = max(max_same_height, c)
        same_heights.append(c)

    pairs = np.zeros(n + 1, np.int64)
    pairs[0] = 1
    for i in range(1, n + 1):
        pairs[i] = pairs[i - 1] * (i * 2 - 1) % MOD

    facts, finvs = precompute_factorials(max_same_height, MOD)

    same_heights.sort()
    pc = -1
    cmbs = [np.ones(1, np.int64)]
    q = [(1, 0)]
    for i, c in enumerate(same_heights):
        if pc == c:
            cmb = cmbs[-1].copy()
        else:
            cmb = np.zeros(c // 2 + 1, np.int64)
            for j in range(c // 2 + 1):
                cmb[j] = facts[c] * finvs[j * 2] % MOD * finvs[c - j * 2] % MOD * pairs[j] % MOD
            pc = c
        cmbs.append(cmb)
        q.append((cmb.size, i + 1))

    while len(q) > 1:
        _, i = heappop(q)
        _, j = heappop(q)
        cmbs[i] = convolution(cmbs[i], cmbs[j], MOD)
        heappush(q, (cmbs[i].size, i))

    result = cmbs[q[0][1]]

    ans = 0
    for i, v in enumerate(result):
        if i & 1:
            ans = (ans - v * pairs[n - i]) % MOD
        else:
            ans = (ans + v * pairs[n - i]) % MOD

    return ans


if sys.argv[-1] == 'ONLINE_JUDGE':
    from numba.pycc import CC

    cc = CC('my_module')
    cc.export('solve', '(i8[:],)')(solve)
    cc.compile()
    exit()

if os.name == 'posix':
    # noinspection PyUnresolvedReferences
    from my_module import solve
else:
    from numba import njit

    solve = njit('(i8[:],)', cache=True)(solve)
    print('compiled', file=sys.stderr)

inp = np.fromstring(sys.stdin.read(), dtype=np.int64, sep=' ')
ans = solve(inp)
print(ans)

Submission status: Accepted

Solution code:
import sys
from heapq import *
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 998_244_353

@njit((i8, i8), cache=True)
def mpow(a, n):
    p = 1
    while n:
        if n & 1:
            p = p * a % MOD
        a = a * a % MOD
        n >>= 1
    return p


@njit((i8, i8), cache=True)
def ntt_precompute(N, primitive_root):
    roots = np.empty(1 << N, np.int64)
    iroots = np.empty(1 << N, np.int64)
    roots[0] = iroots[0] = 1
    for n in range(N):
        x = mpow(primitive_root, (MOD - 1) >> n + 2)
        y = mpow(x, MOD - 2)
        for i in range(1 << n):
            roots[(1 << n) | i] = roots[i] * x % MOD
            iroots[(1 << n) | i] = iroots[i] * y % MOD
    return roots, iroots


@njit((i8[:], i8[:], i8[:], numba.b1), cache=True)
def ntt(roots, iroots, A, inverse):
    N = len(A)
    if not inverse:
        m = N >> 1
        while m:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m] * roots[k] % MOD
                    A[i], A[i + m] = x + y, x - y
            A %= MOD
            m >>= 1
    else:
        m = 1
        while m < N:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m]
                    A[i], A[i + m] = x + y, (x - y) * iroots[k] % MOD
            A %= MOD
            m <<= 1
        invN = mpow(N, MOD - 2)
        for i in range(N):
            A[i] = A[i] * invN % MOD
    return A


@njit((i8[:], i8[:], i8[:], i8[:]), cache=True)
def ntt_convolve(roots, iroots, A, B):
    n = 0
    L = len(A) + len(B) - 1
    while (1 << n) < L:
        n += 1
    A1 = np.zeros(1 << n, np.int64)
    A1[:len(A)] = A
    B1 = np.zeros(1 << n, np.int64)
    B1[:len(B)] = B
    A, B = A1, B1
    ntt(roots, iroots, A, False)
    ntt(roots, iroots, B, False)
    for i in range(len(A)):
        A[i] = A[i] * B[i] % MOD
    ntt(roots, iroots, A, True)
    return A[:L]

@njit((i8, ), cache=True)
def fact_table(N):
    """Compute 0!, 1!, ..., (N-1)! and their inverses modulo MOD. 

    Parameters
    ----------
    N : int
        length of the output array. 

    Returns
    -------
    fact : np.ndarray
        0!, 1!, ..., (N-1)!

    fact_inv : np.ndarray
        inverse of 0!, 1!, ..., (N-1)! modulo MOD. 
    """
    fact = np.empty(N, np.int64)
    fact[0] = 1
    for n in range(1, N):
        fact[n] = n * fact[n - 1] % MOD
    fact_inv = np.empty(N, np.int64)
    fact_inv[N - 1] = mpow(fact[N - 1], MOD - 2)
    for n in range(N - 1, 0, -1):
        fact_inv[n - 1] = fact_inv[n] * n % MOD
    return fact, fact_inv

@njit((i8[:], ), cache=True)
def main(H):
    max_ntt = 18
    fact, fact_inv = fact_table(110_000)
    roots, iroots = ntt_precompute(max_ntt, 3)

    def f(n, k):
        # n 人居て、k 個のペアを作る
        if k < 0 or 2 * k > n:
            return 0
        x = fact[n] * fact_inv[n - 2 * k] % MOD
        x = x * fact_inv[k] % MOD
        x = x * mpow(2, MOD - 1 - k) % MOD
        return x

    counts = np.bincount(H).astype(np.int64)
    counts = counts[counts > 0]
    N = len(counts)

    q = [(0, 0)] * 0  # 次数、インデックス
    polys = [np.zeros(0, np.int64) for _ in range(N)]
    for k in range(N):
        n = counts[k]
        P = np.empty(n // 2 + 1, np.int64)
        for i in range(n // 2 + 1):
            P[i] = f(n, i)
        polys[k] = P
        heappush(q, (len(P), k))

    while len(q) >= 2:
        _, i = heappop(q)
        _, j = heappop(q)
        P = polys[i]
        Q = polys[j]
        R = ntt_convolve(roots, iroots, P, Q)
        polys[i] = R
        heappush(q, (len(R), i))
    F = polys[q[0][1]]
    N = len(H)
    ans = 0
    for n in range(len(F)):
        x = f(N - 2 * n, N // 2 - n) * F[n] % MOD
        if n & 1:
            x = -x
        ans += x
    return ans % MOD

H = np.array(read().split(), np.int64)[1:]

print(main(H))

Submission status: Accepted

Solution code:
ROOT = 3
MOD = 998244353
roots  = [pow(ROOT,(MOD-1)>>i,MOD) for i in range(24)] # 1 の 2^i 乗根
iroots = [pow(x,MOD-2,MOD) for x in roots] # 1 の 2^i 乗根の逆元

def untt(a,n):
    for i in range(n):
        m = 1<<(n-i-1)
        for s in range(1<<i):
            w_N = 1
            s *= m*2
            for p in range(m):
                a[s+p], a[s+p+m] = (a[s+p]+a[s+p+m])%MOD, (a[s+p]-a[s+p+m])*w_N%MOD
                w_N = w_N*roots[n-i]%MOD

def iuntt(a,n):
    for i in range(n):
        m = 1<<i
        for s in range(1<<(n-i-1)):
            w_N = 1
            s *= m*2
            for p in range(m):
                a[s+p], a[s+p+m] = (a[s+p]+a[s+p+m]*w_N)%MOD, (a[s+p]-a[s+p+m]*w_N)%MOD
                w_N = w_N*iroots[i+1]%MOD
            
    inv = pow((MOD+1)//2,n,MOD)
    for i in range(1<<n):
        a[i] = a[i]*inv%MOD

def convolution(a,b):
    la = len(a)
    lb = len(b)
    deg = la+lb-2
    n = deg.bit_length()
    if min(la, lb) <= 50:
        if la < lb:
            la,lb = lb,la
            a,b = b,a
        res = [0]*(la+lb-1)
        for i in range(la):
            for j in range(lb):
                res[i+j] += a[i]*b[j]
                res[i+j] %= MOD
        return res

    N = 1<<n
    a += [0]*(N-len(a))
    b += [0]*(N-len(b))
    untt(a,n)
    untt(b,n)
    for i in range(N):
      a[i] = a[i]*b[i]%MOD
    iuntt(a,n)
    return a[:deg+1]

def convolution_all(polys):
    if not polys: return [1]
    #polys.sort(key=lambda x:len(x))
    N = len(polys)
    height = (N-1).bit_length()
    N0 = 1<<(height)
    data = [None]*(N0*2)
    data[N0:N0+N] = polys
    data[N0+N:] = [[1] for _ in range(N0-N)]
    
    #print(data)
    #print(polys)
    for k in range(N0-1,0,-1):
        data[k] = convolution(data[2*k], data[2*k+1])        

    return data[1]

SIZE=10**5 + 10 #998244353 #ここを変更する

inv = [0]*SIZE  # inv[j] = j^{-1} mod MOD
fac = [0]*SIZE  # fac[j] = j! mod MOD
finv = [0]*SIZE # finv[j] = (j!)^{-1} mod MOD
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2,SIZE):
    fac[i] = fac[i-1]*i%MOD
finv[-1] = pow(fac[-1],MOD-2,MOD)
for i in range(SIZE-1,0,-1):
    finv[i-1] = finv[i]*i%MOD
    inv[i] = finv[i]*fac[i-1]%MOD

def choose(n,r): # nCk mod MOD の計算
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

dfac = [0]*SIZE  # fac[j] = j! mod MOD
dfac[0] = dfac[1] = 1
for i in range(2,SIZE):
    dfac[i] = dfac[i-2]*i%MOD

import sys
readline = sys.stdin.readline
read = sys.stdin.read

n,*h = map(int, read().split())
from collections import Counter

def get(k): # k 個の頂点からiペア作れる場合の数
    res = [1]
    for i in range(k//2):
        res.append(res[-1]*choose(k-2*i,2)%MOD)
    for i in range(2,k//2 + 1):
        res[i] = res[i]*finv[i]%MOD
    return res        


d = Counter(h)
res = []
for k,v in d.items():
    if v >= 2: res.append(get(v))
    if v > n: 
        print(0)
        exit()
    
bad = convolution_all(res)

#print(bad)
#bad = [1,5,7,3]
ans = 0
sgn = 1
for i,ai in enumerate(bad):
    ans += sgn*(dfac[2*n-2*i-1] if i < n else 1)*ai  
    ans %= MOD
    sgn *= -1
    #print(ans)

print(ans%MOD)

Submission status: Accepted

Solution code:
ROOT = 3
MOD = 998244353
roots  = [pow(ROOT,(MOD-1)>>i,MOD) for i in range(24)] # 1 の 2^i 乗根
iroots = [pow(x,MOD-2,MOD) for x in roots] # 1 の 2^i 乗根の逆元

def untt(a,n):
    for i in range(n):
        m = 1<<(n-i-1)
        for s in range(1<<i):
            w_N = 1
            s *= m*2
            for p in range(m):
                a[s+p], a[s+p+m] = (a[s+p]+a[s+p+m])%MOD, (a[s+p]-a[s+p+m])*w_N%MOD
                w_N = w_N*roots[n-i]%MOD

def iuntt(a,n):
    for i in range(n):
        m = 1<<i
        for s in range(1<<(n-i-1)):
            w_N = 1
            s *= m*2
            for p in range(m):
                a[s+p], a[s+p+m] = (a[s+p]+a[s+p+m]*w_N)%MOD, (a[s+p]-a[s+p+m]*w_N)%MOD
                w_N = w_N*iroots[i+1]%MOD
            
    inv = pow((MOD+1)//2,n,MOD)
    for i in range(1<<n):
        a[i] = a[i]*inv%MOD

def convolution(a,b):
    la = len(a)
    lb = len(b)
    deg = la+lb-2
    n = deg.bit_length()
    if min(la, lb) <= 50:
        if la < lb:
            la,lb = lb,la
            a,b = b,a
        res = [0]*(la+lb-1)
        for i in range(la):
            for j in range(lb):
                res[i+j] += a[i]*b[j]
                res[i+j] %= MOD
        return res

    N = 1<<n
    a += [0]*(N-len(a))
    b += [0]*(N-len(b))
    untt(a,n)
    untt(b,n)
    for i in range(N):
      a[i] = a[i]*b[i]%MOD
    iuntt(a,n)
    return a[:deg+1]

def convolution_all(polys):
    if not polys: return [1]
    polys.sort(key=lambda x:len(x))
    N = len(polys)
    height = (N-1).bit_length()
    N0 = 1<<(height)
    data = [None]*(N0*2)
    data[N0:N0+N] = polys
    data[N0+N:] = [[1] for _ in range(N0-N)]
    
    #print(data)
    #print(polys)
    for k in range(N0-1,0,-1):
        data[k] = convolution(data[2*k], data[2*k+1])        

    return data[1]

SIZE=10**5 + 10 #998244353 #ここを変更する

inv = [0]*SIZE  # inv[j] = j^{-1} mod MOD
fac = [0]*SIZE  # fac[j] = j! mod MOD
finv = [0]*SIZE # finv[j] = (j!)^{-1} mod MOD
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2,SIZE):
    fac[i] = fac[i-1]*i%MOD
finv[-1] = pow(fac[-1],MOD-2,MOD)
for i in range(SIZE-1,0,-1):
    finv[i-1] = finv[i]*i%MOD
    inv[i] = finv[i]*fac[i-1]%MOD

def choose(n,r): # nCk mod MOD の計算
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

dfac = [0]*SIZE  # fac[j] = j! mod MOD
dfac[0] = dfac[1] = 1
for i in range(2,SIZE):
    dfac[i] = dfac[i-2]*i%MOD

import sys
readline = sys.stdin.readline
read = sys.stdin.read

n,*h = map(int, read().split())
from collections import Counter

def get(k): # k 個の頂点からiペア作れる場合の数
    res = [1]
    for i in range(k//2):
        res.append(res[-1]*choose(k-2*i,2)%MOD)
    for i in range(2,k//2 + 1):
        res[i] = res[i]*finv[i]%MOD
    return res        


d = Counter(h)
res = []
for k,v in d.items():
    if v >= 2: res.append(get(v))
    if v > n: 
        print(0)
        exit()
    
bad = convolution_all(res)

#print(bad)
#bad = [1,5,7,3]
ans = 0
sgn = 1
for i,ai in enumerate(bad):
    ans += sgn*(dfac[2*n-2*i-1] if i < n else 1)*ai  
    ans %= MOD
    sgn *= -1
    #print(ans)

print(ans%MOD)

Submission status: Accepted

Solution code:
import atexit
import os
import sys
import __pypy__

class NumberTheroemTransform:
  def __init__(self, mod, pr):
    self.mod = mod
    self.pr = pr
    self.M = 2
    self.w = [1]
    self.y = [1]

  def setwy(self, M):
    if M <= self.M:
      return
    self.w += [0] * ((M - self.M) // 2)
    self.y += [0] * ((M - self.M) // 2)
    self.M = M
    z = pow(self.pr, (self.mod - 1) // self.M, self.mod)
    x = pow(z, self.mod - 2, self.mod)
    j = M // 4
    while j:
      self.w[j] = z
      z = z * z % self.mod
      self.y[j] = x
      x = x * x % self.mod
      j //= 2
    self.y[0] = 1
    self.w[0] = 1
    j = self.M // 2
    js = 2
    while js < j:
      z = self.w[js]
      x = self.y[js]
      for k2 in range(js):
        self.w[k2 + js] = self.w[k2] * z % self.mod
        self.y[k2 + js] = self.y[k2] * x % self.mod
      js *= 2

  def fft(self, a):
    mod = self.mod
    self.setwy(len(a))
    u = 1
    v = len(a) >> 1
    while v:
      for j in range(v):
        a[j], a[j + v] = a[j] + a[j + v], a[j] - a[j + v]
        if a[j] >= mod:
          a[j] -= mod
        if a[j + v] < 0:
          a[j + v] += mod
      for jh in range(1, u):
        wj = self.w[jh]
        js = jh * v * 2
        je = js + v
        for j in range(js, je):
          ajv = wj * a[j + v] % mod
          a[j + v] = a[j] - ajv
          if a[j + v] < 0:
            a[j + v] += mod
          a[j] = a[j] + ajv
          if a[j] >= mod:
            a[j] -= mod
      u *= 2
      v >>= 1

  def ifft(self, a):
    mod = self.mod
    self.setwy(len(a))
    u = len(a) >> 1
    v = 1
    while u:
      for j in range(v):
        a[j], a[j + v] = a[j] + a[j + v], a[j] - a[j + v]
        if a[j] >= mod:
          a[j] -= mod
        if a[j + v] < 0:
          a[j + v] += mod
      for jh in range(1, u):
        wj = self.y[jh]
        js = jh * v * 2
        je = js + v
        for j in range(js, je):
          ajv = a[j] - a[j + v]
          if ajv < 0:
            ajv += mod
          a[j] = a[j] + a[j + v]
          if a[j] >= mod:
            a[j] -= mod
          a[j + v] = wj * ajv % mod
      u >>= 1
      v *= 2

  def multiply(self, s, t):
    sl, tl = len(s), len(t)
    L = sl + tl - 1
    M = 2 ** (L - 1).bit_length()
    s += [0] * (M - sl)
    t += [0] * (M - tl)
    self.fft(s)
    self.fft(t)
    for i in range(M):
      s[i] = s[i] * t[i] % self.mod
    self.ifft(s)
    invk = pow(M, self.mod - 2, self.mod)
    for i in range(L):
      s[i] = s[i] * invk % self.mod
    del s[L:]
    return s

ntt = NumberTheroemTransform(998244353, 3)

class Fastio:
  def __init__(self):
    self.ibuf = bytes()
    self.pil = 0
    self.pir = 0
    self.sb = __pypy__.builders.StringBuilder()

  def load(self):
    self.ibuf = self.ibuf[self.pil:]
    self.ibuf += os.read(0, 131072)
    self.pil = 0
    self.pir = len(self.ibuf)

  def flush(self):
    os.write(1, self.sb.build().encode())

  def fastin(self):
    if self.pir - self.pil < 32:
      self.load()
    minus = 0
    x = 0
    while self.ibuf[self.pil] < 45:
      self.pil += 1
    if self.ibuf[self.pil] == 45:
      minus = 1
      self.pil += 1
    while self.ibuf[self.pil] >= 48:
      x = x * 10 + (self.ibuf[self.pil] & 15)
      self.pil += 1
    if minus:
      x = -x
    return x

  def fastout(self, x):
    self.sb.append(str(x))

  def fastoutln(self, x):
    self.sb.append(str(x) + '\n')


fastio = Fastio()
rd = fastio.fastin
wt = fastio.fastout
wtn = fastio.fastoutln
atexit.register(fastio.flush)


MOD = 998244353
M = 100100
f = [0] * M
g = [0] * M
h = [0] * M
f[0] = g[0] = f[1] = g[1] = h[1] = 1
for i in range(2, M):
  f[i] = f[i - 1] * i % MOD
  h[i] = MOD - h[MOD % i] * (MOD // i) % MOD
  g[i] = g[i - 1] * h[i] % MOD

N = rd()
cnt = [0] * M
for i in range(N * 2):
  n = rd()
  cnt[n] += 1

class C:
  def __init__(sl, a):
    sl.a = a

  def __lt__(sl, ot):
    return len(sl.a) < len(ot.a)

lst = [C([1])]
for x in cnt:
  if x < 2:
    continue
  v = [0] * (x // 2 + 1)
  c = 1
  n = x
  m = 0
  for i in range(x // 2 + 1):
    if i != 0:
      c = c * n % MOD * (n - 1) % MOD
      if c & 1:
        c = (c + MOD) // 2
      else:
        c //= 2
      n -= 2
      m += 1
    v[i] = c * g[m] % MOD
  lst.append(C(v))



import heapq

heapq.heapify(lst)
while len(lst) > 1:
  a = heapq.heappop(lst)
  b = heapq.heappop(lst)
  c = ntt.multiply(a.a, b.a)
  heapq.heappush(lst, C(c))
f = lst[0].a

pre = [0] * (N + 1)
n = 2 * N
m = 0
c = 1
for i in range(N + 1):
  if i != 0:
    c = c * (i * 2) % MOD * (i * 2 - 1) % MOD
    c = (c + MOD) >> 1 if c & 1 else c >> 1
    m += 1
  pre[i] = c * g[m] % MOD
ans = 0
for i in range(len(f)):
  c = f[i] * pre[N - i]
  ans += -c if i & 1 else c
print(ans % MOD)

Submission status: Runtime Error

Solution code:
import sys
readline = sys.stdin.readline
from collections import Counter
from heapq import heappop as hpp, heappush as hp


MOD = 998244353
def frac(limit):
    frac = [1]*limit
    for i in range(2,limit):
        frac[i] = i * frac[i-1]%MOD
    fraci = [None]*limit
    fraci[-1] = pow(frac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        fraci[i] = fraci[i+1] * (limit + i + 1) % MOD
    return frac, fraci
frac, fraci = frac(1341398)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return frac[a]*fraci[b]*fraci[a-b]%MOD


pr = 3
LS = 20
class NTT:
    def __init__(self):
        self.N0 = 1<<LS
        omega = pow(pr, (MOD-1)//self.N0, MOD)
        omegainv = pow(omega, MOD-2, MOD)
        self.w = [0]*(self.N0//2)
        self.winv = [0]*(self.N0//2)
        self.w[0] = 1
        self.winv[0] = 1
        self.cache = dict()
        for i in range(1, self.N0//2):
            self.w[i] = (self.w[i-1]*omega)%MOD
            self.winv[i] = (self.winv[i-1]*omegainv)%MOD
        used = set()
        for i in range(self.N0//2):
            if i in used:
                continue
            j = 0
            for k in range(LS-1):
                j |= (i>>k&1) << (LS-2-k)
            used.add(j)
            self.w[i], self.w[j] = self.w[j], self.w[i]
            self.winv[i], self.winv[j] = self.winv[j], self.winv[i]
        
    def _fft(self, A):
        M = len(A)
        bn = 1
        hbs = M>>1
        while hbs:
            for j in range(hbs):
                A[j], A[j+hbs] = A[j] + A[j+hbs], A[j] - A[j+hbs]
                if A[j] > MOD:
                    A[j] -= MOD
                if A[j+hbs] < 0:
                    A[j+hbs] += MOD
            for bi in range(1, bn):
                wi = self.w[bi]
                for j in range(bi*hbs*2, bi*hbs*2+hbs):
                    A[j], A[j+hbs] = (A[j] + wi*A[j+hbs])%MOD, (A[j] - wi*A[j+hbs])%MOD
            bn <<= 1
            hbs >>= 1

    def _ifft(self, A):
        M = len(A)
        bn = M>>1
        hbs = 1
        while bn:
            for j in range(hbs):
                A[j], A[j+hbs] = A[j] + A[j+hbs], A[j] - A[j+hbs]
                if A[j] > MOD:
                    A[j] -= MOD
                if A[j+hbs] < 0:
                    A[j+hbs] += MOD
            for bi in range(1, bn):
                winvi = self.winv[bi]
                for j in range(bi*hbs*2, bi*hbs*2+hbs):
                    A[j], A[j+hbs] = (A[j] + A[j+hbs])%MOD, winvi*(A[j] - A[j+hbs])%MOD
            bn >>= 1
            hbs <<= 1
    def convolve(self, A, B):
        key = tuple(A+[None]+B)
        if key in self.cache:
            return self.cache
        LA = len(A)
        LB = len(B)
        if min(LA, LB) <= 50:
            res = [0]*(LA+LB-1)
            for i in range(LA):
                for j in range(LB):
                    res[i+j] = (res[i+j]+A[i]*B[j])%MOD
            self.cache[key] = res
            return res
            
        
        LC = LA+LB-1
        M = 1<<(LC-1).bit_length()
        A += [0]*(M-LA)
        B += [0]*(M-LB)
        self._fft(A)
        self._fft(B)
        C = [0]*M
        for i in range(M):
            C[i] = A[i]*B[i]%MOD
        self._ifft(C)
        minv = pow(M, MOD-2, MOD)
        for i in range(LC):
            C[i] = C[i]*minv%MOD
        
        res = C[:LC]
        self.cache[key] = res
        return res
        
    
    def inverse(self, A):
        LA = len(A)
        dep = (LA-1).bit_length()
        M = 1<<dep
        A += [0]*(M-LA)
        
        g = [pow(A[0], MOD-2, MOD)]
        for n in range(dep):
            dl = 1<<(n+1)
            f = A[:dl]
            fg = self.convolve(f, g[:])[:dl]
            fgg = self.convolve(fg, g[:])[:dl]
            ng = [None]*dl
            for i in range(dl//2):
                ng[i] = (2*g[i]-fgg[i])%MOD
            for i in range(dl//2, dl):
                ng[i] = MOD-fgg[i]
            g = ng[:]
        return g[:LA]
    
    def integral(self, A):
        A = [1] + [A[i]*fraci[i+2]%MOD for i in range(len(A))]



N = int(readline())
H = [int(readline()) for _ in range(2*N)]
Ch = list(Counter(H).values())
Ch.sort()

table = []
Q = []
M = len(Ch)

for j in range(M):
    c = Ch[j]
    res = [0]*(c//2 + 1)
    for i in range(c//2 + 1):
        res[i] = comb(c, 2*i)*frac[2*i]*fraci[i]%MOD

    table.append(res)
    hp(Q, (c//2, j))


FT = NTT()
for _ in range(M-1):
    c1, i = hpp(Q)
    c2, j = hpp(Q)
    table[i] = FT.convolve(table[i], table[j])
    table[j] = None
    hp(Q, (c1+c2, i))
c, i = hpp(Q)
res = table[i]
ans = 0
for i in range(len(res)):
    ans = (ans + (-1 if i&1 else 1)*res[i]*comb(N, i)*frac[2*N-2*i]*frac[i]) % MOD

ans = ans*pow((MOD+1)//2, N, MOD)*fraci[N]%MOD
print(ans)

Submission status: Runtime Error

Solution code:
import sys
readline = sys.stdin.readline
from collections import Counter, deque


def frac(limit):
    frac = [1]*limit
    for i in range(2,limit):
        frac[i] = i * frac[i-1]%MOD
    fraci = [None]*limit
    fraci[-1] = pow(frac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        fraci[i] = fraci[i+1] * (limit + i + 1) % MOD
    return frac, fraci
frac, fraci = frac(1341398)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return frac[a]*fraci[b]*fraci[a-b]%MOD

#from submission 17076021
MOD = 998244353
ROOT = 3
roots  = [pow(ROOT,(MOD-1)>>i,MOD) for i in range(24)] # 1 の 2^i 乗根
iroots = [pow(x,MOD-2,MOD) for x in roots] # 1 の 2^i 乗根の逆元
 
def untt(a,n):
    for i in range(n):
        m = 1<<(n-i-1)
        for s in range(1<<i):
            w_N = 1
            s *= m*2
            for p in range(m):
                a[s+p], a[s+p+m] = (a[s+p]+a[s+p+m])%MOD, (a[s+p]-a[s+p+m])*w_N%MOD
                w_N = w_N*roots[n-i]%MOD
 
def iuntt(a,n):
    for i in range(n):
        m = 1<<i
        for s in range(1<<(n-i-1)):
            w_N = 1
            s *= m*2
            for p in range(m):
                a[s+p], a[s+p+m] = (a[s+p]+a[s+p+m]*w_N)%MOD, (a[s+p]-a[s+p+m]*w_N)%MOD
                w_N = w_N*iroots[i+1]%MOD
            
    inv = pow((MOD+1)//2,n,MOD)
    for i in range(1<<n):
        a[i] = a[i]*inv%MOD
 
def convolution(a,b):
    la = len(a)
    lb = len(b)
    deg = la+lb-2
    n = deg.bit_length()
    if min(la, lb) <= 50:
        if la < lb:
            la,lb = lb,la
            a,b = b,a
        res = [0]*(la+lb-1)
        for i in range(la):
            for j in range(lb):
                res[i+j] += a[i]*b[j]
                res[i+j] %= MOD
        return res
 
    N = 1<<n
    a += [0]*(N-len(a))
    b += [0]*(N-len(b))
    untt(a,n)
    untt(b,n)
    for i in range(N):
      a[i] = a[i]*b[i]%MOD
    iuntt(a,n)
    return a[:deg+1]

N = int(readline())
H = [int(readline()) for _ in range(2*N)]
Ch = list(Counter(H).values())
Ch.sort()

Q = deque([])
cnt = 0
for c in Ch[::-1]:
    if c == 1:
        continue
    res = [0]*(c//2 + 1)
    for i in range(c//2 + 1):
        res[i] = comb(c, 2*i)*frac[2*i]*fraci[i]%MOD

    Q.append(res)
    cnt += 1

if not cnt:
    print(frac[2*N]*fraci[N]*pow((MOD+1)//2, N, MOD)%MOD)
else:
    M = cnt
    for _ in range(M-1):
        qi = Q.pop()
        qj = Q.pop()
        Q.appendleft(convolution(qi, qj))
    ans = 0
    res = Q[0]
    for i in range(len(res)):
        ans = (ans + (-1 if i&1 else 1)*res[i]*comb(N, i)%MOD*frac[2*N-2*i]*frac[i]) % MOD
    
    ans = ans*pow((MOD+1)//2, N, MOD)*fraci[N]%MOD
    print(ans)