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 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(2 * n)
    doublefact = [1] * (2 * n + 2)
    for i in range(2, 2 * n + 1):
        doublefact[i] = i * doublefact[i - 2] % MOD

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

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

Solution (Python):
import sys
input = lambda: sys.stdin.readline().rstrip()
from collections import Counter
from heapq import heappush, heappop

p, g, ig = 998244353, 3, 332748118
W = [pow(g, (p - 1) >> i, p) for i in range(24)]
iW = [pow(ig, (p - 1) >> i, p) for i in range(24)]
 
def convolve(a, b):
    def fft(f):
        for l in range(k, 0, -1):
            d = 1 << l - 1
            U = [1]
            for i in range(d):
                U.append(U[-1] * W[l] % p)
 
            for i in range(1 << k - l):
                for j in range(d):
                    s = i * 2 * d + j
                    t = s + d
                    f[s], f[t] = (f[s] + f[t]) % p, U[j] * (f[s] - f[t]) % p
 
    def ifft(f):
        for l in range(1, k + 1):
            d = 1 << l - 1
            U = [1]
            for i in range(d):
                U.append(U[-1] * iW[l] % p)
 
            for i in range(1 << k - l):
                for j in range(d):
                    s = i * 2 * d + j
                    t = s + d
                    f[s], f[t] = (f[s] + f[t] * U[j]) % p, (f[s] - f[t] * U[j]) % p
 
    n0 = len(a) + len(b) - 1
    k = (n0).bit_length()
    n = 1 << k
    a = a + [0] * (n - len(a))
    b = b + [0] * (n - len(b))
    fft(a), fft(b)
    for i in range(n):
        a[i] = a[i] * b[i] % p
    ifft(a)
    invn = pow(n, p - 2, p)
    for i in range(n0):
        a[i] = a[i] * invn % p
    del a[n0:]
    return a

N = int(input())
P = 998244353
nn = 2 * N + 10
fa = [1] * (nn+1)
fainv = [1] * (nn+1)
for i in range(nn):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(nn)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P

dfa = [1] * (nn+1)
for i in range(2, nn+1):
    dfa[i] = dfa[i-2] * i % P
dfa += [1]

X = list(Counter([int(input()) for _ in range(2 * N)]).values())
M = len(X)
Y = [[] for _ in range(M)]
for x in X:
    L = [1]
    s = 1
    for i in range(x // 2):
        s = s * ((x - i * 2) * (x - i * 2 - 1) // 2) % P
        L.append(s * fainv[i+1] % P)
    Y.append(L)

for i in range(1, M)[::-1]:
    Y[i] = convolve(Y[i*2], Y[i*2+1])

ans = 0
for i, x in enumerate(Y[1]):
    ans = (ans + x * (-1 if i & 1 else 1) * dfa[N * 2 - i * 2 - 1]) % P
print(ans)

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

    for i in range(len(F) - 1):
        F.append(convolution(F[2 * i], F[2 * i + 1]))

    f = F.pop()
    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()

Solution (Python):
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, deque
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 = deque()
    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:
        F.append(convolution(F.popleft(), F.popleft()))

    f = F.popleft()
    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()

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