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: Wrong Answer

Solution code:
import numpy as np
from collections import Counter
from heapq import *
import sys
input = sys.stdin.buffer.readline
read = sys.stdin.buffer.read
# 10bitごとに3等分して、Toom-3する。
# どうやら50bit弱の整数を答えとするFFTでも精度が怪しいようなので、
# 40bit程度で持てるように10bitごとで3分割した。(このとき結果の最大値は (log2(fft_len) + 10 * 2) bitほど)
# fft_lenが2 ** 30になることは時間計算量の都合からまずないので、10bit分割で精度が十分出る。

import numpy as np

fft_len = 2 ** 20   # == 1048576　答えの配列の長さ以上である必要があり、小さい素数の積であるとよい。(2の累乗など)
MOD = 998244353     # これが2 ** 30程度からそれ以下であることを前提に組んでいる。

def convolve10(As, Bs):
    "As, Bsの要素が10bitに入る(つまり2**10 = 1024以下)なら誤差はまず平気そう。"
    return np.round(np.real(np.fft.irfft(np.fft.rfft(As, fft_len) * np.fft.rfft(Bs, fft_len), fft_len))).astype(np.int64)

def convolve30(As, Bs):
    """
    10bitで3分割してToom-3
    参考: https://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication
    """
    if len(As) + len(Bs) < 100:
        return np.convolve(As, Bs)
    mask = (1 << 10) - 1
    
    m0 = As & mask
    m1 = (As >> 10) & mask
    m2 = As >> 20

    p_tmp = m0 + m2 # Wikipediaではp_0
    p0 = m0
    p1 = p_tmp + m1
    pm1 = p_tmp - m1
    pm2 = ((pm1 + m2) << 1) - m0
    pinf = m2

    n0 = Bs & mask
    n1 = (Bs >> 10) & mask
    n2 = Bs >> 20

    q_tmp = n0 + n2 # Wikiqediaではq_0
    q0 = n0
    q1 = q_tmp + n1
    qm1 = q_tmp - n1
    qm2 = ((qm1 + n2) << 1) - n0
    qinf = n2

    r0   = convolve10(  p0,   q0)
    r1   = convolve10(  p1,   q1)
    rm1  = convolve10( pm1,  qm1)
    rm2  = convolve10( pm2,  qm2)
    rinf = convolve10(pinf, qinf)

    r_0 = r0
    r_4 = rinf
    r_3 = (rm2 - r1) // 3
    r_1 = (r1 - rm1) >> 1
    r_2 = (rm1 - r0)
    r_3 = ((r_2 - r_3) >> 1) + (rinf << 1)
    r_2 += r_1 - r_4
    r_1 -= r_3
    
    ret = ((r_4 << 20) + (r_3 << 10) + r_2) % MOD
    ret = ((ret << 20) + (r_1 << 10) + r_0) % MOD

    return ret


class Arr:

    def __init__(self, arr):
        self.l = len(arr)
        self.arr = arr

    def __lt__(self, other):
        return self.l < other.l


N = int(input())
hs = list(map(int, read().split()))

ch = Counter(hs)
arrays = []
for k, c in ch.items():
    arr = [1]
    while c >= 2:
        arr.append(c * (c - 1) // 2)
        c -= 2
    arrays.append(Arr(np.array(arr)))

heapify(arrays)
while len(arrays) >= 2:
    a1 = heappop(arrays)
    a2 = heappop(arrays)
    heappush(arrays, Arr(convolve30(a1.arr, a2.arr) % MOD))
    
final = arrays[0].arr
# print(final)

MOD = 998244353
table_len = 100010

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)

finv = [0] * table_len
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in range(table_len-1, 0, -1):
    finv[i - 1] = finv[i] * i % MOD
    

ans = 0

