This is the easy version of the problem. The only difference between the easy version and the hard version is the constraints on $n$. You can only make hacks if both versions are solved.

A permutation of $1, 2, \ldots, n$ is a sequence of $n$ integers, where each integer from $1$ to $n$ appears exactly once. For example, $[2,3,1,4]$ is a permutation of $1, 2, 3, 4$, but $[1,4,2,2]$ isn't because $2$ appears twice in it.

Recall that the number of inversions in a permutation $a_1, a_2, \ldots, a_n$ is the number of pairs of indices $(i, j)$ such that $i < j$ and $a_i > a_j$.

Let $p$ and $q$ be two permutations of $1, 2, \ldots, n$. Find the number of permutation pairs $(p,q)$ that satisfy the following conditions:

$p$ is lexicographically smaller than $q$.

the number of inversions in $p$ is greater than the number of inversions in $q$.

Print the number of such pairs modulo $mod$. Note that $mod$ may not be a prime.


-----Input-----

The only line contains two integers $n$ and $mod$ ($1\le n\le 50$, $1\le mod\le 10^9$).


-----Output-----

Print one integer, which is the answer modulo $mod$.


-----Examples-----

Input
4 403458273
Output
17


-----Note-----

The following are all valid pairs $(p,q)$ when $n=4$.

$p=[1,3,4,2]$, $q=[2,1,3,4]$,

$p=[1,4,2,3]$, $q=[2,1,3,4]$,

$p=[1,4,3,2]$, $q=[2,1,3,4]$,

$p=[1,4,3,2]$, $q=[2,1,4,3]$,

$p=[1,4,3,2]$, $q=[2,3,1,4]$,

$p=[1,4,3,2]$, $q=[3,1,2,4]$,

$p=[2,3,4,1]$, $q=[3,1,2,4]$,

$p=[2,4,1,3]$, $q=[3,1,2,4]$,

$p=[2,4,3,1]$, $q=[3,1,2,4]$,

$p=[2,4,3,1]$, $q=[3,1,4,2]$,

$p=[2,4,3,1]$, $q=[3,2,1,4]$,

$p=[2,4,3,1]$, $q=[4,1,2,3]$,

$p=[3,2,4,1]$, $q=[4,1,2,3]$,

$p=[3,4,1,2]$, $q=[4,1,2,3]$,

$p=[3,4,2,1]$, $q=[4,1,2,3]$,

$p=[3,4,2,1]$, $q=[4,1,3,2]$,

$p=[3,4,2,1]$, $q=[4,2,1,3]$.

Source: codeforces | Difficulty: VERY_HARD
Tags: Polynomials and generating functions

Combinatorics

Dynamic programming

Mathematics


Solutions:
import sys
input = sys.stdin.readline

def make_tree(n):
	i = 2
	while True:
		if i >= n * 2:
			tree = [0] * i
			break
		else:
			i *= 2
	return tree

def update(i, x):
	i += len(tree) // 2
	tree[i] = x
	i //= 2
	while True:
		if i == 0:
			break
		tree[i] = tree[2 * i] + tree[2 * i + 1]
		i //= 2
	return

def get_sum(s, t):
	s += len(tree) // 2
	t += len(tree) // 2
	ans = 0
	while True:
		if s > t:
			break
		if s % 2 == 0:
			s //= 2
		else:
			ans += tree[s]
			s = (s + 1) // 2
		if t % 2 == 1:
			t //= 2
		else:
			ans += tree[t]
			t = (t - 1) // 2
	return ans

def prime_factorize(n):
	ans = []
	for i in range(2, int(n ** (1 / 2)) + 1):
		while True:
			if n % i:
				break
			ans.append(i)
			n //= i
		if n == 1:
			break
	if not n == 1:
		ans.append(n)
	return ans

def euler_func(n):
	p = set(prime_factorize(n))
	x = n
	for i in p:
		x *= i - 1
		x //= i
	return x

def f(a, b):
	ans = 1
	for i in range(a + 1, b + 1):
		ans *= i
		ans %= mod
	return ans
(n, mod) = map(int, input().split())
x = euler_func(mod)
l = n + 5
fact = [1] * (l + 1)
for i in range(1, l + 1):
	fact[i] = i * fact[i - 1] % mod
if n <= 3:
	ans = 0
