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: Time Limit Exceeded

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(134139)
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)
        key2 = tuple(B+[None]+A)
        if key in self.cache:
            return self.cache[key]
        if key2 in self.cache:
            return self.cache[key2]
        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: Time Limit Exceeded

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

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

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]    

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

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

Q = []
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 = convolution_all(Q)
    ans = 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)

Submission status: Time Limit Exceeded

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

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


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

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)

Submission status: Time Limit Exceeded

Solution code:
import sys
readline = sys.stdin.readline
from collections import Counter, deque
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

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

Submission status: Time Limit Exceeded

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


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

#from submission 17076021
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)*frac[2*N-2*i]*frac[i]) % MOD
    
    ans = ans*pow((MOD+1)//2, N, MOD)*fraci[N]%MOD
    print(ans)

Submission status: Time Limit Exceeded

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


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

#from submission 17076021
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:
    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)*frac[2*N-2*i]*frac[i]) % MOD
    
    ans = ans*pow((MOD+1)//2, N, MOD)*fraci[N]%MOD
    print(ans)

Submission status: Time Limit Exceeded

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

#from submission 17076021
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 = []
cnt = 0
for c in Ch:
    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

    hp(Q, (c//2, 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):
        c1, i = hpp(Q)
        c2, j = hpp(Q)
        hp(Q, (c1+c2, convolution(i, j)))
    ans = 0
    _, res = Q[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: Time Limit Exceeded

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

#from submission 17076021
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()

table = []
Q = []
cnt = 0
for c in Ch:
    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

    table.append(res)
    hp(Q, (c//2, cnt))
    cnt += 1

if not table:
    print(frac[2*N]*fraci[N]*pow((MOD+1)//2, N, MOD)%MOD)
else:
    M = len(table)
    
    for _ in range(M-1):
        c1, i = hpp(Q)
        c2, j = hpp(Q)
        table[i] = convolution(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: Time Limit Exceeded

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