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


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 = []

for c in Ch:
    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)

M = len(Ch)
FT = NTT()
for i in range(M.bit_length()-1):
    s = 1<<i
    for st in range(0, M-s, 2*s):
        table[st] = FT.convolve(table[st], table[st+s])
        table[st+s] = None

res = [1]
for s in table[::-1]:
    if s is not None:
        res = FT.convolve(s, res)
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
input = sys.stdin.readline
sys.setrecursionlimit(1000000)

from collections import defaultdict
from collections import deque
import heapq
from scipy import fft

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]

import numpy as np

#@njit((i8, i8), cache=True)
def ntt_precompute(N, primitive_root):
    roots = np.empty(1 << N, np.int64)
    iroots = np.empty(1 << N, np.int64)
    roots[0] = iroots[0] = 1
    for n in range(N):
        x = mpow(primitive_root, (MOD - 1) >> n + 2)
        y = mpow(x, MOD - 2)
        for i in range(1 << n):
            roots[(1 << n) | i] = roots[i] * x % MOD
            iroots[(1 << n) | i] = iroots[i] * y % MOD
    return roots, iroots

#@njit((i8, i8), cache=True)
def mpow(a, n):
    p = 1
    while n:
        if n & 1:
            p = p * a % MOD
        a = a * a % MOD
        n >>= 1
    return p

