If we can count the number of \(A\)’s with all elements nonzero, the answer to the original problem follows easily via DP, so we focus on counting these.

Viewing the operations in reverse, the sequences satisfying the condition are those where we can repeatedly perform “choose \(l,r\) such that \(A_i\in \{0,r(r-1)/2+l\}\) for all \(i=l,\ldots,r\), and set \(A_i=0\) (we call this ‘erasing \(r(r-1)/2+l\)’)” to eventually make all elements \(0\).

In this reversed sequence of operations, there exists a partial order “\((l_1,r_1)\) must be done before \((l_2,r_2)\),” and this partial order forms a DAG. Let \(S\) be the set of values that can be erased as the last step (vertices with out-degree \(0\) in the DAG). We want to count the number of valid \(A\)’s for each \(S\), but this is difficult to do directly.

Instead, we consider adding \((-1)^{|T|-1}\) for every nonempty subset \(T\) of \(S\).

For a fixed \(T\), we can count easily. Listing the values in \(T\), we have non-overlapping intervals arranged side by side. For the parts not in \(T\), everything can be erased within those small intervals.

This gives us a DP. Let \(\mathrm{dp}1[i]\) be “the number of sequences \(A\) of length \(i\) that can be constructed with all elements nonzero,” and \(\mathrm{dp}2[i]\) be “the sum of \((-1)^{|T|-1}\) over all sequences \(A\) of length \(i\) that can be constructed with all elements nonzero where position \(i\) contains a value in \(T\).”

The most complex transition, from \(\mathrm{dp}2[l]\) to \(\mathrm{dp}2[r]\), is:

• If the values at \(l\) and \(r\) are the same: since the part between \(l\) and \(r\) can be erased, multiply by \(\mathrm{dp}1[r-1-l]\).
• If the values at \(l\) and \(r\) are different: since the right endpoint of \(l\)’s value and the left endpoint of \(r\)’s value lie in \([l,r]\), multiply by \((r-l)^2\); also multiply by \(\mathrm{dp}1[r-1-l]\) since the part between \(l\) and \(r\) can be erased; and multiply by \(-1\) since the number of distinct values increases.

The time complexity is \(O(N^2)\).

n, mod = map(int, input().split())
dp1, dp2, ans = [0] * (n + 1), [0] * (n + 1), [0] * (n + 1)
dp1[0] = 1
ans[0] = 1
for i in range(1, n + 1):
    dp2[i] = dp1[i - 1] * i
    for j in range(i):
        dp1[i] += dp2[j + 1] * (i - j) % mod * dp1[i - j - 1] % mod
        dp2[i] += dp2[j] * (1 - (i - j) * (i - j)) % mod * dp1[i - j - 1] % mod
        ans[i] += ans[j] * dp1[i - j - 1] % mod
    ans[i] += dp1[i]
    dp1[i] %= mod
    dp2[i] %= mod
    ans[i] %= mod
print(ans[n])