for i in range(N + 1):
    tmp = pow(2, MOD - 1 - N + i, MOD) * fac[2 * (N - i)] % MOD * finv[N - i] % MOD
    # print("t", tmp)
    if len(final) > i:
        tmp *= int(final[i])
        tmp %= MOD
    else:
        break
    if i % 2 == 1:
        ans -= tmp
    else:
        ans += tmp
    ans %= MOD
    # print("a", ans)

print(ans % MOD)

Submission status: Wrong Answer

Solution code:
import numpy as np
from collections import Counter
from heapq import *
import sys
input = sys.stdin.buffer.readline
read = sys.stdin.buffer.read
# 10bitごとに3等分して、Toom-3する。
# どうやら50bit弱の整数を答えとするFFTでも精度が怪しいようなので、
# 40bit程度で持てるように10bitごとで3分割した。(このとき結果の最大値は (log2(fft_len) + 10 * 2) bitほど)
# fft_lenが2 ** 30になることは時間計算量の都合からまずないので、10bit分割で精度が十分出る。

import numpy as np

fft_len = 2 ** 20   # == 1048576　答えの配列の長さ以上である必要があり、小さい素数の積であるとよい。(2の累乗など)
MOD = 998244353     # これが2 ** 30程度からそれ以下であることを前提に組んでいる。

def convolve10(As, Bs):
    "As, Bsの要素が10bitに入る(つまり2**10 = 1024以下)なら誤差はまず平気そう。"
    return np.round(np.real(np.fft.irfft(np.fft.rfft(As, fft_len) * np.fft.rfft(Bs, fft_len), fft_len))).astype(np.int64)

def convolve30(As, Bs):
    """
    10bitで3分割してToom-3
    参考: https://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication
    """
    if len(As) + len(Bs) < 100:
        return np.convolve(As, Bs)
    mask = (1 << 10) - 1
    
    m0 = As & mask
    m1 = (As >> 10) & mask
    m2 = As >> 20

    p_tmp = m0 + m2 # Wikipediaではp_0
    p0 = m0
    p1 = p_tmp + m1
    pm1 = p_tmp - m1
    pm2 = ((pm1 + m2) << 1) - m0
    pinf = m2

    n0 = Bs & mask
    n1 = (Bs >> 10) & mask
    n2 = Bs >> 20

    q_tmp = n0 + n2 # Wikiqediaではq_0
    q0 = n0
    q1 = q_tmp + n1
    qm1 = q_tmp - n1
    qm2 = ((qm1 + n2) << 1) - n0
    qinf = n2

    r0   = convolve10(  p0,   q0)
    r1   = convolve10(  p1,   q1)
    rm1  = convolve10( pm1,  qm1)
    rm2  = convolve10( pm2,  qm2)
    rinf = convolve10(pinf, qinf)

    r_0 = r0
    r_4 = rinf
    r_3 = (rm2 - r1) // 3
    r_1 = (r1 - rm1) >> 1
    r_2 = (rm1 - r0)
    r_3 = ((r_2 - r_3) >> 1) + (rinf << 1)
    r_2 += r_1 - r_4
    r_1 -= r_3
    
    ret = ((r_4 << 20) + (r_3 << 10) + r_2) % MOD
    ret = ((ret << 20) + (r_1 << 10) + r_0) % MOD

    return ret


class Arr:

    def __init__(self, arr):
        self.l = len(arr)
        self.arr = arr

    def __lt__(self, other):
        return self.l < other.l


N = int(input())
hs = list(map(int, read().split()))

ch = Counter(hs)
arrays = []
for k, c in ch.items():
    arr = [1]
    while c >= 2:
        arr.append(c * (c - 1) // 2)
        c -= 2
    arrays.append(Arr(np.array(arr)))

heapify(arrays)
while len(arrays) >= 2:
    a1 = heappop(arrays)
    a2 = heappop(arrays)
    heappush(arrays, Arr(convolve30(a1.arr, a2.arr) % MOD))
    
final = arrays[0].arr
print(final)

MOD = 998244353
table_len = 100010

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)

