Problem:
 Score : 600 points    Problem Statement Let us define the oddness of a permutation p = { p_1,\ p_2,\ ...,\ p_n } of { 1,\ 2,\ ...,\ n } as \sum_{i = 1}^n |i - p_i| .  Find the number of permutations of { 1,\ 2,\ ...,\ n } of oddness k , modulo 10^9+7 .      Constraints  All values in input are integers.  1 \leq n \leq 50  0 \leq k \leq n^2         Input Input is given from Standard Input in the following format:  n  k       Output Print the number of permutations of { 1,\ 2,\ ...,\ n } of oddness k , modulo 10^9+7 .        Sample Input 1 3 2      Sample Output 1 2  There are six permutations of { 1,\ 2,\ 3 }. Among them, two have oddness of 2 :  { 2,\ 1,\ 3 } and { 1,\ 3,\ 2 }.       Sample Input 2 39 14      Sample Output 2 74764168

Solution (Python):
import numpy as np
n, k = map(int, input().split())
m = n**2
mod = 10**9+7
dp = np.zeros((n+1, n+1, m+1), dtype=np.int64)
dp[0,0,0] = 1
for i in range(1, n+1):
  for j in range(i+1):
    dp[i,j,2*j:] = (2*j+1)*dp[i-1,j,:m+1-2*j]
    if j+1 <= i-1:
      dp[i,j,2*j:] += (j+1)*(j+1)*dp[i-1,j+1,:m+1-2*j]
    if j:
      dp[i,j,2*j:] += dp[i-1,j-1,:m+1-2*j]
    dp[i,j,2*j:] %= mod
ans = dp[n,0,k]
print(ans)

Solution (Python):
n, k = map(int, input().split())
m = n**2
mod = 10**9+7
dp = [[[0]*(m+1) for j in range(i+1)] for i in range(n+1)]
dp[0][0][0] = 1
for i in range(n):
  for j in range(i+1):
    for l in range(m+1):
      dp[i][j][l] %= mod
      if j and l+2*(j-1) <= m:
        dp[i+1][j-1][l+2*(j-1)] += dp[i][j][l]*j**2
      if l+2*j <= m:
        dp[i+1][j][l+2*j] += dp[i][j][l]*(2*j+1)
        if l+2*(j+1) <= m: 
          dp[i+1][j+1][l+2*(j+1)] += dp[i][j][l]
ans = dp[n][0][k]%mod
print(ans)

Solution (Python):
N, K = map(int,input().split())
MOD = 10**9+7




dp = []
for _ in range(55):
    dp.append([[0] * (55*55) for _ in range(55)])

dp[0][0][0] = 1

for i in range(N+1):
    for j in range(N+1):
        for k in range(N*N+1):
            if k-2*j < 0 or i-1 < 0: continue

            dp[i][j][k] += (2*j+1) * dp[i-1][j][k-2*j] % MOD
            dp[i][j][k] %= MOD

            dp[i][j][k] += (j+1)*(j+1) * dp[i-1][j+1][k-2*j] % MOD
            dp[i][j][k] %= MOD

            if j > 0:
                dp[i][j][k] += dp[i-1][j-1][k-2*j]
                dp[i][j][k] %= MOD
            
print(dp[N][0][K])

Solution (Python):

"""

https://atcoder.jp/contests/abc134/tasks/abc134_f

少ないのは何？
一応マッチングの問題といえる

dp？
数字 >= index の位置と、数字 < index の場合がある。
前の場合の数字の和を An
indexの和をAi
後者を Bn / Biとする

(An-Ai) + (Bi-Bn) = k となるようにすればいい


箱根的に左から見ていく場合を考える。
そこのindexと、そこに置いてあった数字をどうするか考える

●そこのindex
    取っておいた数字を置く→数字 < indexとなるので Biに加算 (正の寄与)
        →とっておいた数字が1以上ないとダメ
    開けておく　→ Aiに加算 (負の寄与)・後で埋めないとダメ

●そこの数字
    左の空いてるところに入れる　→　Anに加算 (正の寄与)・空きがないとダメ
    取っておく　→ Bnに加算 (負の寄与)

dp[i][score][x][y] = i番目まで見て左に空きマスがx個あり、y個数字を取ってある場合の場合の数
dp[n][score][0][0] が答え
空きマスの数ととってある数字の個数は同じなので添え字を減らせる

→ dp[i][score][x] = i番目まで見て左に空きマスがx個あり、x個数字を取ってある場合の場合の数
dp[n][k][0] が答え
scoreの取りうる範囲は -n**2 ～ n**2でいいよね

"""