#@njit((i8[:], i8[:], i8[:], numba.b1), cache=True)
def ntt(roots, iroots, A, inverse):
    N = len(A)
    if not inverse:
        m = N >> 1
        while m:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m] * roots[k] % MOD
                    A[i], A[i + m] = x + y, x - y
            A %= MOD
            m >>= 1
    else:
        m = 1
        while m < N:
            for k in range(N // m >> 1):
                s = 2 * m * k
                for i in range(s, s + m):
                    x, y = A[i], A[i + m]
                    A[i], A[i + m] = x + y, (x - y) * iroots[k] % MOD
            A %= MOD
            m <<= 1
        invN = mpow(N, MOD - 2)
        for i in range(N):
            A[i] = A[i] * invN % MOD
    return A
 
 
#@njit((i8[:], i8[:], i8[:], i8[:]), cache=True)
def ntt_convolve(roots, iroots, A, B):
    n = 0
    L = len(A) + len(B) - 1
    while (1 << n) < L:
        n += 1
    A1 = np.zeros(1 << n, np.int64)
    A1[:len(A)] = A
    B1 = np.zeros(1 << n, np.int64)
    B1[:len(B)] = B
    A, B = A1, B1
    ntt(roots, iroots, A, False)
    ntt(roots, iroots, B, False)
    for i in range(len(A)):
        A[i] = A[i] * B[i] % MOD
    ntt(roots, iroots, A, True)
    return A[:L]
  
N = int(input())

max_ntt = 18
roots, iroots = ntt_precompute(max_ntt, 3)

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


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

while len(Q)>1:
  _,p = heapq.heappop(Q)
  _,q = heapq.heappop(Q)
  temp = ntt_convolve(roots, iroots, data[p], data[q])
  data.append(temp)
  heapq.heappush(Q,(len(temp),len(data)-1))
Q = list(data[-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: Time Limit Exceeded

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

from collections import defaultdict
from collections import deque
import heapq
from scipy import fft

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]

import numpy as np

def fft_convolve(f, g):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    Lf, Lg = f.shape[-1], g.shape[-1]
    Lg = g.shape[-1]
    L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fh, fl = f >> 15, f & (1 << 15) - 1
    gh, gl = g >> 15, g & (1 << 15) - 1
 
    def conv(f, g):
        Ff = fft.rfft(f, fft_len)
        Fg = fft.rfft(g, fft_len)
        h = fft.irfft(Ff * Fg)
        return np.rint(h)[..., :L].astype(np.int64) % MOD
 
    x = conv(fl, gl)
    z = conv(fh, gh)
    y = conv(fl + fh, gl + gh) - x - z
    return (x + (y << 15) + (z << 30)) % MOD
  
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


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

while len(Q)>1:
  _,p = heapq.heappop(Q)
  _,q = heapq.heappop(Q)
  temp = fft_convolve(data[p], data[q])
  data.append(temp)
  heapq.heappush(Q,(len(temp),len(data)-1))
Q = list(data[-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: Time Limit Exceeded

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

from collections import defaultdict
from collections import deque
import heapq
import numpy as np

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 NumberTheoreticTransformation():
    # 2^n-th root of unity
    ntt_prime_data = {
        17: (3, 4),
        998244353: (31, 23)
    }
    
    def __init__(self, ntt_len, MOD = 998244353):
        """
        ntt_len：2べきが必要
        MOD：準備してあるいくつか限定
        """
        self.N = ntt_len
        N = ntt_len
        self.lgN = N.bit_length() - 1
        self.bit_reverse_index = self.bit_reverse(self.lgN)
        self.MOD = MOD
        self.Ninv = pow(N, MOD - 2, MOD)
        r, d = NumberTheoreticTransformation.ntt_prime_data[MOD]
        assert d >= self.lgN
        r = pow(r, 2 ** (d - self.lgN), MOD) # 1の原子 N 乗根
        powr = self.make_zeta_power(r)
        self.Zeta = [powr[N//2::N//(1<<(i+1))] for i in range(self.lgN)]
        self.ZetaInv = [powr[N//2:0:-N//(1<<(i+1))] for i in range(self.lgN)]
        
    def bit_reverse(self, lgN):
        """
        N = 2 ** lgN までの lgN 桁bit反転
        bit_reverse(3) = np.array([000,100,010,110,001,101,011,111]) = np.array([0,4,2,6,1,5,3,7])
        """
        N = 1 << lgN
        A = np.empty(N, np.int32)
        A[0] = 0
        for i in range(lgN):
            x = 1 << i
            A[:x] <<= 1
            A[x:x+x] = A[:x] + 1
        return A
    
    def make_zeta_power(self, r):
        MOD = self.MOD
        x = np.empty(self.N, np.int64)
        x[0] = 1; x[1] = r
        for n in range(self.lgN):
            x[1 << n: 1 << (n+1)] = x[:1 << n] * (r * x[(1 << n) - 1] % MOD) % MOD
        return x
 
    def ntt(self, A, inverse):
            MOD = self.MOD
            A = A[self.bit_reverse_index]
            N = self.N; lgN = self.lgN
            M = 1
            
            Zeta = self.Zeta if not inverse else self.ZetaInv
            for Z in Zeta:
                A = A.reshape(-1,M+M)
                u = A[:,:M]
                t = A[:,M:]
                t *= Z[None,:]
                t += u
                u <<= 1
                u -= t
                M <<= 1
                A %= MOD
            A = A.ravel()
            if inverse:
                A *= self.Ninv
            A %= MOD
            return A
    
    def convolve(self, A, B):
        MOD = self.MOD
        LA = len(A); LB = len(B)
        N = self.N
        A = np.resize(A, N); A[LA:] = 0
        B = np.resize(B, N); B[LB:] = 0
        fA = self.ntt(A, inverse = False)
        fB = self.ntt(B, inverse = False)
        fA *= fB; fA %= MOD
        return self.ntt(fA, inverse = True)
      
      
N = int(input())
ntt_len = 1 << (N + N).bit_length()
ntt = NumberTheoreticTransformation(ntt_len, MOD)

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

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

while len(Q)>1:
  a,p = heapq.heappop(Q)
  b,q = heapq.heappop(Q)
  temp = ntt.convolve(data[p], data[q])[:a+b-1]
  data.append(temp)
  heapq.heappush(Q,(len(temp),len(data)-1))
Q = list(data[-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: Time Limit Exceeded

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]

import numpy as np

def fft_convolve(f, g):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    Lf, Lg = f.shape[-1], g.shape[-1]
    Lg = g.shape[-1]
    L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fh, fl = f >> 15, f & (1 << 15) - 1
    gh, gl = g >> 15, g & (1 << 15) - 1
 
    def conv(f, g):
        Ff = np.fft.rfft(f, fft_len)
        Fg = np.fft.rfft(g, fft_len)
        h = np.fft.irfft(Ff * Fg)
        return np.rint(h)[..., :L].astype(np.int64) % MOD
 
    x = conv(fl, gl)
    z = conv(fh, gh)
    y = conv(fl + fh, gl + gh) - x - z
    return (x + (y << 15) + (z << 30)) % MOD
  
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


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

while len(Q)>1:
  _,p = heapq.heappop(Q)
  _,q = heapq.heappop(Q)
  temp = fft_convolve(data[p], data[q])
  data.append(temp)
  heapq.heappush(Q,(len(temp),len(data)-1))
Q = list(data[-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: Time Limit Exceeded

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

import numpy as np

def fft_convolve(f, g):
    """
    数列 (多項式) f, g の畳み込みの計算．上下 15 bitずつ分けて計算することで，
    30 bit以下の整数，長さ 250000 程度の数列での計算が正確に行える．
    """
    Lf, Lg = f.shape[-1], g.shape[-1]
    Lg = g.shape[-1]
    L = Lf + Lg - 1
    fft_len = 1 << L.bit_length()
    fh, fl = f >> 15, f & (1 << 15) - 1
    gh, gl = g >> 15, g & (1 << 15) - 1
 
    def conv(f, g):
        Ff = np.fft.rfft(f, fft_len)
        Fg = np.fft.rfft(g, fft_len)
        h = np.fft.irfft(Ff * Fg)
        return np.rint(h)[..., :L].astype(np.int64) % MOD
 
    x = conv(fl, gl)
    z = conv(fh, gh)
    y = conv(fl + fh, gl + gh) - x - z
    return (x + (y << 15) + (z << 30)) % MOD
  
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([np.array(p) for p in P])
while len(Q)>1:
  p = Q.popleft()
  q = Q.popleft()
  Q.append(fft_convolve(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: Time Limit Exceeded

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

N = int(input())
k = 96
P = 998244353
m = int(("1" * 32 + "0" * 64) * 50050, 2)
pa = (1 << 64) - ((1 << 64) % P)
kk = 12
pa1 = (1 << kk * 4 + 16) - ((1 << kk * 4 + 16) % P)
pa2 = (1 << kk * 2 + 24) - ((1 << kk * 2 + 24) % P)
pa3 = (1 << kk + 28) - ((1 << kk + 28) % P)
m1 = int(("1" * (kk * 4 - 16) + "0" * (kk * 4 + 16)) * 50050, 2)
m2 = int(("1" * (kk * 6 - 24) + "0" * (kk * 2 + 24)) * 50050, 2)
m3 = int(("1" * (kk * 7 - 28) + "0" * (kk + 28)) * 50050, 2)
def modP(x):
    x -= ((x & m1) >> kk * 4 + 16) * pa1
    x -= ((x & m2) >> kk * 2 + 24) * pa2
    x -= ((x & m3) >> kk + 28) * pa3
    return x
k = 96
K = 1 << k
nu = lambda L: int("".join([bin(K+a)[-k:] for a in L[::-1]]), 2)
st = lambda n: bin(n)[2:] + "0"
li = lambda s, le: [int(a, 2) if len(a) else 0 for a in [s[-(i+1)*k-1:-i*k-1] for i in range(le)]]

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]
# C = lambda a, b: fa[a] * fainv[b] % P * fainv[a-b] % P if 0 <= b <= a else 0


X = list(Counter([int(input()) for _ in range(2 * N)]).values())
M = len(X)
# print("X =", X)
Y = [0] * 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)
    # print("L =", L)
    # print("nu(L) =", nu(L), li(st(nu(L)), len(L)))
    Y.append(nu(L))
for i in range(1, M)[::-1]:
    Y[i] = modP(Y[i*2] * Y[i*2+1])
# print("Y =", Y)
# print("Y[1] =", [a % P for a in li(st(Y[1]), N+1)])
ans = 0
for i, x in enumerate(li(st(Y[1]), N+1)):
    ans = (ans + x * (-1 if i & 1 else 1) * dfa[N * 2 - i * 2 - 1]) % P
print(ans)

Submission status: Time Limit Exceeded

Solution code:
def main():
    import sys
    from heapq import heappop, heappush, heapify
    sys.setrecursionlimit(10 ** 7)

    def I():
        return int(sys.stdin.readline().rstrip())

    # 998244353 = 119*2**23+1

    mod = 998244353
    primitive_root = 3  # mod の原始根
    roots = [pow(primitive_root, (mod - 1) >> i, mod) for i in range(24)]
    inv_roots = [pow(r, mod - 2, mod) for r in roots]

    # roots[i] = 1 の 2**i 乗根、inv_roots[i] = 1 の 2**i 乗根の逆元

    # 順番は変わる

    def ntt(A, n):
        for i in range(n):
            m = 1 << (n - i - 1)
            for start in range(1 << i):
                w = 1
                start *= m * 2
                for j in range(m):
                    A[start + j], A[start + j + m] = (A[start + j] + A[start + j + m]) % mod, (
                                A[start + j] - A[start + j + m]) * w % mod
                    w *= roots[n - i]
                    w %= mod
        return A

    def inv_ntt(A, n):
        for i in range(n):
            m = 1 << i
            for start in range(1 << (n - i - 1)):
                w = 1
                start *= m * 2
                for j in range(m):
                    A[start + j], A[start + j + m] = (A[start + j] + A[start + j + m] * w) % mod, (
                                A[start + j] - A[start + j + m] * w) % mod
                    w *= inv_roots[i + 1]
                    w %= mod
        a = pow(2, n * (mod - 2), mod)
        for i in range(1 << n):
            A[i] *= a
            A[i] %= mod
        return A

    def convolution(A, B):
        a, b = len(A), len(B)
        deg = a + b - 2
        n = deg.bit_length()
        N = 1 << n
        A += [0] * (N - a)  # A の次数を 2冪-1 にする
        B += [0] * (N - b)  # B の次数を 2冪-1 にする
        A = ntt(A, n)
        B = ntt(B, n)
        C = [(A[i] * B[i]) % mod for i in range(N)]
        C = inv_ntt(C, n)
        return C[:deg + 1]

    N = I()
    A = [0] * 100001  # A[i] = 身長 i の人が何人いるか
    for _ in range(2 * N):
        A[I()] += 1

    fac = [1]  # fac[i] = i!
    for i in range(1, N + 1):
        fac.append((fac[-1] * i) % mod)

    fac_inv = []  # fac_inv[i] = (i!)**-1
    for i in range(N + 1):
        fac_inv.append(pow(fac[i], mod - 2, mod))

    double_fac = [1]  # double[i] = (2*i-1)!!
    for i in range(1, N + 1):
        double_fac.append((double_fac[-1] * (2 * i - 1)) % mod)

    def nCr(n, r):
        if n < r:
            return 0
        return (fac[n] * fac_inv[r] * fac_inv[n - r]) % mod

    B = [(0, [1])]
    heapify(B)

    for i in range(100001):
        a = A[i]
        if a < 2:
            continue
        elif a > N:
            print(0)
            exit()
        else:
            heappush(B, (a // 2, [(nCr(a, 2 * j) * double_fac[j]) % mod for j in range(a // 2 + 1)]))

    while len(B) > 1:
        n0, X0 = heappop(B)
        n1, X1 = heappop(B)
        heappush(B, (n0 + n1, convolution(X0, X1)))

    product = heappop(B)[1]

    ans = 0
    r = 1
    for i in range(min(len(product), N + 1)):
        ans += product[i] * double_fac[N - i] * r
        ans %= mod
        r = -r

    print(ans)


if __name__ == '__main__':
    main()

Submission status: Time Limit Exceeded

Solution code:
import sys
from heapq import heappop,heappush,heapify
sys.setrecursionlimit(10**7)
def I(): return int(sys.stdin.readline().rstrip())


# 998244353 = 119*2**23+1

mod = 998244353
primitive_root = 3  # mod の原始根
roots = [pow(primitive_root,(mod-1) >> i,mod) for i in range(24)]
inv_roots = [pow(r,mod-2,mod) for r in roots]
# roots[i] = 1 の 2**i 乗根、inv_roots[i] = 1 の 2**i 乗根の逆元

# 順番は変わる


def ntt(A,n):
    for i in range(n):
        m = 1 << (n-i-1)
        for start in range(1 << i):
            w = 1
            start *= m*2
            for j in range(m):
                A[start+j],A[start+j+m] = (A[start+j]+A[start+j+m]) % mod,(A[start+j]-A[start+j+m])*w % mod
                w *= roots[n-i]
                w %= mod
    return A


def inv_ntt(A,n):
    for i in range(n):
        m = 1 << i
        for start in range(1 << (n-i-1)):
            w = 1
            start *= m*2
            for j in range(m):
                A[start+j],A[start+j+m] = (A[start+j]+A[start+j+m]*w) % mod,(A[start+j]-A[start+j+m]*w) % mod
                w *= inv_roots[i+1]
                w %= mod
    a = pow(2,n*(mod-2),mod)
    for i in range(1 << n):
        A[i] *= a
        A[i] %= mod
    return A


def convolution(A,B):
    a,b = len(A),len(B)
    deg = a+b-2
    n = deg.bit_length()
    N = 1 << n
    A += [0]*(N-a)  # A の次数を 2冪-1 にする
    B += [0]*(N-b)  # B の次数を 2冪-1 にする
    A = ntt(A,n)
    B = ntt(B,n)
    C = [(A[i]*B[i]) % mod for i in range(N)]
    C = inv_ntt(C,n)
    return C[:deg+1]


N = I()
A = [0]*100001  # A[i] = 身長 i の人が何人いるか
for _ in range(2*N):
    A[I()] += 1


fac = [1]  # fac[i] = i!
for i in range(1,N+1):
    fac.append((fac[-1]*i) % mod)

fac_inv = []  # fac_inv[i] = (i!)**-1
for i in range(N+1):
    fac_inv.append(pow(fac[i],mod-2,mod))

double_fac = [1]  # double[i] = (2*i-1)!!
for i in range(1,N+1):
    double_fac.append((double_fac[-1]*(2*i-1)) % mod)


def nCr(n,r):
    if n < r:
        return 0
    return (fac[n]*fac_inv[r]*fac_inv[n-r]) % mod


B = [(0,[1])]
heapify(B)

for i in range(100001):
    a = A[i]
    if a < 2:
        continue
    elif a > N:
        print(0)
        exit()
    else:
        heappush(B,(a//2,[(nCr(a,2*j)*double_fac[j]) % mod for j in range(a//2+1)]))


while len(B) > 1:
    n0,X0 = heappop(B)
    n1,X1 = heappop(B)
    heappush(B,(n0+n1,convolution(X0,X1)))

product = heappop(B)[1]

ans = 0
r = 1
for i in range(min(len(product),N+1)):
    ans += product[i]*double_fac[N-i]*r
    ans %= mod
    r = -r

print(ans)