finv = [0] * table_len
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in range(table_len-1, 0, -1):
    finv[i - 1] = finv[i] * i % MOD
    

ans = 0

for i in range(N + 1):
    tmp = pow(2, MOD - 1 - N + i, MOD) * fac[2 * (N - i)] % MOD * finv[N - i] % MOD
    # print("t", tmp)
    if len(final) > i:
        tmp *= int(final[i])
        tmp %= MOD
    else:
        break
    if i % 2 == 1:
        ans -= tmp
    else:
        ans += tmp
    ans %= MOD
    # print("a", ans)

print(ans % MOD)

Submission status: Wrong Answer

Solution code:
import numpy as np
from collections import Counter
from heapq import *
import sys
input = sys.stdin.buffer.readline
read = sys.stdin.buffer.read


def convolve(a, b):
    l = len(a) + len(b)
    return np.round(np.real(np.fft.irfft(np.fft.rfft(a, l) * np.fft.rfft(b, l), l))).astype(np.int64)


class Arr:

    def __init__(self, arr):
        self.l = len(arr)
        self.arr = arr

    def __lt__(self, other):
        return self.l < other.l


N = int(input())
hs = list(map(int, read().split()))

ch = Counter(hs)
arrays = []
for k, c in ch.items():
    arr = [1]
    while c >= 2:
        arr.append(c * (c - 1) // 2)
        c -= 2
    arrays.append(Arr(np.array(arr)))

heapify(arrays)
while len(arrays) >= 2:
    a1 = heappop(arrays)
    a2 = heappop(arrays)
    heappush(arrays, Arr(convolve(a1.arr, a2.arr)))
    
final = arrays[0].arr

MOD = 998244353
table_len = 100010

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)

finv = [0] * table_len
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in range(table_len-1, 0, -1):
    finv[i - 1] = finv[i] * i % MOD
    

ans = 0

for i in range(N + 1):
    tmp = pow(2, MOD - 1 - N + i, MOD) * fac[2 * (N - i)] % MOD * finv[N - i] % MOD
    # print("t", tmp)
    if len(final) > i:
        tmp *= int(final[i])
        tmp %= MOD
    else:
        break
    if i % 2 == 1:
        ans -= tmp
    else:
        ans += tmp
    ans %= MOD
    # print("a", ans)

print(ans % MOD)

Submission status: Wrong Answer

Solution code:
import numpy as np
from collections import Counter
from heapq import *
import sys
input = sys.stdin.buffer.readline
read = sys.stdin.buffer.read


def convolve(a, b):
    l = len(a) + len(b)
    return np.round(np.real(np.fft.irfft(np.fft.rfft(a, l) * np.fft.rfft(b, l), l))).astype(np.int64)


class Arr:

    def __init__(self, arr):
        self.l = len(arr)
        self.arr = arr

    def __lt__(self, other):
        return self.l < other.l


N = int(input())
hs = list(map(int, read().split()))

ch = Counter(hs)
arrays = []
for k, c in ch.items():
    arr = [1]
    while c >= 2:
        arr.append(c * (c - 1) // 2)
        c -= 2
    arrays.append(Arr(np.array(arr)))

heapify(arrays)
while len(arrays) >= 2:
    a1 = heappop(arrays)
    a2 = heappop(arrays)
    heappush(arrays, Arr(convolve(a1.arr, a2.arr)))
    
final = arrays[0].arr
print(final)

MOD = 998244353
table_len = 100010

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)

finv = [0] * table_len
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in range(table_len-1, 0, -1):
    finv[i - 1] = finv[i] * i % MOD
    

ans = 0

for i in range(N + 1):
    tmp = pow(2, MOD - 1 - N + i, MOD) * fac[2 * (N - i)] % MOD * finv[N - i] % MOD
    # print("t", tmp)
    if len(final) > i:
        tmp *= int(final[i])
        tmp %= MOD
    else:
        break
    if i % 2 == 1:
        ans -= tmp
    else:
        ans += tmp
    ans %= MOD
    # print("a", ans)

