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

Solution (Python):
import os
import sys
from heapq import heappop, heappush

import numpy as np




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':
    
    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)

Solution (Python):
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):
        
        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))

Solution (Python):
ROOT = 3
MOD = 998244353
roots  = [pow(ROOT,(MOD-1)>>i,MOD) for i in range(24)] 
iroots = [pow(x,MOD-2,MOD) for x in roots] 

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]
    
    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)]
    
    
    
    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 

inv = [0]*SIZE  
fac = [0]*SIZE  
finv = [0]*SIZE 
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): 
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

dfac = [0]*SIZE  
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): 
    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)



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%MOD)

Solution (Python):
ROOT = 3
MOD = 998244353
roots  = [pow(ROOT,(MOD-1)>>i,MOD) for i in range(24)] 
iroots = [pow(x,MOD-2,MOD) for x in roots] 

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)]
    
    
    
    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 

inv = [0]*SIZE  
fac = [0]*SIZE  
finv = [0]*SIZE 
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): 
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

dfac = [0]*SIZE  
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): 
    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)



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%MOD)

Solution (Python):
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)