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
readline = sys.stdin.readline
from collections import Counter

MOD = 998244353
def fac(limit):
    fac = [1]*limit
    for i in range(2,limit):
        fac[i] = i * fac[i-1]%MOD
    faci = [None]*limit
    faci[-1] = pow(fac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        faci[i] = faci[i+1] * (limit + i + 1) % MOD
    return fac, faci
fac, faci = fac(234139)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[b]*faci[a-b]%MOD
def perm(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[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
        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):
        LA = len(A)
        LB = len(B)
        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
        return C[:LC]
        
        return C
    
    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]*faci[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())
M = len(Ch)
table = [[] for _ in range(M)]


for c in Ch:
    res = [0]*(c//2 + 1)
    for i in range(c//2 + 1):
        res[i] = perm(c, 2*i)*faci[i]%MOD
        hoge = comb(c, 2*i)*fac[2*i]*faci[i]%MOD
    table.append(res)
 
FT = NTT()


for i in range(M-1, 0, -1):
    table[i] = FT.convolve(table[2*i], table[2*i+1])

ans = 0
for i in range(len(table[1])):
    ans = (ans + (-1 if i&1 else 1)*table[1][i]*comb(N, i)*fac[2*N-2*i]*fac[i]) % MOD
 
ans = ans*pow((MOD+1)//2, N, MOD)*faci[N]%MOD
print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
from collections import Counter

MOD = 998244353
def fac(limit):
    fac = [1]*limit
    for i in range(2,limit):
        fac[i] = i * fac[i-1]%MOD
    faci = [None]*limit
    faci[-1] = pow(fac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        faci[i] = faci[i+1] * (limit + i + 1) % MOD
    return fac, faci
fac, faci = fac(234139)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[b]%MOD*faci[a-b]%MOD
    
def perm(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[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
        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):
        LA = len(A)
        LB = len(B)
        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
        return C[:LC]
        
        return C
    
    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]*faci[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())
M = len(Ch)
table = [[] for _ in range(M)]


for c in Ch:
    res = [0]*(c//2 + 1)
    for i in range(c//2 + 1):
        
        
        res[i] = comb(c, 2*i)*fac[2*i]%MOD*faci[i]%MOD
    table.append(res)
 
FT = NTT()


for i in range(M-1, 0, -1):
    table[i] = FT.convolve(table[2*i], table[2*i+1])

ans = 0
for i in range(len(table[1])):
    ans = (ans + (-1 if i&1 else 1)*table[1][i]*comb(N, i)*fac[2*N-2*i]*fac[i]) % MOD
 
ans = ans*pow((MOD+1)//2, N, MOD)*faci[N]%MOD
print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
from collections import Counter

MOD = 998244353
def fac(limit):
    fac = [1]*limit
    for i in range(2,limit):
        fac[i] = i * fac[i-1]%MOD
    faci = [None]*limit
    faci[-1] = pow(fac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        faci[i] = faci[i+1] * (limit + i + 1) % MOD
    return fac, faci
fac, faci = fac(234139)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[b]*faci[a-b]%MOD
def perm(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[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
        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):
        LA = len(A)
        LB = len(B)
        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
        return C[:LC]
        
        return C
    
    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]*faci[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())
M = len(Ch)
table = [[] for _ in range(M)]


for c in Ch:
    res = [0]*(c//2 + 1)
    for i in range(c//2 + 1):
        res[i] = perm(c, 2*i)*faci[i]%MOD
        
    table.append(res)
 
FT = NTT()


for i in range(M-1, 0, -1):
    table[i] = FT.convolve(table[2*i], table[2*i+1])

ans = 0
for i in range(len(table[1])):
    ans = (ans + (-1 if i&1 else 1)*table[1][i]*comb(N, i)*fac[2*N-2*i]*fac[i]) % MOD
 
ans = ans*pow((MOD+1)//2, N, MOD)*faci[N]%MOD
print(ans)

Solution (Python):
import sys
readline = sys.stdin.readline
from collections import Counter

MOD = 998244353
def fac(limit):
    fac = [1]*limit
    for i in range(2,limit):
        fac[i] = i * fac[i-1]%MOD
    faci = [None]*limit
    faci[-1] = pow(fac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        faci[i] = faci[i+1] * (limit + i + 1) % MOD
    return fac, faci
fac, faci = fac(234139)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[b]*faci[a-b]%MOD
def perm(a, b):
    if not a >= b >= 0:
        return 0
    return fac[a]*faci[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
        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):
        LA = len(A)
        LB = len(B)
        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
        return C[:LC]
        
        return C
    
    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]*faci[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())
M = len(Ch)
table = [[] for _ in range(M)]


for c in Ch:
    res = [0]*(c//2 + 1)
    for i in range(c//2 + 1):
        res[i] = perm(c, 2*i)*faci[i]%MOD
    table.append(res)
 
FT = NTT()


for i in range(M-1, 0, -1):
    table[i] = FT.convolve(table[2*i], table[2*i+1])

ans = 0
for i in range(len(table[1])):
    ans = (ans + (-1 if i&1 else 1)*table[1][i]*comb(N, i)*fac[2*N-2*i]*fac[i]) % MOD
 
ans = ans*pow((MOD+1)//2, N, MOD)*faci[N]%MOD
print(ans)

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

    from collections import deque
    q = deque(polys)
    for _ in range(len(polys)-1):
        a = q.popleft()
        b = q.popleft()
        q.append(convolution(a,b))
    return q[0]        
        

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 sys
input = sys.stdin.readline
sys.setrecursionlimit(1000000)

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 998244353
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]



def DD(arg): return defaultdict(arg)


def inv(n): return pow(n, MOD-2, MOD)


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count





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]
  
N = I()

nokori = [1]
for n in range(1,2*N+10):
  nokori.append((nokori[-1]*(2*n-1)*(2*n)*inv(2)*inv(n))%MOD)

dic = DD(int)
for _ in range(N*2):
  dic[I()] += 1
A = [dic[k] for k in dic]
  
P = []
for a in A:
  temp = 1
  count = 0
  P.append([temp])
  while a >= 2:
    temp *= a*(a-1)//2
    temp %= MOD
    count += 1
    P[-1].append((temp*gyaku_kaijo(count))%MOD)
    a -= 2

Q = deque(P)
while len(Q)>1:
  p = Q.popleft()
  q = Q.popleft()
  Q.append(convolution(p, q))
Q = list(Q[0])
  
ans = 0
M = len(Q)

def sign(x):
  if x%2: return -1
  else: return 1

for i in range(M):
  ans += sign(i) * Q[i] * nokori[N-i]
  ans %= MOD
  
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(1000000)

from collections import defaultdict
from collections import deque
import heapq

MOD = 998244353

def DD(arg): return defaultdict(arg)

def inv(n): return pow(n, MOD-2, MOD)

kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

class NTT_convolution:
  def __init__(self, ROOT=3, MOD=998244353):
    self.ROOT = 3
    self.MOD = 998244353
    self.roots  = [pow(self.ROOT,(self.MOD-1)>>i,self.MOD) for i in range(24)] 
    self.iroots = [pow(x,self.MOD-2,self.MOD) for x in self.roots] 
 
  def untt(self,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])%self.MOD, (a[s+p]-a[s+p+m])*w_N%self.MOD
          w_N = w_N*self.roots[n-i]%self.MOD
 
  def iuntt(self,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)%self.MOD, (a[s+p]-a[s+p+m]*w_N)%self.MOD
          w_N = w_N*self.iroots[i+1]%self.MOD
      
    inv = pow((self.MOD+1)//2,n,self.MOD)
    for i in range(1<<n):
      a[i] = a[i]*inv%self.MOD
 
  def convolution(self,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] %= self.MOD
      return res
 
    N = 1<<n
    a += [0]*(N-len(a))
    b += [0]*(N-len(b))
    self.untt(a,n)
    self.untt(b,n)
    for i in range(N):
      a[i] = a[i]*b[i]%self.MOD
    self.iuntt(a,n)
    return a[:deg+1]

N = int(input())
NTT = NTT_convolution()

nokori = [1]
for n in range(1,2*N+10):
  nokori.append((nokori[-1]*(2*n-1)*(2*n)*inv(2)*inv(n))%MOD)

dic = DD(int)
for _ in range(N*2):
  dic[int(input())] += 1
A = [dic[k] for k in dic]
  
P = []
for a in A:
  temp = 1
  count = 0
  P.append([temp])
  while a >= 2:
    temp *= a*(a-1)//2
    temp %= MOD
    count += 1
    P[-1].append((temp*gyaku_kaijo(count))%MOD)
    a -= 2

Q = []
for i in range(len(P)):
  heapq.heappush(Q,(len(P[i]),i))

while len(Q)>1:
  _,p = heapq.heappop(Q)
  _,q = heapq.heappop(Q)
  temp = NTT.convolution(P[p], P[q])
  P.append(temp)
  heapq.heappush(Q,(len(temp),len(P)-1))
Q = list(P[-1])

ans = 0
M = len(Q)

def sign(x):
  if x%2: return -1
  else: return 1

for i in range(M):
  ans += sign(i) * Q[i] * nokori[N-i]
  ans %= MOD
  
print(ans)

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(1000000)

from collections import defaultdict
from collections import deque
import heapq

MOD = 998244353

def DD(arg): return defaultdict(arg)

def inv(n): return pow(n, MOD-2, MOD)

kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

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]

N = int(input())

nokori = [1]
for n in range(1,2*N+10):
  nokori.append((nokori[-1]*(2*n-1)*(2*n)*inv(2)*inv(n))%MOD)

dic = DD(int)
for _ in range(N*2):
  dic[int(input())] += 1
A = [dic[k] for k in dic]
  
P = []
for a in A:
  temp = 1
  count = 0
  P.append([temp])
  while a >= 2:
    temp *= a*(a-1)//2
    temp %= MOD
    count += 1
    P[-1].append((temp*gyaku_kaijo(count))%MOD)
    a -= 2

Q = []
for i in range(len(P)):
  heapq.heappush(Q,(len(P[i]),i))

while len(Q)>1:
  _,p = heapq.heappop(Q)
  _,q = heapq.heappop(Q)
  temp = convolution(P[p], P[q])
  P.append(temp)
  heapq.heappush(Q,(len(temp),len(P)-1))
Q = list(P[-1])

ans = 0
M = len(Q)

def sign(x):
  if x%2: return -1
  else: return 1

for i in range(M):
  ans += sign(i) * Q[i] * nokori[N-i]
  ans %= MOD
  
print(ans)