print(ans % MOD)

Submission status: Wrong Answer

Solution code:
import os
import sys

import numpy as np


# https://github.com/atcoder/ac-library/blob/master/atcoder/convolution.hpp

def solve(inp):
    def bit_length(n):
        x = 0
        while 1 << x < n:
            x += 1
        return x

    def bit_scan_forward(n):
        x = 0
        while n & 1 == 0:
            n >>= 1
            x += 1
        return x

    def pow_mod(x, n, m):
        r = 1
        y = x % m
        while n:
            if n & 1:
                r = (r * y) % m
            y = y * y % m
            n >>= 1
        return r

    def get_primitive_root(m):
        if m == 2: return 1
        if m == 167772161: return 3
        if m == 469762049: return 3
        if m == 754974721: return 11
        if m == 998244353: return 3
        divs = [2]
        x = (m - 1) // 2
        while x & 1 == 0:
            x >>= 1
        i = 3
        while i * i <= x:
            if x % i == 0:
                divs.append(i)
                x //= i
                while x % i == 0:
                    x //= i
            i += 2
        if x > 1:
            divs.append(x)
        g = 2
        while True:
            ok = True
            for d in divs:
                if pow_mod(g, (m - 1) // d, m) == 1:
                    ok = False
                    break
            if ok:
                return g

    def butterfly_prepare(mod, primitive_root):
        sum_e = np.zeros(30, np.int64)
        sum_ie = np.zeros(30, np.int64)
        es = np.zeros(30, np.int64)
        ies = np.zeros(30, np.int64)
        cnt2 = bit_scan_forward(mod - 1)
        e = pow_mod(primitive_root, (mod - 1) >> cnt2, mod)
        ie = pow_mod(e, mod - 2, mod)
        for i in range(cnt2, 1, -1):
            es[i - 2] = e
            ies[i - 2] = ie
            e = e * e % mod
            ie = ie * ie % mod
        now_e = 1
        now_ie = 1
        for i in range(cnt2 - 1):
            sum_e[i] = es[i] * now_e % mod
            sum_ie[i] = ies[i] * now_ie % mod
            now_e = now_e * ies[i] % mod
            now_ie = now_ie * es[i] % mod
        return sum_e, sum_ie

    def butterfly(aaa, mod, sum_e):
        n = aaa.size
        h = bit_length(n)

        for ph in range(1, h + 1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            now = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = aaa[i + offset]
                    r = aaa[i + offset + p] * now % mod
                    aaa[i + offset] = (l + r) % mod
                    aaa[i + offset + p] = (l - r) % mod
                now = now * sum_e[bit_scan_forward(~s)] % mod

    def butterfly_inv(aaa, mod, sum_ie):
        n = aaa.size
        h = bit_length(n)

        for ph in range(h, 0, -1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            inow = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = aaa[i + offset]
                    r = aaa[i + offset + p]
                    aaa[i + offset] = (l + r) % mod
                    aaa[i + offset + p] = ((l - r) * inow) % mod
                inow = inow * sum_ie[bit_scan_forward(~s)] % mod

    def convolution(aaa, bbb, MOD):
        n = aaa.size
        m = bbb.size
        z = 1 << bit_length(n + m - 1)
        raaa = np.zeros(z, np.int64)
        rbbb = np.zeros(z, np.int64)
        raaa[:n] = aaa
        rbbb[:m] = bbb
        primitive_root = get_primitive_root(MOD)
        sum_e, sum_ie = butterfly_prepare(MOD, primitive_root)
        butterfly(raaa, MOD, sum_e)
        butterfly(rbbb, MOD, sum_e)
        ccc = raaa * rbbb % MOD
        butterfly_inv(ccc, MOD, sum_ie)
        iz = pow_mod(z, MOD - 2, MOD)
        result = ccc[:n + m - 1] * iz % MOD
        return result

    def mod_pow(x, a, MOD):
        ret = 1
        cur = x
        while a:
            if a & 1:
                ret = ret * cur % MOD
            cur = cur * cur % MOD
            a >>= 1
        return ret

    def precompute_factorials(n, MOD):
        factorials = np.ones(n + 1, dtype=np.int64)
        for m in range(2, n + 1):
            factorials[m] = factorials[m - 1] * m % MOD
        inversions = np.ones(n + 1, dtype=np.int64)
        inversions[n] = mod_pow(factorials[n], MOD - 2, MOD)
        for m in range(n, 2, -1):
            inversions[m - 1] = inversions[m] * m % MOD
        return factorials, inversions

    MOD = 998244353
    n = inp[0]
    hhh = inp[1:]
    cnt = {0: 0}
    for h in hhh:
        if h in cnt:
            cnt[h] += 1
        else:
            cnt[h] = 1

    same_heights = []
    max_pair = 0
    max_cnt = 0
    for c in cnt.values():
        if c <= 1:
            continue
        max_pair += c // 2
        max_cnt = max(max_cnt, c)
        same_heights.append(c)

    facts, finvs = precompute_factorials(max_cnt, MOD)

    l = len(same_heights)
    z = 1 << bit_length(max_pair + 1)

    dp = np.zeros(z, np.int64)
    dp[0] = 1
    primitive_root = get_primitive_root(MOD)
    sum_e, sum_ie = butterfly_prepare(MOD, primitive_root)
    butterfly(dp, MOD, sum_e)

    for i in range(l):
        c = same_heights[i]
        cmb = np.zeros(z, np.int64)
        for j in range(c // 2 + 1):
            cmb[j] = facts[c] * finvs[j * 2] % MOD * finvs[c - j * 2] % MOD
        butterfly(cmb, MOD, sum_e)
        dp = dp * cmb % MOD

    butterfly_inv(dp, MOD, sum_ie)
    iz = pow_mod(z, MOD - 2, MOD)
    result = dp * iz % MOD

    pairs = np.zeros(n + 1, np.int64)
    pairs[0] = 1
    for i in range(1, n + 1):
        pairs[i] = pairs[i - 1] * (i * 2 - 1) % MOD

    ans = 0
    for i, v in enumerate(result):
        if i & 1:
            ans = (ans - v * pairs[n - i]) % MOD
        else:
            ans = (ans + v * pairs[n - i]) % MOD

    return ans


if sys.argv[-1] == 'ONLINE_JUDGE':
    from numba.pycc import CC

    cc = CC('my_module')
    cc.export('solve', '(i8[:],)')(solve)
    cc.compile()
    exit()

if os.name == 'posix':
    # noinspection PyUnresolvedReferences
    from my_module import solve
else:
    from numba import njit

    solve = njit('(i8[:],)', cache=True)(solve)
    print('compiled', file=sys.stderr)

inp = np.fromstring(sys.stdin.read(), dtype=np.int64, sep=' ')
ans = solve(inp)
print(ans)

Submission status: Wrong Answer

Solution code:
import os
import sys

import numpy as np


# https://github.com/atcoder/ac-library/blob/master/atcoder/convolution.hpp

def solve(inp):
    def bit_length(n):
        x = 0
        while 1 << x < n:
            x += 1
        return x

    def bit_scan_forward(n):
        x = 0
        while n & 1 == 0:
            n >>= 1
            x += 1
        return x

    def pow_mod(x, n, m):
        r = 1
        y = x % m
        while n:
            if n & 1:
                r = (r * y) % m
            y = y * y % m
            n >>= 1
        return r

    def get_primitive_root(m):
        if m == 2: return 1
        if m == 167772161: return 3
        if m == 469762049: return 3
        if m == 754974721: return 11
        if m == 998244353: return 3
        divs = [2]
        x = (m - 1) // 2
        while x & 1 == 0:
            x >>= 1
        i = 3
        while i * i <= x:
            if x % i == 0:
                divs.append(i)
                x //= i
                while x % i == 0:
                    x //= i
            i += 2
        if x > 1:
            divs.append(x)
        g = 2
        while True:
            ok = True
            for d in divs:
                if pow_mod(g, (m - 1) // d, m) == 1:
                    ok = False
                    break
            if ok:
                return g

    def butterfly_prepare(mod, primitive_root):
        sum_e = np.zeros(30, np.int64)
        sum_ie = np.zeros(30, np.int64)
        es = np.zeros(30, np.int64)
        ies = np.zeros(30, np.int64)
        cnt2 = bit_scan_forward(mod - 1)
        e = pow_mod(primitive_root, (mod - 1) >> cnt2, mod)
        ie = pow_mod(e, mod - 2, mod)
        for i in range(cnt2, 1, -1):
            es[i - 2] = e
            ies[i - 2] = ie
            e = e * e % mod
            ie = ie * ie % mod
        now_e = 1
        now_ie = 1
        for i in range(cnt2 - 1):
            sum_e[i] = es[i] * now_e % mod
            sum_ie[i] = ies[i] * now_ie % mod
            now_e = now_e * ies[i] % mod
            now_ie = now_ie * es[i] % mod
        return sum_e, sum_ie

    def butterfly(aaa, mod, sum_e):
        n = aaa.size
        h = bit_length(n)

        for ph in range(1, h + 1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            now = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = aaa[i + offset]
                    r = aaa[i + offset + p] * now % mod
                    aaa[i + offset] = (l + r) % mod
                    aaa[i + offset + p] = (l - r) % mod
                now = now * sum_e[bit_scan_forward(~s)] % mod

    def butterfly_inv(aaa, mod, sum_ie):
        n = aaa.size
        h = bit_length(n)

        for ph in range(h, 0, -1):
            w = 1 << (ph - 1)
            p = 1 << (h - ph)
            inow = 1
            for s in range(w):
                offset = s << (h - ph + 1)
                for i in range(p):
                    l = aaa[i + offset]
                    r = aaa[i + offset + p]
                    aaa[i + offset] = (l + r) % mod
                    aaa[i + offset + p] = ((l - r) * inow) % mod
                inow = inow * sum_ie[bit_scan_forward(~s)] % mod

    def convolution(aaa, bbb, MOD):
        n = aaa.size
        m = bbb.size
        z = 1 << bit_length(n + m - 1)
        raaa = np.zeros(z, np.int64)
        rbbb = np.zeros(z, np.int64)
        raaa[:n] = aaa
        rbbb[:m] = bbb
        primitive_root = get_primitive_root(MOD)
        sum_e, sum_ie = butterfly_prepare(MOD, primitive_root)
        butterfly(raaa, MOD, sum_e)
        butterfly(rbbb, MOD, sum_e)
        ccc = raaa * rbbb % MOD
        butterfly_inv(ccc, MOD, sum_ie)
        iz = pow_mod(z, MOD - 2, MOD)
        result = ccc[:n + m - 1] * iz % MOD
        return result

    def mod_pow(x, a, MOD):
        ret = 1
        cur = x
        while a:
            if a & 1:
                ret = ret * cur % MOD
            cur = cur * cur % MOD
            a >>= 1
        return ret

    def precompute_factorials(n, MOD):
        factorials = np.ones(n + 1, dtype=np.int64)
        for m in range(2, n + 1):
            factorials[m] = factorials[m - 1] * m % MOD
        inversions = np.ones(n + 1, dtype=np.int64)
        inversions[n] = mod_pow(factorials[n], MOD - 2, MOD)
        for m in range(n, 2, -1):
            inversions[m - 1] = inversions[m] * m % MOD
        return factorials, inversions

    MOD = 998244353
    n = inp[0]
    hhh = inp[1:]
    cnt = {0: 0}
    for h in hhh:
        if h in cnt:
            cnt[h] += 1
        else:
            cnt[h] = 1

    same_heights = []
    max_pair = 0
    max_cnt = 0
    for c in cnt.values():
        if c <= 1:
            continue
        max_pair += c // 2
        max_cnt = max(max_cnt, c)
        same_heights.append(c)

    facts, finvs = precompute_factorials(max_cnt, MOD)
    l = len(same_heights)
    dp = np.ones(1, np.int64)
    for i in range(l):
        c = same_heights[i]
        cmb = np.zeros(c // 2 + 1, np.int64)
        for j in range(c // 2 + 1):
            cmb[j] = facts[c] * finvs[j * 2] % MOD * finvs[c - j * 2] % MOD
        dp = convolution(dp, cmb, MOD)

    pairs = np.zeros(n + 1, np.int64)
    pairs[0] = 1
    for i in range(1, n + 1):
        pairs[i] = pairs[i - 1] * (i * 2 - 1) % MOD

    ans = 0
    for i, v in enumerate(dp):
        if i & 1:
            ans = (ans - v * pairs[n - i]) % MOD
        else:
            ans = (ans + v * pairs[n - i]) % MOD

    return ans


if sys.argv[-1] == 'ONLINE_JUDGE':
    from numba.pycc import CC

    cc = CC('my_module')
    cc.export('solve', '(i8[:],)')(solve)
    cc.compile()
    exit()

if os.name == 'posix':
    # noinspection PyUnresolvedReferences
    from my_module import solve
else:
    from numba import njit

    solve = njit('(i8[:],)', cache=True)(solve)
    print('compiled', file=sys.stderr)

inp = np.fromstring(sys.stdin.read(), dtype=np.int64, sep=' ')
ans = solve(inp)
print(ans)

Submission status: Wrong Answer

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

def convolution_all(polys):
    if not polys: return [1]
    polys.sort(key=lambda x:len(x))
    N = len(polys)
    height = (N-1).bit_length()
    N0 = 1<<(height)
    data = [None]*(N0*2)
    data[N0:N0+N] = polys
    data[N0+N:] = [[1] for _ in range(N0-N)]
    
    #print(data)
    #print(polys)
    for k in range(N0-1,0,-1):
        data[k] = convolution(data[2*k], data[2*k+1])        

    return data[1]

SIZE=10**5 + 10 #998244353 #ここを変更する

inv = [0]*SIZE  # inv[j] = j^{-1} mod MOD
fac = [0]*SIZE  # fac[j] = j! mod MOD
finv = [0]*SIZE # finv[j] = (j!)^{-1} mod MOD
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2,SIZE):
    fac[i] = fac[i-1]*i%MOD
finv[-1] = pow(fac[-1],MOD-2,MOD)
for i in range(SIZE-1,0,-1):
    finv[i-1] = finv[i]*i%MOD
    inv[i] = finv[i]*fac[i-1]%MOD

def choose(n,r): # nCk mod MOD の計算
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

dfac = [0]*SIZE  # fac[j] = j! mod MOD
dfac[0] = dfac[1] = 1
for i in range(2,SIZE):
    dfac[i] = dfac[i-2]*i%MOD

import sys
readline = sys.stdin.readline
read = sys.stdin.read

n,*h = map(int, read().split())
from collections import Counter

def get(k): # k 個の頂点からiペア作れる場合の数
    res = [1]
    for i in range(k//2):
        res.append(res[-1]*choose(k-2*i,2)%MOD)
    for i in range(2,k//2 + 1):
        res[i] = res[i]*finv[i]%MOD
    return res        


d = Counter(h)
res = []
for k,v in d.items():
    if v >= 2: res.append(get(v))
    if v > n: 
        print(0)
        exit()
    
bad = convolution_all(res)

#print(bad)
#bad = [1,5,7,3]
ans = 0
sgn = 1
for i,ai in enumerate(bad):
    ans += sgn*(dfac[2*n-2*i-1] if i < n else 1)*ai  
    ans %= MOD
    sgn *= -1
    #print(ans)

print(ans%MOD)
print(1)

Submission status: Wrong Answer

Solution code:
from operator import mul
from functools import reduce
import math

def combinations_count(n, r):
    r = min(r, n - r)
    numer = reduce(mul, range(n, n - r, -1), 1)
    denom = reduce(mul, range(1, r + 1), 1)
    return numer // denom

n = int(input())
dic = {}
for i in range(2*n):
  tmp = int(input())
  dic.setdefault(tmp,0)
  dic[tmp] += 1
print(dic)
dic_max = max(dic.values())
dic_min = min(dic.values())
print(dic_max,dic_min)
print(combinations_count(dic_max, dic_min))
ans = 1
for vl in dic.values():
  print(dic_max,vl)
  ans *= combinations_count(dic_max, vl) * math.factorial(vl)
ans = ans//math.factorial(dic_max)
print(ans)

Submission status: Wrong Answer

Solution code:
import sys
from heapq import *
import numpy as np
import numba
from numba import njit, b1, i4, i8, f8

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 998_244_353

@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), 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[:], 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]

@njit((i8, ), cache=True)
def fact_table(N):
    """Compute 0!, 1!, ..., (N-1)! and their inverses modulo MOD. 

    Parameters
    ----------
    N : int
        length of the output array. 

    Returns
    -------
    fact : np.ndarray
        0!, 1!, ..., (N-1)!

    fact_inv : np.ndarray
        inverse of 0!, 1!, ..., (N-1)! modulo MOD. 
    """
    fact = np.empty(N, np.int64)
    fact[0] = 1
    for n in range(1, N):
        fact[n] = n * fact[n - 1] % MOD
    fact_inv = np.empty(N, np.int64)
    fact_inv[N - 1] = mpow(fact[N - 1], MOD - 2)
    for n in range(N - 1, 0, -1):
        fact_inv[n - 1] = fact_inv[n] * n % MOD
    return fact, fact_inv

@njit((i8[:], ), cache=True)
def main(H):
    max_ntt = 18
    fact, fact_inv = fact_table(100_000)
    roots, iroots = ntt_precompute(max_ntt, 3)

    def f(n, k):
        # n 人居て、k 個のペアを作る
        if k < 0 or 2 * k > n:
            return 0
        x = fact[n] * fact_inv[n - 2 * k] % MOD
        x = x * fact_inv[k] % MOD
        x = x * mpow(2, MOD - 1 - k) % MOD
        return x

    counts = np.bincount(H).astype(np.int64)
    counts = counts[counts > 0]
    N = len(counts)

    q = [(0, 0)] * 0  # 次数、インデックス
    polys = [np.zeros(0, np.int64) for _ in range(N)]
    for k in range(N):
        n = counts[k]
        P = np.empty(n // 2 + 1, np.int64)
        for i in range(n // 2 + 1):
            P[i] = f(n, i)
        polys[k] = P
        heappush(q, (len(P), k))

    while len(q) >= 2:
        _, i = heappop(q)
        _, j = heappop(q)
        P = polys[i]
        Q = polys[j]
        R = ntt_convolve(roots, iroots, P, Q)
        polys[i] = R
        heappush(q, (len(R), i))
    F = polys[q[0][1]]

    N = len(H)
    ans = 0
    for n in range(len(F)):
        x = f(N - 2 * n, N // 2 - n) * F[n] % MOD
        if n & 1:
            x = -x
        ans += x
    return ans % MOD

H = np.array(read().split(), np.int64)[1:]

print(main(H))