n,k = map(int,input().split())

scr = 10000

dp = [[[0] * (n+1) for i in range(scr)] for j in range(n+1)]
dp[0][0][0] = 1
mod = 10**9+7

for i in range(0,n):
    for nsc in range(-scr//2,scr//2):
        for x in range(n):

            
            nownum = i+1

            
            dp[i+1][nsc][x] += dp[i][nsc][x]
            dp[i+1][nsc][x] %= mod

            
            
            
            if x > 0:
                dp[i+1][(nsc+2*nownum)%scr][x-1] += dp[i][nsc][x] * x * x
                dp[i+1][(nsc+2*nownum)%scr][x-1] %= mod

            
            
            
            if x > 0:
                dp[i+1][nsc][x] += dp[i][nsc][x] * x
                dp[i+1][nsc][x] %= mod

            
            
            
            if x > 0:
                dp[i+1][nsc][x] += dp[i][nsc][x] * x
                dp[i+1][nsc][x] %= mod

            
            
            
            if x < n:
                dp[i+1][(nsc-2*nownum)%scr][x+1] += dp[i][nsc][x]
                dp[i+1][(nsc-2*nownum)%scr][x+1] %= mod

print (dp[n][k][0] % mod)

Solution (Python):
MOD=10**9+7
n,kk=map(int,input().split())

dp=[[[0]*(n**2+1) for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0]=1

for i in range(1,n+1):
  for j in range(i+1):
    for k in range(n**2+1):
      if k-2*j<0:
        continue
      dp[i][j][k]=(2*j+1)*dp[i-1][j][k-2*j]  
      if j+1<=n:
        dp[i][j][k]+=(j+1)**2*dp[i-1][j+1][k-2*j]
      if j-1>=0:
        dp[i][j][k]+=dp[i-1][j-1][k-2*j]
      dp[i][j][k]%=MOD
      
print(dp[n][0][kk])

Solution (Python):
import sys
def input():
	return sys.stdin.readline()[:-1]
MOD = 10**9+7
n, K = map(int, input().split())
dp = [[[0 for _ in range(K+1)] for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0] = 1

for i in range(1, n+1):
	for j in range(n+1):
		for k in range(K+1):
			if k >= 2*j:
				dp[i][j][k] += (2*j+1) * dp[i-1][j][k-2*j]
			if k >= 2*j and j < n:
				dp[i][j][k] += (j+1) * (j+1) * dp[i-1][j+1][k-2*j]
			if k >= 2*j and j > 0:
				dp[i][j][k] += dp[i-1][j-1][k-2*j]
			dp[i][j][k] %= MOD

print(dp[n][0][K])

Solution (Python):
import sys
def input():
	return sys.stdin.readline()[:-1]
MOD = 10**9+7
n, K = map(int, input().split())
dp = [[[0 for _ in range(K+1)] for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0] = 1

for i in range(1, n+1):
	for j in range(n+1):
		for k in range(K+1):
			if k >= 2*j:
				dp[i][j][k] += (2*j+1) * dp[i-1][j][k-2*j]
			if k >= 2*j and j < n:
				dp[i][j][k] += (j+1) * (j+1) * dp[i-1][j+1][k-2*j]
			if k >= 2*j and j > 0:
				dp[i][j][k] += dp[i-1][j-1][k-2*j]
			dp[i][j][k] %= MOD

print(dp[n][0][K])

Solution (Python):
import sys

input = sys.stdin.readline
sys.setrecursionlimit(100000)
mod = 10 ** 9 + 7


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


def read_index():
    return map(lambda x: int(x) - 1, input().split())


def read_list():
    return list(read_values())


def read_lists(N):
    return [read_list() for n in range(N)]


class V:
    def __init__(self, f, v=None):
        self.f = f
        self.v = v

    def __str__(self):
        return str(self.v)

    def ud(self, n):
        if n is None:
            return

        if self.v is None:
            self.v = n
            return
        self.v = self.f(self.v, n)


def main():
    N, K = read_values()
    dp = [
        [[0 for i in range(N ** 2 + 1)] for ii in range(N + 1)] for iii in range(N + 1)
    ]
    dp[0][0][0] = 1

    for i in range(1, N + 1):
        for j in range(N + 1):
            for k in range(K + 1):
                if k - 2 * j < 0:
                    continue

                dp[i][j][k] = (
                    (2 * j + 1) * dp[i - 1][j][k - 2 * j]
                    + (
                        (j + 1) * (j + 1) * dp[i - 1][j + 1][k - 2 * j]
                        if j + 1 <= N
                        else 0
                    )
                    + dp[i - 1][j - 1][k - 2 * j]
                ) % mod

    print(dp[N][0][K])


if __name__ == "__main__":
    main()

Solution (Python):
import copy
n,m=map(int,input().split())
a=[[0 for _ in range(m+3)] for _ in range(n+2)]
a[0][0]=1
a[1][2]=1
for i in range(1,n):
  b=copy.deepcopy(a)
  for j in range(0,i+2):
    for k in range(0,m+1):
      if k>=2*j:
        a[j][k]=(2*j+1)*b[j][k-2*j]+(j+1)**2*b[j+1][k-2*j]
        if j>0:
          a[j][k]+=b[j-1][k-2*j]
        a[j][k]%=(10**9+7)
      else:
        a[j][k]=0
print(a[0][m])

Solution (Python):
def main():
    import sys
    input = sys.stdin.readline

    mod = 1000000007
    N, K = map(int, input().split())
    n = N//2
    dp = [[[0] * (N**2+1) for _ in range(n+1)] for _ in range(N+1)]
    dp[0][0][0] = 1

    for i in range(N):
        for j in range(n+1):
            for k in range(N**2+1):
                
                if k+2*j <= N**2:
                    dp[i+1][j][k+2*j] = (dp[i+1][j][k+2*j] + dp[i][j][k])%mod
                    if j+1 <= n:
                        dp[i+1][j+1][k+2*j] = (dp[i+1][j+1][k+2*j] + dp[i][j][k])%mod
                    dp[i + 1][j][k + 2 * j] = (dp[i + 1][j][k + 2 * j] + (dp[i][j][k]*2 * j)%mod) % mod
                    dp[i+1][j-1][k+2*j] = (dp[i+1][j-1][k+2*j] + (dp[i][j][k]*j*j)%mod)%mod
    print(dp[-1][0][K])


if __name__ == '__main__':
    main()

Solution (Python):
mod = 10 ** 9 + 7

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

dp = [[[0 for _ in range(2610)] for _ in range(51)] for _ in range(51)]
dp[0][0][0] = 1

for i in range(N):
    for j in range(N):
        for k in range(K + 1):
            dp[i + 1][j + 1][k + 2 * j + 2] += dp[i][j][k]
            if j > 0:
                dp[i + 1][j][k + 2 * j] += dp[i][j][k] * j
                dp[i + 1][j][k + 2 * j] += dp[i][j][k] * j

            dp[i + 1][j][k + 2 * j] += dp[i][j][k]

            if j > 0:
                dp[i + 1][j - 1][k + 2 * j - 2] += dp[i][j][k] * j * j


print(dp[N][0][K] % mod)

Solution (Python):
def permutation_oddness(n: int, k: int)->int:
    MOD = 10**9 + 7
    MAX_K = n*n//2

    if k > MAX_K or k % 2 == 1:
        return 0

    dp = [[[0] * (MAX_K+1) for _ in range(n+1)] for _ in range(n+1)]
    dp[0][0][0] = 1
    for i in range(n):
        for j in range(n+1):
            for odd in range(0, MAX_K+1, 2):
                if odd < 2*j:
                    continue

                dp[i+1][j][odd] = (2*j+1) * dp[i][j][odd-2*j] % MOD
                if j+1 <= n:
                    dp[i+1][j][odd] += (j+1)*(j+1)*dp[i][j+1][odd-2*j] % MOD
                if j-1 >= 0:
                    dp[i+1][j][odd] += dp[i][j-1][odd-2*j] % MOD
                dp[i+1][j][odd] %= MOD

    return dp[n][0][k]


if __name__ == "__main__":
    n, k = map(int, input().split())
    ans = permutation_oddness(n, k)
    print(ans)

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=10**9+7
input=lambda :sys.stdin.readline().rstrip()
def resolve():
    n,K=map(int,input().split())
    dp=[[[0]*2700 for _ in range(55)] for _ in range(55)]
    dp[0][0][0]=1
    for i in range(n):
        for j in range(i+1):
            for k in range(K+1):
                dp[i+1][j][k+2*j]+=dp[i][j][k]*(2*j+1)
                dp[i+1][j][k+2*j]%=MOD
                dp[i+1][j+1][k+2*(j+1)]+=dp[i][j][k]
                dp[i+1][j+1][k+2*(j+1)]%=MOD
                if(j):
                    dp[i+1][j-1][k+2*(j-1)]+=dp[i][j][k]*j*j
                    dp[i+1][j-1][k+2*(j-1)]%=MOD
    print(dp[n][0][K])
resolve()

Solution (Python):
import sys
sys.setrecursionlimit(2147483647)
INF=float("inf")
MOD=10**9+7
input=lambda :sys.stdin.readline().rstrip()
def resolve():
    n,K=map(int,input().split())
    dp=[[[0]*2700 for _ in range(55)] for _ in range(55)]
    dp[0][0][0]=1
    for i in range(n):
        for j in range(i+1):
            for k in range(K+1):
                dp[i+1][j][k+2*j]+=(2*j+1)*dp[i][j][k]
                dp[i+1][j][k+2*j]%=MOD
                dp[i+1][j+1][k+2*(j+1)]+=dp[i][j][k]
                dp[i+1][j+1][k+2*(j+1)]%=MOD
                if(j):
                    dp[i+1][j-1][k+2*(j-1)]+=j*j*dp[i][j][k]
                    dp[i+1][j-1][k+2*(j-1)]%=MOD
    print(dp[n][0][K])
resolve()

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
from bisect import bisect_left, bisect_right
import random
from itertools import permutations, accumulate, combinations
import sys
import string



INF = float('inf')
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.readline().split()
def S(): return sys.stdin.readline().strip()
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 10 ** 9 + 7


n, k = LI()
dp = [[[0]*(k+2*(n+1)) for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0] = 1
for i in range(n):
    for j in range(i+1):
        for l in range(k+1):
            dp[i+1][j+1][l+2*(j+1)] += dp[i][j][l]
            dp[i+1][j+1][l+2*(j+1)]%=mod
            dp[i+1][j][l+2*j] += dp[i][j][l]*(j*2+1)
            dp[i+1][j][l+2 * j] %= mod
            if j!= 0:
                dp[i+1][j-1][l+(j-1)*2] += dp[i][j][l]*j*j
                dp[i+1][j-1][l+(j-1)*2]%=mod



print(dp[n][0][k]%mod)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import math
from bisect import bisect_left, bisect_right
import random
from itertools import permutations, accumulate, combinations
import sys
import string



INF = float('inf')
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS(): return sys.stdin.readline().split()
def S(): return sys.stdin.readline().strip()
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 10 ** 9 + 7


n, k = LI()
dp = [[[0]*n**2*4 for _ in range(n+1)] for _ in range(n+1)]
dp[0][0][0] = 1
for i in range(n):
    for j in range(i+1):
        for l in range(k+1):
            dp[i+1][j+1][l+2*(j+1)] += dp[i][j][l]
            dp[i+1][j+1][l+2*(j+1)]%=mod
            dp[i+1][j][l+2*j] += dp[i][j][l]*(j*2+1)
            dp[i+1][j][l+2 * j] %= mod
            if j!= 0:
                dp[i+1][j-1][l+(j-1)*2] += dp[i][j][l]*j*j
                dp[i+1][j-1][l+(j-1)*2]%=mod



print(dp[n][0][k]%mod)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

MOD = 10**9 + 7

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


L = N*(N-1)//2+1
dp = np.zeros((N,L),dtype=np.int64)
dp[0,0] = 1

for _ in range(1,N+1):
    prev = dp
    dp = np.zeros_like(prev)
    for n in range(N):
        x = prev[n]
        
        dp[n,n:] += x[:L-n] * (n+n+1)
        
        if n > 0: dp[n-1,n-1:] += x[:L-n+1] * n
        
        if n+1 < N: dp[n+1,n+1:] += x[:L-n-1] * (n+1)
    dp %= MOD

x = dp[0]
if K&1:
    answer = 0
elif K//2 < L:
    answer = x[K//2]
else:
    answer = 0
print(answer)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

import numpy as np

MOD = 10**9 + 7

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


L = N*(N-1)//2+1
dp = np.zeros((N,L),dtype=np.int64)
dp[0,0] = 1

for _ in range(1,N+1):
    prev = dp
    dp = np.zeros_like(prev)
    for n in range(N):
        
        dp[n,n:] += prev[n,:L-n]
        
        if n > 0: dp[n-1,n-1:] += prev[n,:L-n+1] * n
        
        dp[n,n:] += prev[n,:L-n] * n
        
        dp[n,n:] += prev[n,:L-n] * n
        
        if n+1 < N: dp[n+1,n+1:] += prev[n,:L-n-1] * (n+1)
    dp %= MOD

x = dp[0]
if K&1:
    answer = 0
elif K//2 < L:
    answer = x[K//2]
else:
    answer = 0
print(answer)

Solution (Python):
n, K = map(int, input().split())
dp = [[0]*(n**2+1) for _ in range(n+2)]
dp[0][0] = 1
for i in range(n):
    ndp = [[0]*(n**2+1) for _ in range(n+2)]
    for j in range(i+2):
        for k in range(n**2+1):
            if k-2*j<0:
                continue
            ndp[j][k] = ((2*j+1)*dp[j][k-2*j]+(j+1)*(j+1)*dp[j+1][k-2*j]+dp[j-1][k-2*j])%(10**9+7)
    dp = ndp
print(dp[0][K])

Solution (Python):
import os
import sys

import numpy as np

if os.getenv("LOCAL"):
    sys.stdin = open("_in.txt", "r")

sys.setrecursionlimit(2147483647)
INF = float("inf")
IINF = 10 ** 18
MOD = 10 ** 9 + 7

N, K = list(map(int, sys.stdin.readline().split()))






















































dp = np.zeros((N + 1, K + 1, N + 1), dtype=int)
dp[0, 0, 0] = 1
for i in range(N):
    for a in range(N + 1):
        k = np.arange(K + 1)
        
        if a == 0:
            dp[i + 1, :, a] += dp[i, :, a]
        else:
            dp[i + 1, a * 2:, a] += dp[i, :-(a * 2), a]

        
        
        if a * 2 - 2 == 0:
            dp[i + 1, :, a - 1] += dp[i, :, a] * (a ** 2)
        elif a * 2 - 2 > 0:
            dp[i + 1, a * 2 - 2:, a - 1] += dp[i, :-(a * 2 - 2), a] * (a ** 2)

        
        if a > 0:
            dp[i + 1, a * 2:, a] += dp[i, :-(a * 2), a] * a * 2

        
        if a + 1 < N:
            dp[i + 1, a * 2 + 2:, a + 1] += dp[i, :-(a * 2 + 2), a]
    dp[i + 1] %= MOD
print(dp[-1][K][0] % MOD)

Solution (Python):
import numpy as np


def solve(n, k):
    if k % 2 == 1:
        return 0
    k //= 2

    MOD = 10 ** 9 + 7

    dp = np.zeros((1, k + 1), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(1, n + 1):
        max_d = min(i + 1, n - i + 1, k + 1)
        ndp = np.zeros((max_d, k + 1), dtype=np.int64)
        for d, ks in enumerate(dp):
            base = ks[:k - d + 1]
            if d > 0:
                ndp[d - 1, d:] += base * d ** 2
            if max_d > d:
                ndp[d, d:] += base * (2 * d + 1)
            if max_d > d + 1:
                ndp[d + 1, d:] += base
        dp = ndp % MOD

    return dp[0, k]


n, k = map(int, input().split())
print(solve(n, k))