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

heapify(B)

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)

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

heapify(B)
product = [1]
while B:
    n,X = heappop(B)
    product = convolution(product,X)

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)

Submission status: Time Limit Exceeded

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
    elif a > N:
        print(0)
        exit()
    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: Time Limit Exceeded

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)

    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

    facts, finvs = precompute_factorials(max_cnt, MOD)

    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)
    cmb = np.zeros(z, np.int64)
    same_heights.sort()
    pc = -1

    for c in same_heights:
        if pc != c:
            cmb.fill(0)
            for j in range(c // 2 + 1):
                cmb[j] = facts[c] * finvs[j * 2] % MOD * finvs[c - j * 2] % MOD * pairs[j] % MOD
            butterfly(cmb, MOD, sum_e)
            pc = c
        dp *= cmb
        dp %= MOD

    butterfly_inv(dp, MOD, sum_ie)
    iz = pow_mod(z, MOD - 2, MOD)
    result = dp * iz % 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)