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: Runtime Error

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

M = len(table)

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: Runtime Error

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
import numba
from numba import njit, b1, i4, i8, f8

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: Runtime Error

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: Runtime Error

Solution code:

## https://atcoder.jp/contests/abc170/submissions/14416787
class SegmentTree():
    def __init__(self,n,init):
        self.offset=2**((n-1).bit_length())
        self.tree=[init]*(2*self.offset)
        self.init=init
    def update(self,pos,val):
        pos+=self.offset
        self.tree[pos]=val
        while pos>1:
            pos=pos//2
            self.tree[pos]=max(self.tree[2*pos],self.tree[2*pos+1])
    def query(self,l,r):
        l+=self.offset
        r+=self.offset
        ret=self.init
        while l<=r:
            ret=max(ret,self.tree[r])
            r=(r-1)//2
            ret=max(ret,self.tree[l])
            l=(l+1)//2
        return ret

import sys
input=sys.stdin.readline
n,k=map(int,input().split())
a=[int(input()) for _ in range(n)]

st=SegmentTree(300001,0)
##
for ai in a:
    l=max(0,ai-k)
    r=min(300001,ai+k)
    m=st.query(l,r)
    st.update(ai,m+1)

print(st.query(0,300001))

Submission status: Runtime Error

Solution code:

## https://atcoder.jp/contests/abc170/submissions/14416787
class SegmentTree():
    def __init__(self,n,init):
        self.offset=2**((n-1).bit_length())
        self.tree=[init]*(2*self.offset)
        self.init=init
    def update(self,pos,val):
        pos+=self.offset
        self.tree[pos]=val
        while pos>1:
            pos=pos//2
            self.tree[pos]=max(self.tree[2*pos],self.tree[2*pos+1])
    def query(self,l,r):
        l+=self.offset
        r+=self.offset
        ret=self.init
        while l<=r:
            ret=max(ret,self.tree[r])
            r=(r-1)//2
            ret=max(ret,self.tree[l])
            l=(l+1)//2
        return ret

import sys
input=sys.stdin.readline
n,k=map(int,input().split())
a=[int(input()) for _ in range(n)]

st=SegmentTree(300001,0)
##この木の query は開区間でない
for ai in a:
    l=max(0,ai-k)
    r=min(300001,ai+k)
    m=st.query(l,r)
    st.update(ai,m+1)

print(st.query(0,300001))

Submission status: Runtime Error

Solution code:
import sys
INF = 1 << 60
MOD = 10**9 + 7 # 998244353
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

class NTT(object):
    def __init__(self, prime, root):
        self._prime = prime
        self._root = root

    def _fmt(self, A, inverse = False):
        N = len(A)
        logN = (N - 1).bit_length()
        prime = self._prime
        base = pow(self._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(self, f, g):
        N = 1 << (len(f) + len(g) - 2).bit_length()
        prime = self._prime
        Ff, Fg = self._fmt(f + [0] * (N - len(f))), self._fmt(g + [0] * (N - len(g)))
        N_inv = pow(N, prime - 2, prime)
        fg = self._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

ntt = NTT(998244353, 3)

from collections import Counter
def resolve():
    N = int(input())
    H = [int(input()) for _ in range(2 * N)]
    C = Counter(H)

    mf = modfact(N)
    doublefact = [1] * (2 * N + 10)
    doublefact[1] = 1
    doublefact[2] = 2
    for i in range(3, 2 * N + 1):
        doublefact[i] = i * doublefact[i - 2] % MOD

    F = []
    for n in C.values():
        f = []
        for i in range(n // 2 + 1):
            f.append(mf.comb(n, 2 * i) * doublefact[2 * i - 1] % MOD)
        F.append(f)

    while len(F) > 1:
        nF = []
        for i in range(len(F) // 2):
            nF.append(ntt.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 - i) - 1] % MOD
        if i & 1:
            ans -= f[i]
        else:
            ans += f[i]
    ans %= MOD
    print(ans)
resolve()

Submission status: Runtime Error

Solution code:
import sys
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 = []
for i in range(100001):
    a = A[i]
    if a < 2:
        continue
    else:
        B.append([(nCr(a,2*j)*double_fac[j]) % mod for j in range(a//2+1)])

b = len(B)

r = 1
while r < b:
    for i in range(0,b,r*2):
        if i+r < b:
            B[i] = convolution(B[i],B[i+r])
    r *= 2

if b > 0:
    C = B[0]  # B 内の多項式の積
else:
    C = [1]

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

print(ans)

Submission status: Runtime Error

Solution code:

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


import sys
sys.setrecursionlimit(10**7)
def I(): return int(sys.stdin.readline().rstrip())


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 = []
for i in range(100001):
    a = A[i]
    if a < 2:
        continue
    else:
        B.append([(nCr(a,2*j)*double_fac[j]) % mod for j in range(a//2+1)])

b = len(B)

r = 1
while r < b:
    for i in range(0,b,r*2):
        if i+r < b:
            B[i] = convolution(B[i],B[i+r])
    r *= 2

if b > 0:
    C = B[0]  # B 内の多項式の積
else:
    C = [1]

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

print(ans)

Submission status: Runtime Error

Solution code:

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


import sys
sys.setrecursionlimit(10**7)
def I(): return int(sys.stdin.readline().rstrip())


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 = []
for i in range(2*N+1):
    a = A[i]
    if a < 2:
        continue
    else:
        B.append([(nCr(a,2*j)*double_fac[j]) % mod for j in range(a//2+1)])

b = len(B)

r = 1
while r < b:
    for i in range(0,b,r*2):
        if i+r < b:
            B[i] = convolution(B[i],B[i+r])
    r *= 2

if b > 0:
    C = B[0]  # B 内の多項式の積
else:
    C = [1]

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

print(ans)