else:
	dp = [0] * 2000
	(dp[0], dp[1], dp[2], dp[3]) = (1, 2, 2, 1)
	tree = make_tree(2005)
	for i in range(4):
		update(i, dp[i])
	m = 3
	ans = 0
	for i in range(n - 3):
		dp0 = list(dp)
		ans0 = 0
		for j in range(i + 3):
			for k in range(1, m + 2):
				ans0 += get_sum(j + k + 1, 2004) * dp0[k - 1] % mod
				ans0 %= mod
				dp[j + k] += dp0[k - 1]
				dp[j + k] %= mod
				update(j + k, dp[j + k])
		m += i + 3
		ans += f(i + 4, n) * ans0 % mod
		ans %= mod
print(ans)


import sys
sys.stderr = sys.stdout

def egcd(a, b):
	r0 = a
	r1 = b
	s0 = 1
	s1 = 0
	t0 = 0
	t1 = 1
	while r1:
		q = r0 // r1
		(r0, r1) = (r1, r0 - q * r1)
		(s0, s1) = (s1, s0 - q * s1)
		(t0, t1) = (t1, t0 - q * t1)
	return (r0, s0)

def tf(n, M):
	C = [0] * (n + 1)
	C[n] = 1
	for k in reversed(range(n)):
		C[k] = C[k + 1] * (k + 1) % M
		if not C[k]:
			break
	return C

def app(n, M):
	K = [[[1]]]
	for i in range(n - 1):
		Ki = K[i]
		ni = len(Ki)
		Ki_ = [[0] * (i + 2) for _ in range(ni + i + 1)]
		K.append(Ki_)
		for (j, Kij) in enumerate(Ki):
			for (k, v) in enumerate(Kij):
				for h in range(i + 1):
					Ki_[j + h][k] += v
					Ki_[j + h][k] %= M
				Ki_[j + i + 1][i + 1] += v
				Ki_[j + i + 1][i + 1] %= M
	C = tf(n, M)
	s = 0
	for i in range(1, n):
		Ki = K[i]
		ni = len(Ki)
		mi = len(Ki[0])
		Si = [[None] * mi for _ in range(ni)]
		for j in range(ni):
			Si[j][mi - 1] = Ki[j][mi - 1]
			for k in reversed(range(mi - 1)):
				Si[j][k] = (Ki[j][k] + Si[j][k + 1]) % M
		for j in range(1, ni):
			for k in range(mi):
				Si[j][k] += Si[j - 1][k]
				Si[j][k] %= M
		si = 0
		for j in range(1, ni):
			for k in range(mi - 1):
				si += Ki[j][k] * Si[j - 1][k + 1] % M
				si %= M
		c = C[i + 1]
		si *= c
		si %= M
		s += si
		s %= M
	return s

def main():
	(n, m) = readinti()
	print(app(n, m))

def readint():
	return int(input())

def readinti():
	return map(int, input().split())

def readintt():
	return tuple(readinti())

def readintl():
	return list(readinti())

def readinttl(k):
	return [readintt() for _ in range(k)]

def readintll(k):
	return [readintl() for _ in range(k)]

def log(*args, **kwargs):
	print(*args, **kwargs, file=sys.__stderr__)
main()


import sys

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

def MI():
	return map(int, sys.stdin.readline().rstrip().split())

def LI():
	return list(map(int, sys.stdin.readline().rstrip().split()))

def LI2():
	return list(map(int, sys.stdin.readline().rstrip()))

def S():
	return sys.stdin.readline().rstrip()

def LS():
	return list(sys.stdin.readline().rstrip().split())

def LS2():
	return list(sys.stdin.readline().rstrip())
(n, mod) = MI()
m = n * (n - 1) // 2
dp = [[[0] * (n + 1) for _ in range(m + 1)] for _ in range(n + 1)]
dp[1][0][1] = 1
for i in range(1, n):
	for j in range(m + 1):
		for k in range(n + 1):
			a = dp[i][j][k]
			if not a:
				continue
			dp[i + 1][j + i][i + 1] += a
			dp[i + 1][j + i][i + 1] %= mod
			for l in range(i):
				dp[i + 1][j + i - 1 - l][k] += a
				dp[i + 1][j + i - 1 - l][k] %= mod
S_dp = [[[0] * (n + 1) for k in range(m + 1)] for i in range(n + 1)]
for i in range(n + 1):
	for k in range(n + 1):
		a = 0
		for j in range(m + 1):
			a += dp[i][j][k]
			a %= mod
			S_dp[i][j][k] = a
