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: Accepted

Solution code:
# -*- coding: utf-8 -*-

#############
# Libraries #
#############

import sys
input = sys.stdin.readline
sys.setrecursionlimit(1000000)

import math
#from math import gcd
import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache

#############
# Constants #
#############

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

#############
# Functions #
#############

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


#####Shorten#####
def DD(arg): return defaultdict(arg)

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

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

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

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

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

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

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

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

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

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

#############
# Main Code #
#############

ROOT = 3
MOD = 998244353
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 = 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)

Submission status: Accepted

Solution code:
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)] # 1 の 2^i 乗根
    self.iroots = [pow(x,self.MOD-2,self.MOD) for x in self.roots] # 1 の 2^i 乗根の逆元
 
  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)

Submission status: Accepted

Solution code:
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)] # 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(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)

Submission status: Accepted

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

Submission status: Accepted

Solution code:
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)

Submission status: Accepted

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

Submission status: Accepted

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

Submission status: Accepted

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