Problem:
 Score : 1200 points    Problem Statement Given is a positive integer N .
Find the number of permutations (P_1,P_2,\cdots,P_{3N}) of (1,2,\cdots,3N) that can be generated through the procedure below.
This number can be enormous, so print it modulo a prime number M .   Make N sequences A_1,A_2,\cdots,A_N of length 3 each, using each of the integers 1 through 3N exactly once.  Let P be an empty sequence, and do the following operation 3N times.  Among the elements that are at the beginning of one of the sequences A_i that is non-empty, let the smallest be x .  Remove x from the sequence, and add x at the end of P .         Constraints  1 \leq N \leq 2000  10^8 \leq M \leq 10^9+7  M is a prime number.  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M       Output Print the number of permutations modulo M .        Sample Input 1 1 998244353      Sample Output 1 6  All permutations of length 3 count.       Sample Input 2 2 998244353      Sample Output 2 261       Sample Input 3 314 1000000007      Sample Output 3 182908545

Solution (Python):
N,M=map(int,input().split())

dp=[[0 for i in range(2*N+1)] for j in range(3*N+1)]

for j in range(N+1):
    dp[0][j]=1

for i in range(1,3*N+1):
    for j in range(2*N+1):
        if 2*N-1>=j>=1:
            dp[i][j]=(dp[i-1][j-1]+dp[i-2][j+1]*(i-1)+dp[i-3][j]*(i-1)*(i-2))%M
        elif j==0:
            dp[i][j]=(dp[i-1][j]+dp[i-2][j+1]*(i-1)+dp[i-3][j]*(i-1)*(i-2))%M
        else:
            dp[i][j]=(dp[i-1][j-1]+dp[i-3][j]*(i-1)*(i-2))%M

print(dp[3*N][N]%M)

Solution (Python):
N,m=map(int,input().split())
M=3*N
DP=[[0]*(N+5) for i in range(M+5)]
DP[0][0]=1
for i in range(1,M+1):
    for j in range(N+1):
        DP[i][j]=(DP[i-1][j]+(M-i+1)*DP[i-2][j-1]+(M-i+1)*(M-i+2)*DP[i-3][j-1])%m
print(sum(DP[3*N][:N+1])%m)

Solution (Python):
N,mod=map(int,input().split())

DP=[[0]*(N+5) for i in range(3*N+5)]
DP[0][0]=1

for i in range(3*N):
    for j in range(N+1):
        DP[i][j]%=mod
        DP[i+1][j]+=DP[i][j]
        DP[i+2][j+1]+=((3*N-i)-1)*DP[i][j]
        DP[i+3][j+1]+=((3*N-i)-1)*((3*N-i)-2)*DP[i][j]

ANS=0
for i in range(N+1):
    ANS+=DP[3*N][i]
print(ANS%mod)

Solution (Python):
import sys
import numpy as np

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

N, MOD = map(int, read().split())

N *= 3


dp = np.zeros((N + 10, N + N + 1), np.int64)
dp[0, N] = 1
for n in range(N + 1):
    dp[n] %= MOD
    dp[n + 1, 1:] += dp[n, :-1]
    dp[n + 2, :-1] += dp[n, 1:] * (n + 1)
    dp[n + 3] += dp[n] * ((n + 1) * (n + 2))

x = dp[N, N:].sum() % MOD
print(x)

Solution (Python):
import sys
readline = sys.stdin.readline
N, MOD = map(int, readline().split())

dp = [[0]*(N+5) for _ in range(3*N+4)]
dp[0][0] = 1

for i in range(3*N):
    for j in range(N+1):
        d = dp[i][j]
        dp[i+1][j] = (dp[i+1][j]+d)%MOD
        dp[i+2][j+1] = (dp[i+2][j+1]+(i+1)*d)%MOD
        dp[i+3][j+1] = (dp[i+3][j+1]+(i+2)*(i+1)*d)%MOD

res = 0
for j in range(N+1):
    res = (res+dp[3*N][j])%MOD
print(res)

Solution (Python):
from numpy import*
n,M=map(int,input().split())
l=n*3+1
d=zeros((l,n*5),int64)
d[0][0]=1
for i in range(1,l):j,k=i-1,i-2;d[i]=((d[i-3]*k+roll(d[k],-1))*j+roll(d[j],1))%M
print(sum(d[-1][:l])%M)

Solution (Python):
from numpy import*;n,M=map(int,input().split())
l=n*3+1
d=zeros((l,n*5),int64)
d[0][0]=1
for i in range(1,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+roll(d[k],-1)*j+roll(d[j],1))%M
print(sum(d[-1][:l])%M)

