Read problem statements in [Hindi], [Bengali], [Mandarin Chinese], [Russian], and [Vietnamese] as well.

You are given an integer sequence $A_{1}, A_{2}, \ldots, A_{N}$. Find the number of integer sequences $X_{1}, X_{2}, \ldots, X_{N}$ such that
$0 ≤ X_{i} < 239$ for each valid $i$
$X_{i} \neq X_{j}$ for each valid $i \neq j$
$(X_{1} \cdot A_{1}) + (X_{2} \cdot A_{2}) + \ldots + (X_{N} \cdot A_{N})$ is divisible by $239$

Since the answer may be large, compute it modulo $998,244,353$.

------  Input ------
The first line of the input contains a single integer $N$.
The second line contains $N$ space-separated integers $A_{1}, A_{2}, \ldots, A_{N}$.

------  Output ------
Print a single line containing one integer ― the number of sequences satisfying all conditions modulo $998,244,353$.

------  Constraints ------
$1 ≤ N ≤ 20$
$0 ≤ A_{i} < 239$ for each valid $i$

------  Subtasks ------
Subtask #1 (10 points): $N ≤ 10$

Subtask #2 (10 points): $N ≤ 15$

Subtask #3 (80 points): original constraints

----- Sample Input 1 ------ 
3

1 1 1
----- Sample Output 1 ------ 
56406

Difficulty: VERY_HARD
Tags: Dynamic programming

Bit manipulation

Combinatorics

Mathematics
Skill types: Dynamic programming

Bit manipulation

Starter code:

Solutions:
from functools import reduce
from operator import mul

def main():
	n = int(input())
	print(count_vecs(list(map(int, input().split()))))
p = 998244353

def count_vecs(A, q=239):
	if A == []:
		return 1
	if A.count(0) > 0:
		nonzeroes = [a for a in A if a != 0]
		return count_vecs(nonzeroes, q) * fall(q - len(nonzeroes), A.count(0)) % p
	ans = fall(q - 1, len(A) - 1)
	if sum(A) % q == 0:
		A_dict = dict()
		for a in A:
			if a not in A_dict:
				A_dict[a] = 0
			A_dict[a] += 1
		A_unique = sorted(A_dict.keys())
		M = [A_dict[a] for a in A_unique]
		for (P, mult) in zp(A_unique, M, q):
			m = len(P)
			mu_coef = (-1) ** (len(A) - m) * prod([fac(k - 1) for k in P])
			ans += mu_coef * pow_mod(q, m - 1) * (q - 1) * mult
			ans %= p
	return ans

def zp(A, M, q=239, accum=[], sym_facs=[]):
	if sum(M) == 0:
		repeats = 1
		for L in zip(*accum):
			repeats *= multinomial_mod_p(L)
			repeats %= p
		for (num, den) in sym_facs:
			repeats *= num * invert_mod(den, p)
			repeats %= p
		yield (list(map(sum, accum)), repeats)
		return
	k = len(M)
	start = 0
	while M[start] == 0:
		start += 1
	M[start] -= 1
	curr = [0] * k
	dirs = [1] * k
	subset_sum = A[start]
	while True:
		i = start
		while i < k and (not 0 <= curr[i] + dirs[i] <= M[i]):
			dirs[i] *= -1
			i += 1
		if i == k:
			break
		curr[i] += dirs[i]
		subset_sum += dirs[i] * A[i]
		subset_sum %= q
		if subset_sum == 0:
			rest = [M[i] - curr[i] for i in range(k)]
			curr[start] += 1
			for P in zp(A, rest, q, accum + [curr], sym_facs + [(curr[start], M[start] + 1)]):
				yield P
			curr[start] -= 1

def multinomial_mod_p(L):
	n = sum(L)
	if n == 0:
		return 1
	L_sorted = sorted(L)
	accum = 1
	pos = 0
	j = 0
	denom = 1
	while pos < len(L) - 1:
		if denom > L_sorted[pos]:
			denom = 1
			pos += 1
		else:
			accum = accum * (n - j) * invert_mod(denom, p)
			accum %= p
			denom += 1
			j += 1
	return accum

def pow_mod(b, n):
	if n <= 1:
		return b ** n % p
	sqrt = pow_mod(b, n // 2)
	ans = sqrt * sqrt % p
	if n % 2:
		ans = ans * b % p
	return ans

def prod(L):
	return reduce(mul, L, 1)

def fall(x, n):
	accum = 1
	for i in range(x, x - n, -1):
		accum *= i
		accum %= p
	return accum

def fac(n):
	return fall(n, n)

def bezout_solver(a, b):
	if a < b:
		(y, x) = bezout_solver(b, a)
		return (x, y)
	elif b + b == b:
		return (a // a, b - b)
	(q, r) = (a // b, a % b)
	(z, x) = bezout_solver(b, r)
	return (x, z - q * x)

def invert_mod(num, modulus):
	return bezout_solver(modulus, num)[1] % modulus
main()