SS_dp = [[[0] * (n + 1) for j in range(m + 1)] for i in range(n + 1)]
for i in range(n + 1):
	for j in range(m + 1):
		a = 0
		for k in range(n + 1):
			a += S_dp[i][j][k]
			a %= mod
			SS_dp[i][j][k] = a
ANS = [0] * (n + 1)
for i in range(1, n + 1):
	ans = i * ANS[i - 1] % mod
	for k in range(n + 1):
		for j in range(1, m + 1):
			if not dp[i][j][k]:
				continue
			ans += dp[i][j][k] * (SS_dp[i][j - 1][-1] - SS_dp[i][j - 1][k])
			ans %= mod
	ANS[i] = ans
print(ANS[-1])


def divisors(M):
	d = []
	i = 1
	while M >= i ** 2:
		if M % i == 0:
			d.append(i)
			if i ** 2 != M:
				d.append(M // i)
		i = i + 1
	return d

def popcount(x):
	x = x - (x >> 1 & 1431655765)
	x = (x & 858993459) + (x >> 2 & 858993459)
	x = x + (x >> 4) & 252645135
	x = x + (x >> 8)
	x = x + (x >> 16)
	return x & 127

def eratosthenes(n):
	res = [0 for i in range(n + 1)]
	prime = set([])
	for i in range(2, n + 1):
		if not res[i]:
			prime.add(i)
			for j in range(1, n // i + 1):
				res[i * j] = 1
	return prime

def factorization(n):
	res = []
	for p in prime:
		if n % p == 0:
			while n % p == 0:
				n //= p
			res.append(p)
	if n != 1:
		res.append(n)
	return res

def euler_phi(n):
	res = n
	for x in range(2, n + 1):
		if x ** 2 > n:
			break
		if n % x == 0:
			res = res // x * (x - 1)
			while n % x == 0:
				n //= x
	if n != 1:
		res = res // n * (n - 1)
	return res

def ind(b, n):
	res = 0
	while n % b == 0:
		res += 1
		n //= b
	return res

def isPrimeMR(n):
	d = n - 1
	d = d // (d & -d)
	L = [2, 3, 5, 7, 11, 13, 17]
	for a in L:
		t = d
		y = pow(a, t, n)
		if y == 1:
			continue
		while y != n - 1:
			y = y * y % n
			if y == 1 or t == n - 1:
				return 0
			t <<= 1
	return 1

def findFactorRho(n):
	from math import gcd
	m = 1 << n.bit_length() // 8
	for c in range(1, 99):
		f = lambda x: (x * x + c) % n
		(y, r, q, g) = (2, 1, 1, 1)
		while g == 1:
			x = y
			for i in range(r):
				y = f(y)
			k = 0
			while k < r and g == 1:
				ys = y
				for i in range(min(m, r - k)):
					y = f(y)
					q = q * abs(x - y) % n
				g = gcd(q, n)
				k += m
			r <<= 1
		if g == n:
			g = 1
			while g == 1:
				ys = f(ys)
				g = gcd(abs(x - ys), n)
		if g < n:
			if isPrimeMR(g):
				return g
			elif isPrimeMR(n // g):
				return n // g
			return findFactorRho(g)

def primeFactor(n):
	i = 2
	ret = {}
	rhoFlg = 0
	while i * i <= n:
		k = 0
		while n % i == 0:
			n //= i
			k += 1
		if k:
			ret[i] = k
		i += 1 + i % 2
		if i == 101 and n >= 2 ** 20:
			while n > 1:
				if isPrimeMR(n):
					(ret[n], n) = (1, 1)
				else:
					rhoFlg = 1
					j = findFactorRho(n)
					k = 0
					while n % j == 0:
						n //= j
						k += 1
					ret[j] = k
	if n > 1:
		ret[n] = 1
	if rhoFlg:
		ret = {x: ret[x] for x in sorted(ret)}
	return ret

def divisors(n):
	res = [1]
	prime = primeFactor(n)
	for p in prime:
		newres = []
		for d in res:
			for j in range(prime[p] + 1):
				newres.append(d * p ** j)
		res = newres
	res.sort()
	return res

def xorfactorial(num):
	if num == 0:
		return 0
	elif num == 1:
		return 1
	elif num == 2:
		return 3
	elif num == 3:
		return 0
	else:
		x = baseorder(num)
		return 2 ** x * ((num - 2 ** x + 1) % 2) + function(num - 2 ** x)

def xorconv(n, X, Y):
	if n == 0:
		res = [X[0] * Y[0] % mod]
		return res
	x = [X[i] + X[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	y = [Y[i] + Y[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	z = [X[i] - X[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	w = [Y[i] - Y[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	res1 = xorconv(n - 1, x, y)
	res2 = xorconv(n - 1, z, w)
	former = [(res1[i] + res2[i]) * inv for i in range(2 ** (n - 1))]
	latter = [(res1[i] - res2[i]) * inv for i in range(2 ** (n - 1))]
	former = list(map(lambda x: x % mod, former))
	latter = list(map(lambda x: x % mod, latter))
	return former + latter

def merge_sort(A, B):
	(pos_A, pos_B) = (0, 0)
	(n, m) = (len(A), len(B))
	res = []
	while pos_A < n and pos_B < m:
		(a, b) = (A[pos_A], B[pos_B])
		if a < b:
			res.append(a)
			pos_A += 1
		else:
			res.append(b)
			pos_B += 1
	res += A[pos_A:]
	res += B[pos_B:]
	return res

class UnionFindVerSize:

	def __init__(self, N):
		self._parent = [n for n in range(0, N)]
		self._size = [1] * N
		self.group = N

	def find_root(self, x):
		if self._parent[x] == x:
			return x
		self._parent[x] = self.find_root(self._parent[x])
		stack = [x]
		while self._parent[stack[-1]] != stack[-1]:
			stack.append(self._parent[stack[-1]])
		for v in stack:
			self._parent[v] = stack[-1]
		return self._parent[x]

	def unite(self, x, y):
		gx = self.find_root(x)
		gy = self.find_root(y)
		if gx == gy:
			return
		self.group -= 1
		if self._size[gx] < self._size[gy]:
			self._parent[gx] = gy
			self._size[gy] += self._size[gx]
		else:
			self._parent[gy] = gx
			self._size[gx] += self._size[gy]

	def get_size(self, x):
		return self._size[self.find_root(x)]

	def is_same_group(self, x, y):
		return self.find_root(x) == self.find_root(y)

class WeightedUnionFind:

	def __init__(self, N):
		self.parent = [i for i in range(N)]
		self.size = [1 for i in range(N)]
		self.val = [0 for i in range(N)]
		self.flag = True
		self.edge = [[] for i in range(N)]

	def dfs(self, v, pv):
		stack = [(v, pv)]
		new_parent = self.parent[pv]
		while stack:
			(v, pv) = stack.pop()
			self.parent[v] = new_parent
			for (nv, w) in self.edge[v]:
				if nv != pv:
					self.val[nv] = self.val[v] + w
					stack.append((nv, v))

	def unite(self, x, y, w):
		if not self.flag:
			return
		if self.parent[x] == self.parent[y]:
			self.flag = self.val[x] - self.val[y] == w
			return
		if self.size[self.parent[x]] > self.size[self.parent[y]]:
			self.edge[x].append((y, -w))
			self.edge[y].append((x, w))
			self.size[x] += self.size[y]
			self.val[y] = self.val[x] - w
			self.dfs(y, x)
		else:
			self.edge[x].append((y, -w))
			self.edge[y].append((x, w))
			self.size[y] += self.size[x]
			self.val[x] = self.val[y] + w
			self.dfs(x, y)

class Dijkstra:

	class Edge:

		def __init__(self, _to, _cost):
			self.to = _to
			self.cost = _cost

	def __init__(self, V):
		self.G = [[] for i in range(V)]
		self._E = 0
		self._V = V

	@property
	def E(self):
		return self._E

	@property
	def V(self):
		return self._V

	def add_edge(self, _from, _to, _cost):
		self.G[_from].append(self.Edge(_to, _cost))
		self._E += 1

	def shortest_path(self, s):
		import heapq
		que = []
		d = [10 ** 15] * self.V
		d[s] = 0
		heapq.heappush(que, (0, s))
		while len(que) != 0:
			(cost, v) = heapq.heappop(que)
			if d[v] < cost:
				continue
			for i in range(len(self.G[v])):
				e = self.G[v][i]
				if d[e.to] > d[v] + e.cost:
					d[e.to] = d[v] + e.cost
					heapq.heappush(que, (d[e.to], e.to))
		return d

def Z_algorithm(s):
	N = len(s)
	Z_alg = [0] * N
	Z_alg[0] = N
	i = 1
	j = 0
	while i < N:
		while i + j < N and s[j] == s[i + j]:
			j += 1
		Z_alg[i] = j
		if j == 0:
			i += 1
			continue
		k = 1
		while i + k < N and k + Z_alg[k] < j:
			Z_alg[i + k] = Z_alg[k]
			k += 1
		i += k
		j -= k
	return Z_alg

class BIT:

	def __init__(self, n, mod=0):
		self.BIT = [0] * (n + 1)
		self.num = n
		self.mod = mod

	def query(self, idx):
		res_sum = 0
		mod = self.mod
		while idx > 0:
			res_sum += self.BIT[idx]
			if mod:
				res_sum %= mod
			idx -= idx & -idx
		return res_sum

	def update(self, idx, x):
		mod = self.mod
		while idx <= self.num:
			self.BIT[idx] += x
			if mod:
				self.BIT[idx] %= mod
			idx += idx & -idx
		return

class dancinglink:

	def __init__(self, n, debug=False):
		self.n = n
		self.debug = debug
		self._left = [i - 1 for i in range(n)]
		self._right = [i + 1 for i in range(n)]
		self.exist = [True for i in range(n)]

	def pop(self, k):
		if self.debug:
			assert self.exist[k]
		L = self._left[k]
		R = self._right[k]
		if L != -1:
			if R != self.n:
				(self._right[L], self._left[R]) = (R, L)
			else:
				self._right[L] = self.n
		elif R != self.n:
			self._left[R] = -1
		self.exist[k] = False

	def left(self, idx, k=1):
		if self.debug:
			assert self.exist[idx]
		res = idx
		while k:
			res = self._left[res]
			if res == -1:
				break
			k -= 1
		return res

	def right(self, idx, k=1):
		if self.debug:
			assert self.exist[idx]
		res = idx
		while k:
			res = self._right[res]
			if res == self.n:
				break
			k -= 1
		return res

class SparseTable:

	def __init__(self, A, merge_func, ide_ele):
		N = len(A)
		n = N.bit_length()
		self.table = [[ide_ele for i in range(n)] for i in range(N)]
		self.merge_func = merge_func
		for i in range(N):
			self.table[i][0] = A[i]
		for j in range(1, n):
			for i in range(0, N - 2 ** j + 1):
				f = self.table[i][j - 1]
				s = self.table[i + 2 ** (j - 1)][j - 1]
				self.table[i][j] = self.merge_func(f, s)

	def query(self, s, t):
		b = t - s + 1
		m = b.bit_length() - 1
		return self.merge_func(self.table[s][m], self.table[t - 2 ** m + 1][m])

class BinaryTrie:

	class node:

		def __init__(self, val):
			self.left = None
			self.right = None
			self.max = val

	def __init__(self):
		self.root = self.node(-10 ** 15)

	def append(self, key, val):
		pos = self.root
		for i in range(29, -1, -1):
			pos.max = max(pos.max, val)
			if key >> i & 1:
				if pos.right is None:
					pos.right = self.node(val)
					pos = pos.right
				else:
					pos = pos.right
			elif pos.left is None:
				pos.left = self.node(val)
				pos = pos.left
			else:
				pos = pos.left
		pos.max = max(pos.max, val)

	def search(self, M, xor):
		res = -10 ** 15
		pos = self.root
		for i in range(29, -1, -1):
			if pos is None:
				break
			if M >> i & 1:
				if xor >> i & 1:
					if pos.right:
						res = max(res, pos.right.max)
					pos = pos.left
				else:
					if pos.left:
						res = max(res, pos.left.max)
					pos = pos.right
			elif xor >> i & 1:
				pos = pos.right
			else:
				pos = pos.left
		if pos:
			res = max(res, pos.max)
		return res

def solveequation(edge, ans, n, m):
	x = [0] * m
	used = [False] * n
	for v in range(n):
		if used[v]:
			continue
		y = dfs(v)
		if y != 0:
			return False
	return x

	def dfs(v):
		used[v] = True
		r = ans[v]
		for (to, dire, id) in edge[v]:
			if used[to]:
				continue
			y = dfs(to)
			if dire == -1:
				x[id] = y
			else:
				x[id] = -y
			r += y
		return r

class SegmentTree:

	def __init__(self, init_val, segfunc, ide_ele):
		n = len(init_val)
		self.segfunc = segfunc
		self.ide_ele = ide_ele
		self.num = 1 << (n - 1).bit_length()
		self.tree = [ide_ele] * 2 * self.num
		self.size = n
		for i in range(n):
			self.tree[self.num + i] = init_val[i]
		for i in range(self.num - 1, 0, -1):
			self.tree[i] = self.segfunc(self.tree[2 * i], self.tree[2 * i + 1])

	def update(self, k, x):
		k += self.num
		self.tree[k] = x
		while k > 1:
			self.tree[k >> 1] = self.segfunc(self.tree[k], self.tree[k ^ 1])
			k >>= 1

	def query(self, l, r):
		if r == self.size:
			r = self.num
		res = self.ide_ele
		l += self.num
		r += self.num
		while l < r:
			if l & 1:
				res = self.segfunc(res, self.tree[l])
				l += 1
			if r & 1:
				res = self.segfunc(res, self.tree[r - 1])
			l >>= 1
			r >>= 1
		return res

	def bisect_l(self, l, r, x):
		l += self.num
		r += self.num
		Lmin = -1
		Rmin = -1
		while l < r:
			if l & 1:
				if self.tree[l] <= x and Lmin == -1:
					Lmin = l
				l += 1
			if r & 1:
				if self.tree[r - 1] <= x:
					Rmin = r - 1
			l >>= 1
			r >>= 1
		if Lmin != -1:
			pos = Lmin
			while pos < self.num:
				if self.tree[2 * pos] <= x:
					pos = 2 * pos
				else:
					pos = 2 * pos + 1
			return pos - self.num
		elif Rmin != -1:
			pos = Rmin
			while pos < self.num:
				if self.tree[2 * pos] <= x:
					pos = 2 * pos
				else:
					pos = 2 * pos + 1
			return pos - self.num
		else:
			return -1
import sys, random, bisect
from collections import deque, defaultdict
from heapq import heapify, heappop, heappush
from itertools import permutations
from math import gcd, log
input = lambda : sys.stdin.readline().rstrip()
mi = lambda : map(int, input().split())
li = lambda : list(mi())
(N, mod) = mi()
g = [1 for i in range(N + 1)]
for i in range(1, N + 1):
	g[i] = g[i - 1] * i % mod
COMB = [[0 for j in range(N + 1)] for i in range(N + 1)]
COMB[0][0] = 1
for i in range(1, N + 1):
	COMB[i][0] = 1
	for j in range(1, i + 1):
		COMB[i][j] = COMB[i - 1][j] + COMB[i - 1][j - 1]
		COMB[i][j] %= mod
dp = [[0 for j in range(N * (N - 1) // 2 + 1)] for i in range(N + 1)]
dp[1][0] = 1
for i in range(2, N + 1):
	for k in range(N * (N - 1) // 2 + 1):
		if not dp[i - 1][k]:
			continue
		dp[i][k] += dp[i - 1][k]
		dp[i][k] %= mod
		if k + i <= N * (N - 1) // 2:
			dp[i][k + i] -= dp[i - 1][k]
			dp[i][k + i] %= mod
	for k in range(1, N * (N - 1) // 2 + 1):
		dp[i][k] += dp[i][k - 1]
		dp[i][k] %= mod
imos = [[dp[i][j] for j in range(N * (N - 1) // 2 + 1)] for i in range(N + 1)]
for i in range(1, N + 1):
	for j in range(1, N * (N - 1) // 2 + 1):
		imos[i][j] += imos[i][j - 1]
		imos[i][j] %= mod
res = 0
for i in range(1, N):
	cnt = N - i
	inverse = [0 for inv in range(-(N - i + 1), 0)]
	for j in range(1, N - i + 1):
		for k in range(j + 1, N - i + 2):
			inv = j - k
			inverse[inv] += 1
			inverse[inv] %= mod
	for inv in range(-(N - i + 1), 0):
		for ip in range(cnt * (cnt - 1) // 2 + 1):
			if ip + inv <= 0:
				continue
			R = ip + inv - 1
			res += COMB[N][i - 1] * dp[cnt][ip] % mod * (imos[cnt][R] * inverse[inv] % mod) * g[i - 1] % mod
			res %= mod
print(res)

def inversion(p):
	res = 0
	for i in range(N):
		for j in range(i + 1, N):
			if p[i] > p[j]:
				res += 1
	return res

def solve(N):
	L = list(permutations([i + 1 for i in range(N)]))
	res = 0
	check = [0 for i in range(N + 1)]
	for p in L:
		for q in L:
			if p < q and inversion(p) > inversion(q):
				res += 1
				for i in range(N):
					if p[i] < q[i]:
						check[i + 1] += 1
						break
	for i in range(1, N + 1):
		print(i, check[i])
	return res


def solve():
	(n, mod) = map(int, input().split())
	m = n * (n - 1) // 2
	dp = [[0] * (n * (n - 1) + 1) for i in range(n + 1)]
	dp1 = [[0] * (n * (n - 1) + 1) for i in range(n + 1)]
	dp[1][m] = 1
	for i in range(2, n + 1):
		for x in range(m - (i - 1) * (i - 2) // 2, m + (i - 1) * (i - 2) // 2 + 1):
			for d in range(-i + 1, i):
				dp[i][x + d] = (dp[i][x + d] + dp[i - 1][x] * (i - abs(d))) % mod
	for i in range(2, n + 1):
		for x in range(m - (i - 1) * (i - 2) // 2, m + (i - 1) * (i - 2) // 2 + 1):
			dp1[i][x] = dp1[i - 1][x] * i % mod
		for x in range(m - (i - 1) * (i - 2) // 2, m + (i - 1) * (i - 2) // 2 + 1):
			for d in range(1, i):
				dp1[i][x - d] = (dp1[i][x - d] + dp[i - 1][x] * (i - d)) % mod
	return sum(dp1[-1][m + 1:]) % mod
import sys
input = lambda : sys.stdin.readline().rstrip()
print(solve())


def solve(n, mod):
	if n <= 3:
		return 0
	else:
		last_ans = 17
		ans = 17
		size = (n - 1) * (n - 2) // 2 + 1
		inv = 1
		record = [[0 for i in range(size + 1)] for i in range(2)]
		record[0][0:4] = [1, 2, 2, 1]
		for m in range(4, n):
			tp_size = m * (m - 1) // 2
			ans = (m + 1) * last_ans
			tp = 0
			for i in range(m):
				tp += record[inv ^ 1][i]
				record[inv][i] = tp % mod
			for i in range(m, tp_size // 2 + 1):
				tp += record[inv ^ 1][i]
				tp -= record[inv ^ 1][i - m]
				record[inv][i] = tp % mod
			for i in range(tp_size // 2 + 1):
				record[inv][tp_size - i] = record[inv][i]
			total = 0
			for i in range(tp_size - 1):
				tp = 0
				for j in range(max(1, -tp_size + i + m + 2), m + 1):
					tp += j * int(record[inv][m + 2 - j + i])
				total += int(record[inv][i])
				total %= mod
				tp *= total
				tp %= mod
				ans += tp
				ans %= mod
			last_ans = ans
			inv ^= 1
		return int(ans) % mod

def main():
	(n, mod) = map(int, input().split())
	print(solve(n, mod))
main()


from functools import lru_cache
import itertools

def invnum(seq):
	N = len(seq)
	res = 0
	for i in range(N):
		for j in range(i + 1, N):
			if seq[i] > seq[j]:
				res += 1
	return res

def naive(N, k, flag):
	res = 0
	for p in itertools.permutations(range(N), N):
		for q in itertools.permutations(range(N), N):
			if (p < q or not flag) and invnum(p) - invnum(q) > k:
				res += 1
	return res
(n, mod) = map(int, input().split())
memo = [[{}, {}] for _ in range(n + 1)]

def solve(N, k, flag):
	if k in memo[N][flag]:
		return memo[N][flag][k]
	if N == 1:
		res = -1
		if flag:
			res = 0
		elif k < 0:
			res = 1
		else:
			res = 0
		memo[N][flag][k] = res
		return res
	invnum_INF = N * (N - 1) // 2
	if k < -invnum_INF - 1:
		res = solve(N, -invnum_INF - 1, flag)
		memo[N][flag][k] = res
		return res
	if k >= invnum_INF:
		memo[N][flag][k] = 0
		return 0
	res = 0
	for delta in range(-N + 1, N):
		if flag and delta > 0:
			continue
		wt = min(N - delta, delta + N)
		res += solve(N - 1, k - delta, (delta == 0) & flag) * wt
		res %= mod
	memo[N][flag][k] = res % mod
	return res % mod
ans = solve(n, 0, 1) % mod
print(ans)