Buggy python submission (verdict: Runtime Error):
import sys
INF = 1 << 60
MOD = 10**9 + 7 # 998244353
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

class NTT(object):
    def __init__(self, prime, root):
        self._prime = prime
        self._root = root

    def _fmt(self, A, inverse = False):
        N = len(A)
        logN = (N - 1).bit_length()
        prime = self._prime
        base = pow(self._root, (prime - 1) // N * (1 - 2 * inverse) % (prime - 1), prime)
        step = N
        for k in range(logN):
            step >>= 1
            w = pow(base, step, prime)
            wj = 1
            nA = [0] * N
            for j in range(1 << k):
                for i in range(1 << logN - k - 1):
                    s, t = i + j * step, i + j * step + (N >> 1)
                    ps, pt = i + j * step * 2, i + j * step * 2 + step
                    nA[s], nA[t] = (A[ps] + A[pt] * wj) % prime, (A[ps] - A[pt] * wj) % prime
                wj = (wj * w) % prime
            A = nA
        return A

    def convolution(self, f, g):
        N = 1 << (len(f) + len(g) - 2).bit_length()
        prime = self._prime
        Ff, Fg = self._fmt(f + [0] * (N - len(f))), self._fmt(g + [0] * (N - len(g)))
        N_inv = pow(N, prime - 2, prime)
        fg = self._fmt([a * b % prime * N_inv % prime for a, b in zip(Ff, Fg)], inverse = True)
        del fg[len(f) + len(g) - 1:]
        return fg

class modfact(object):
    def __init__(self, n):
        fact, invfact = [1] * (n + 1), [1] * (n + 1)
        for i in range(1, n + 1): fact[i] = i * fact[i - 1] % MOD
        invfact[n] = pow(fact[n], MOD - 2, MOD)
        for i in range(n - 1, -1, -1): invfact[i] = invfact[i + 1] * (i + 1) % MOD
        self._fact, self._invfact = fact, invfact

    def fact(self, n):
        return self._fact[n]

    def invfact(self, n):
        return self._invfact[n]

    def comb(self, n, k):
        return self._fact[n] * self._invfact[k] % MOD * self._invfact[n - k] % MOD if 0 <= k <= n else 0

    def perm(self, n, k):
        return self._fact[n] * self._invfact[n - k] % MOD if 0 <= k <= n else 0

ntt = NTT(998244353, 3)

from collections import Counter
def resolve():
    N = int(input())
    H = [int(input()) for _ in range(2 * N)]
    C = Counter(H)

    mf = modfact(N)
    doublefact = [1] * (2 * N + 10)
    doublefact[1] = 1
    doublefact[2] = 2
    for i in range(3, 2 * N + 1):
        doublefact[i] = i * doublefact[i - 2] % MOD

    F = []
    for n in C.values():
        f = []
        for i in range(n // 2 + 1):
            f.append(mf.comb(n, 2 * i) * doublefact[2 * i - 1] % MOD)
        F.append(f)

    while len(F) > 1:
        nF = []
        for i in range(len(F) // 2):
            nF.append(ntt.convolution(F[2 * i], F[2 * i + 1]))
        if len(F) & 1:
            nF.append(F[-1])
        F = nF

    f = F[0]
    ans = 0
    for i in range(len(f)):
        f[i] = f[i] * doublefact[2 * (N - i) - 1] % MOD
        if i & 1:
            ans -= f[i]
        else:
            ans += f[i]
    ans %= MOD
    print(ans)
resolve()

Fixed submission:
import sys
INF = 1 << 60
MOD = 998244353
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

prime = 998244353
root = 3
def _fmt(A, inverse = False):
    N = len(A)
    logN = (N - 1).bit_length()
    base = pow(root, (prime - 1) // N * (1 - 2 * inverse) % (prime - 1), prime)
    step = N
    for k in range(logN):
        step >>= 1
        w = pow(base, step, prime)
        wj = 1
        nA = [0] * N
        for j in range(1 << k):
            for i in range(1 << logN - k - 1):
                s, t = i + j * step, i + j * step + (N >> 1)
                ps, pt = i + j * step * 2, i + j * step * 2 + step
                nA[s], nA[t] = (A[ps] + A[pt] * wj) % prime, (A[ps] - A[pt] * wj) % prime
            wj = (wj * w) % prime
        A = nA
    return A

def convolution(f, g):
    N = 1 << (len(f) + len(g) - 2).bit_length()
    Ff, Fg = _fmt(f + [0] * (N - len(f))), _fmt(g + [0] * (N - len(g)))
    N_inv = pow(N, prime - 2, prime)
    fg = _fmt([a * b % prime * N_inv % prime for a, b in zip(Ff, Fg)], inverse = True)
    del fg[len(f) + len(g) - 1:]
    return fg

class modfact(object):
    def __init__(self, n):
        fact, invfact = [1] * (n + 1), [1] * (n + 1)
        for i in range(1, n + 1): fact[i] = i * fact[i - 1] % MOD
        invfact[n] = pow(fact[n], MOD - 2, MOD)
        for i in range(n - 1, -1, -1): invfact[i] = invfact[i + 1] * (i + 1) % MOD
        self._fact, self._invfact = fact, invfact

    def fact(self, n):
        return self._fact[n]

    def invfact(self, n):
        return self._invfact[n]

    def comb(self, n, k):
        return self._fact[n] * self._invfact[k] % MOD * self._invfact[n - k] % MOD if 0 <= k <= n else 0

    def perm(self, n, k):
        return self._fact[n] * self._invfact[n - k] % MOD if 0 <= k <= n else 0

from collections import Counter
def resolve():
    n = int(input())
    H = [int(input()) for _ in range(2 * n)]

    mf = modfact(n << 1)
    doublefact = [1] * (2 * n + 10)
    for i in range(2, 2 * n + 1):
        doublefact[i] = i * doublefact[i - 2] % MOD

    F = []
    for k in Counter(H).values():
        f = []
        for i in range(k // 2 + 1):
            f.append(mf.comb(k, 2 * i) * doublefact[2 * i - 1] % MOD)
        F.append(f)

    while len(F) > 1:
        nF = []
        for i in range(len(F) // 2):
            nF.append(convolution(F[2 * i], F[2 * i + 1]))
        if len(F) & 1:
            nF.append(F[-1])
        F = nF

    f = F[0]
    ans = 0
    for i in range(len(f)):
        f[i] = f[i] * doublefact[2 * n - 2 * i - 1] % MOD
        if i & 1:
            ans -= f[i]
        else:
            ans += f[i]
    ans %= MOD
    print(ans)
resolve()