Solution (Python):
from numpy import *
n,M=map(int,input().split())
l=n*3+1
d=zeros((l,n*5),int64)
d[0][0]=1
for i in range(1,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+roll(d[k],-1)*j+roll(d[j],1))%M
print(sum(d[-1][:l])%M)

Solution (Python):
from numpy import *
n,M=map(int,input().split())
l=n*3+1
d=zeros((l,n*5),int64)
d[0][0]=1
for i in range(1,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+roll(d[k],-1)*j+roll(d[j],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
from numpy import *
n,M=map(int,input().split())
l=n*3+1
r=roll
d=zeros((l,n*5),int64)
d[0][0]=1
for i in range(1,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+r(d[k],-1)*j+r(d[j],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
import numpy as p
n,M=map(int,input().split())
l=n*3+1
r=p.roll
d=p.zeros((l,n*5),p.int64)
d[0][0]=1
for i in range(1,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+r(d[k],-1)*j+r(d[j],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
import numpy as p
n,M=map(int,input().split())
l=n*3+1
r=p.roll
d=p.zeros((l,n*5),p.int64)
d[0][0]=d[1][1]=d[2][2]=d[2][-1]=1
for i in range(3,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+r(d[k],-1)*j+r(d[j],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
import numpy as p
n,M=map(int,input().split())
l=n*3+1
d=p.zeros((l,n*5),p.int64)
d[0][0]=d[1][1]=d[2][2]=d[2][-1]=1
for i in range(3,l):j,k=i-1,i-2;d[i]=(d[i-3]*k*j+p.roll(d[k],-1)*j+p.roll(d[j],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
import numpy as p
n,M=map(int,input().split())
l=n*3+1
d=p.zeros((l,n*5),p.int64)
d[0][0]=d[1][1]=d[2][2]=d[2][-1]=1
for i in range(3,l):d[i]=(d[i-3]*(i-2)*(i-1)+p.roll(d[i-2],-1)*(i-1)+p.roll(d[i-1],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
import numpy as p
n,M=map(int,input().split())
l=n*3+1
d=p.zeros((l,n*5),p.int64)
d[0][0]=d[1][1]=d[2][2]=d[2][-1]=1
for i in range(3,l):
  d[i]=(d[i-3]*(i-2)*(i-1)+p.roll(d[i-2],-1)*(i-1)+p.roll(d[i-1],1))%M
print(d[-1][:l].sum()%M)

Solution (Python):
n,m=map(int,input().split())
N=3*n
d=[[0]*(N+1) for i in range(N+4)]
d[0][0]=1
for i in range(N):
  a,b,c,e=d[i:i+4]
  for j in range(N):J=j+1;b[j]=(b[j]+a[j])%m;c[J]=(c[J]+a[j]*(i+1))%m;e[J]=(e[J]+a[j]*(i+1)*(i+2))%m
print(sum(d[N][:n+1])%m)

Solution (Python):
n,m=map(int,input().split())
N=3*n
d=[[0]*(N+1) for i in range(N+4)]
d[0][0]=1
for i in range(N):
  a,b,c,e=d[i:i+4]
  for j in range(N):J=j+1;b[j]=(b[j]+a[j])%m;c[J]=(c[J]+a[j]*(i+1))%m;e[J]=(e[J]+a[j]*(i+1)*(i+2))%m
print(sum(d[N][:n+1])%m)

Solution (Python):
n,m=map(int,input().split())
N=3*n
d=[[0]*(N+1) for i in range(N+5)]
d[0][0]=1
for i in range(N):
  a,b,c,e,*_=d[i:]
  for j in range(N):
    J=j+1
    b[j]=(b[j]+a[j])%m
    c[J]=(c[J]+a[j]*(i+1))%m
    e[J]=(e[J]+a[j]*(i+1)*(i+2))%m
print(sum(d[N][:n+1])%m)

Solution (Python):
n,m=map(int,input().split())
N=3*n
d=[[0]*(N+1) for i in range(N+4)]
d[0][0]=1
for i in range(N):
  a,b,c,e=d[i:i+4]
  for j in range(N):
    J=j+1
    b[j]=(b[j]+a[j])%m
    c[J]=(c[J]+a[j]*(i+1))%m
    e[J]=(e[J]+a[j]*(i+1)*(i+2))%m
print(sum(d[N][:n+1])%m)

Solution (Python):
n,m=map(int,input().split())
n3=3*n
x=lambda a,b:(a+b)%mod
d=[[0]*(n+2) for i in range(n3+4)]
d[0][0]=1
for i in range(n3):
  d0,d1,d2,d3=d[i:i+4]
  for j in range(n+1):
    d1[j]=(d1[j]+d0[j])%m
    d2[j+1]=(d2[j+1]+d0[j]*(i+1))%m
    d3[j+1]=(d3[j+1]+d0[j]*(i+1)*(i+2))%m
print(sum(d[n3][:n+1])%m)

Solution (Python):
import numpy as np
n,MOD=map(int,input().split())
n3=n*3+1
dp=np.zeros((n3,n*5),np.int64)
dp[0][0]=dp[1][1]=dp[2][2]=dp[2][-1]=1
for i in range(3,n3):
  dp[i]=(dp[i-3]*(i-2)*(i-1)+np.roll(dp[i-2],-1)*(i-1)+np.roll(dp[i-1],1))%MOD
print(dp[-1][:n3].sum()%MOD)

Solution (Python):
import numpy as np

n, MOD = map(int, input().split())
n3 = n * 3
dp = np.zeros((n3 + 1, int(n * 4.5) + 1), np.int64)
dp[0][0] = dp[1][1] = dp[2][2] = dp[2][-1] = 1

for i in range(3, n3 + 1):
    dp[i] += dp[i - 3] * (i - 2) * (i - 1)
    dp[i] += np.roll(dp[i - 2], -1) * (i - 1)
    dp[i] += np.roll(dp[i - 1], 1)
    dp[i] %= MOD

print(dp[-1][:n3 + 1].sum() % MOD)

Solution (Python):
n, MOD = map(int, input().split())

n3 = n * 3
dp = [[0] * int(i * 1.5 + 5) for i in range(n3 + 1)]
dp[0][0] = dp[1][1] = dp[2][2] = dp[2][-1] = 1

for i in range(3, n3 + 1):
    i12 = (i - 1) * (i - 2)
    dpi3, dpi2, dpi1, dpi = dp[i - 3:i + 1]
    for j in range(-(i // 2), i + 1):
        tmp = dpi3[j] * i12
        tmp += dpi2[j + 1] * (i - 1)
        tmp += dpi1[j - 1]
        dpi[j] = tmp % MOD

print(sum(dp[-1][:n3 + 1]) % MOD)

Solution (Python):
from functools import lru_cache

N, M = map(int, input().split())

@lru_cache(maxsize=None)
def mod_inv(x):
  x1, y1, z1 = 1, 0, x
  x2, y2, z2 = 0, 1, M
  while z1 != 1:
    d, m = divmod(z2, z1)
    x1, x2 = x2-d*x1, x1
    y1, y2 = y2-d*y1, y1
    z1, z2 = m, z1
  return x1%M

def gen_Y(i):
  
  s = 3*(N-i)
  r = s*(s-1)>>1
  d_r = (s<<1)-3
  for j in range(1, N-i+1):
    yield r * mod_inv(j)
    r -= d_r
    d_r -= 4

def gen_X():
  
  a = N
  b = 3*N - 1
  for i in range(1, N+1):
    yield a * b * (b-1) * mod_inv(i)
    a -= 1
    b -= 3

def acc_mulmod(it, init=1):
  x = init
  yield x
  for y in it:
    x = x*y % M
    yield x

ans = sum(sum(acc_mulmod(gen_Y(i), A)) for i, A in enumerate(acc_mulmod(gen_X())))%M

print(ans)

Solution (Python):
from functools import lru_cache

N, M = map(int, input().split())

@lru_cache(maxsize=None)
def mod_inv(x):
  x1, y1, z1 = 1, 0, x
  x2, y2, z2 = 0, 1, M
  while z1 != 1:
    d, m = divmod(z2, z1)
    x1, x2 = x2-d*x1, x1
    y1, y2 = y2-d*y1, y1
    z1, z2 = m, z1
  return x1%M

def gen_Y(i):
  
  s = 3*(N-i)
  r = s*(s-1)>>1
  d_r = (s<<1)-3
  for j in range(1, N-i+1):
    yield r * mod_inv(j)
    r -= d_r
    d_r -= 4

def gen_X():
  
  a = N
  b = 3*N - 1
  for i in range(1, N+1):
    yield a * b * (b-1) * mod_inv(i)
    a -= 1
    b -= 3

def acc_mulmod(it, init=1):
  x = init
  yield x
  for y in it:
    x = x*y % M
    yield x

ans = sum(sum(acc_mulmod(gen_Y(i), A)) for i, A in enumerate(acc_mulmod(gen_X())))%M

print(ans)

Solution (Python):
from functools import lru_cache, reduce
from itertools import accumulate

N, M = map(int, input().split())

@lru_cache(maxsize=None)
def mod_inv(x):
  x1, y1, z1 = 1, 0, x
  x2, y2, z2 = 0, 1, M
  while z1 != 1:
    d, m = divmod(z2, z1)
    x1, x2 = x2-d*x1, x1
    y1, y2 = y2-d*y1, y1
    z1, z2 = m, z1
  return x1%M

def gen_Y(i, A):
  yield A
  
  s = 3*(N-i)
  r = s*(s-1)>>1
  d_r = (s<<1)-3
  for j in range(1, N-i+1):
    yield r * mod_inv(j)
    r -= d_r
    d_r -= 4

def gen_X():
  
  yield 1
  a = N
  b = 3*N - 1
  for i in range(1, N+1):
    yield a * b * (b-1) * mod_inv(i)
    a -= 1
    b -= 3

def mul_mod(x, y):
  return x * y % M

def acc_mod(it):
  return accumulate(it, mul_mod)

ans = sum(sum(acc_mod(gen_Y(i, A))) for i, A in enumerate(acc_mod(gen_X())))%M

print(ans)

Solution (Python):
N, P = map(int, input().split())
K = 96
m = int(("1" * 32 + "0" * 64) * 8192, 2)
pa = (1 << 64) - ((1 << 64) % P)
modP = lambda x: x - ((x & m) >> 64) * pa

def getsum(x):
    for i in range(13):
        x += x >> (K << i)
    return x & ((1 << K) - 1)

a, b, c = 1 << K * N, 0, 0
fa = 1
for i in range(1, 3 * N + 1):
    a, b, c = modP(((a << K) + (b >> K) + c) * pow(i, P-2, P)), a, b
    fa = fa * i % P

print(getsum(a >> K * N) * fa % P)

Solution (Python):
N, M = map(int, input().split())
dp = [[0]*(6*N+1) for _ in range(3*N+1)]
dp[0][3*N] = 1
for i in range(3*N):
    for j in range(-3*N+1, 3*N):
        if i+1 <= 3*N:
            dp[i+1][j+1+3*N] += dp[i][j+3*N]
            dp[i+1][j+1+3*N] %= M
        if i+2 <= 3*N:
            dp[i+2][j-1+3*N] += (i+1)*dp[i][j+3*N]
            dp[i+2][j-1+3*N] %= M
        if i+3 <= 3*N:
            dp[i+3][j+3*N] += (i+1)*(i+2)*dp[i][j+3*N]
            dp[i+3][j+3*N] %= M
r = 0
for i in range(3*N, 6*N+1):
    r += dp[3*N][i]
print(r%M)

Solution (Python):

"""
Writer: SPD_9X2
https://atcoder.jp/contests/agc043/tasks/agc043_d

結局は、N個をX個の集合に、3個を最大として分ける通り数？
それを集合内降順ソートして、全体で昇順ソートしたものは構成可能

3こ取る集合の数、2個とる集合の数を決めれば、1この集合の数もわかる
→2重ループで解けそう

1ことるやつは、 全体C取る数
2ことるやつは、 残りC取るやつ　* (Π(i=1~N) 2i C 2) / (N!)
3ことるやつは　  ( Π(i=1~N) 3i C 3 ) / (N!)

Π(i=1~N) 2i C 2 = TX[i]
Π(i=1~N) 3i C 3 = TT[i]を各Nについて前計算しておく必要性あり

グループ内の最大が先頭である必要があるけど、残りは自由
→なんか違くね…

→2の個数が1の個数を超えたらダメなのか！！


もうちょっとちゃんと考察しよう
ある順列があっったとする。これが構築可能かの判定
→大小の塊をセットにする。

塊の大きさは3以下、かつ2の塊の数が1より多かったらダメ
このような分け方のとき、かならず構築できる？

1,2の大きさを組み合わせ、先頭が昇順になるように合わせればおｋ

"""
def inverse(a,mod): 
    return pow(a,mod-2,mod)







def modfac(n, MOD):
 
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(f)
    inv = pow(f, MOD - 2, MOD)
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv
    return factorials, invs


def modnCr(n,r,mod,fac,inv): 

    return fac[n] * inv[n-r] * inv[r] % mod


N,mod = map(int,input().split())

fac,inv = modfac(10*N+30,mod)

TX = [1]

for i in range(1,2 * N+10):
    
    TX.append(TX[-1] * modnCr(2*i,2,mod,fac,inv) % mod)

TT = [1]

for i in range(1,2 * N+10):
    
    TT.append(TT[-1] * modnCr(3*i,3,mod,fac,inv) % mod)

ans = 0

for tri in range(N+1):

    for sec in range(2*N):

        one = 3*N - 3*tri - 2*sec
        if one < 0 or one < sec:
            continue

        now = modnCr(3*N,one,mod,fac,inv)
        now *= modnCr(3*N-one,2*sec,mod,fac,inv) * TX[sec] * inv[sec]
        now *= TT[tri] * inv[tri] * pow(2,tri,mod)

        now %= mod
        ans += now
        
        ans %= mod

print (ans)