Problem:
 Score : 100 points    Problem Statement There are N men and N women, both numbered 1, 2, \ldots, N .  For each i, j ( 1 \leq i, j \leq N ), the compatibility of Man i and Woman j is given as an integer a_{i, j} .
If a_{i, j} = 1 , Man i and Woman j are compatible; if a_{i, j} = 0 , they are not.  Taro is trying to make N pairs, each consisting of a man and a woman who are compatible.
Here, each man and each woman must belong to exactly one pair.  Find the number of ways in which Taro can make N pairs, modulo 10^9 + 7 .      Constraints  All values in input are integers.  1 \leq N \leq 21  a_{i, j} is 0 or 1 .         Input Input is given from Standard Input in the following format:  N  a_{1, 1}  \ldots  a_{1, N}  :  a_{N, 1}  \ldots  a_{N, N}       Output Print the number of ways in which Taro can make N pairs, modulo 10^9 + 7 .        Sample Input 1 3
0 1 1
1 0 1
1 1 1      Sample Output 1 3  There are three ways to make pairs, as follows ( (i, j) denotes a pair of Man i and Woman j ):   (1, 2), (2, 1), (3, 3)  (1, 2), (2, 3), (3, 1)  (1, 3), (2, 1), (3, 2)        Sample Input 2 4
0 1 0 0
0 0 0 1
1 0 0 0
0 0 1 0      Sample Output 2 1  There is one way to make pairs, as follows:   (1, 2), (2, 4), (3, 1), (4, 3)        Sample Input 3 1
0      Sample Output 3 0       Sample Input 4 21
0 0 0 0 0 0 0 1 1 0 1 1 1 1 0 0 0 1 0 0 1
1 1 1 0 0 1 0 0 0 1 0 0 0 0 1 1 1 0 1 1 0
0 0 1 1 1 1 0 1 1 0 0 1 0 0 1 1 0 0 0 1 1
0 1 1 0 1 1 0 1 0 1 0 0 1 0 0 0 0 0 1 1 0
1 1 0 0 1 0 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0
0 1 1 0 1 1 1 0 1 1 1 0 0 0 1 1 1 1 0 0 1
0 1 0 0 0 1 0 1 0 0 0 1 1 1 0 0 1 1 0 1 0
0 0 0 0 1 1 0 0 1 1 0 0 0 0 0 1 1 1 1 1 1
0 0 1 0 0 1 0 0 1 0 1 1 0 0 1 0 1 0 1 1 1
0 0 0 0 1 1 0 0 1 1 1 0 0 0 0 1 1 0 0 0 1
0 1 1 0 1 1 0 0 1 1 0 0 0 1 1 1 1 0 1 1 0
0 0 1 0 0 1 1 1 1 0 1 1 0 1 1 1 0 0 0 0 1
0 1 1 0 0 1 1 1 1 0 0 0 1 0 1 1 0 1 0 1 1
1 1 1 1 1 0 0 0 0 1 0 0 1 1 0 1 1 1 0 0 1
0 0 0 1 1 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 1
1 0 1 1 0 1 0 1 0 0 1 0 0 1 1 0 1 0 1 1 0
0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 1 0 0 1
0 0 0 1 0 0 1 1 0 1 0 1 0 1 1 0 0 1 1 0 1
0 0 0 0 1 1 1 0 1 0 1 1 1 0 1 1 0 0 1 1 0
1 1 0 1 1 0 0 1 1 0 1 1 0 1 1 1 1 1 0 1 0
1 0 0 1 1 0 1 1 1 1 1 0 1 0 1 1 0 0 0 0 0      Sample Output 4 102515160  Be sure to print the number modulo 10^9 + 7 .

Solution (Python):
n = int(input())
g=[]
for i in range(n):
    
    g.append(list(map(int,input().split())))
mod=10**9+7
dp=[0]*(2**n)
dp[0]=1
for i in range(1,2**n):
    e=bin(i)[2:]
    y=e.count('1')
    for j in range(len(e)):
        if e[-j-1]=='1':
            if g[y-1][j]:
                    dp[i]+=dp[i-(2**j)]
                    dp[i]%=mod


print(dp[-1]%mod)

Solution (Python):
N = int(input())
Match = [0] * N
mod = 10**9 + 7
for i in range(N):
    A = list(map(int, input().split()))
    m = 0
    for j, a in enumerate(A):
        if a:
            m |= (1 << j)
    Match[i] = m

DP = [0] * (1 << N)
DP[0] = 1

bit_count = [0]
for _ in range(N):
    bit_count += [x+1 for x in bit_count]

for i in range((1 << N) - 1):
    c = bit_count[i]
    if DP[i] == 0:
        continue
    for j in range(N):
        if (Match[c] >> j) & 1 == 1:
            if (i >> j) & 1 == 0:
                DP[i| (1 << j)] += DP[i]
                DP[i| (1 << j)] %= mod
print(DP[-1] % mod)

Solution (Python):
def edu_dp_d_matching():
    n = int(input())
    a = []
    for _ in range(n):
        a.append(list(map(int, input().split())))
    dp = [[0] * (1 << n) for _ in range(n)]
    for i in range(n):
        if a[0][i] == 1: dp[0][(1 << i)] = 1
    for i in range(1, n):
        for k in range(1 << n):
            if dp[i - 1][k] == 0: continue
            for j in range(n):
                if a[i][j] == 0 or k & (1 << j) != 0: continue
                dp[i][k | (1 << j)] += dp[i - 1][k]

    print(dp[-1][-1] % (pow(10, 9) + 7))


edu_dp_d_matching()

Solution (Python):
n=int(input())
a=[list(map(int,input().split())) for i in range(n)]
mod=10**9+7

dp=[0]*(2**n)
dp[0]=1

for i in range(1,2**n):
    x=bin(i)[2:]
    y=x.count('1')
    for j in range(len(x)):
      if x[-j-1]=='1':
        if a[y-1][j]==1:
            dp[i]+=dp[i-(2**j)]
            dp[i]%=mod


print(dp[-1])

Solution (Python):
N = int(input())
A = [list(map(int, input().split())) for i in range(N)]
MOD = 10**9 + 7

memo = [-1]*(1 << N) 
memo[0] = 1
def dfs(state, c):
    if memo[state] != -1:
        return memo[state]
    r = 0
    for i in range(N):
        if state & (1 << i) and A[c][i]:
            r += dfs(state ^ (1 << i), c+1) 
    memo[state] = r = r % MOD
    return r
print(dfs((1 << N)-1, 0))

Solution (Python):
MOD = 1000000007

N = int(input())
A = [tuple(map(int, input().split())) for _ in range(N)]

dp = [0 for _ in range(1 << N)]
dp[0] = 1

for i in range(1, 1 << N):
    k = 0
    for j in range(N):
        if (i >> j) & 1:
            k += 1
    for j in range(N):
        if (i >> j) & 1 and A[k - 1][j]:
            dp[i] += dp[i - (1 << j)]
            dp[i] %= MOD

print(dp[-1])

Solution (Python):
MOD = 1000000007

N = int(input())
A = [tuple(map(int, input().split())) for _ in range(N)]

dp = [0 for _ in range(1 << N)]
dp[0] = 1

for i in range(1, 1 << N):
    k = 0
    for j in range(N):
        if (i >> j) & 1:
            k += 1
    for j in range(N):
        if (i >> j) & 1 and A[k - 1][j]:
            dp[i] += dp[i ^ (1 << j)]
            dp[i] %= MOD

print(dp[-1])

Solution (Python):
n=int(input())
a=[]
for i in range(n):
    a.append(list(map(int,input().split())))

def b(n,k):
    return (n>>k)&1
MOD=10**9+7

dp=[[0 for j in range(1<<n)] for i in range(n+1)]
dp[0][0]=1
for i in range(n):
    for S in range(1<<n):
        if (dp[i][S]):
            for j in range(n):
                if (b(S,j)==0 and a[i][j]==1):
                    dp[i+1][S^(1<<j)]+=(dp[i][S^(1<<j)]+dp[i][S])%MOD
                    dp[i+1][S^(1<<j)]%=MOD
print(dp[n][-1])

Solution (Python):
N=int(input())
A=[list(map(int,input().split())) for i in range(N)]
DP=[-1]*(1<<N)
def f(x,k):
  if x==0:
    DP[x]=1
  if DP[x]!=-1:
    return DP[x]
  DP[x]=0
  for i in range(N):
    if x&(1<<i):
      if A[k][i]:
        DP[x]+=f(x^(1<<i),k+1)
  return DP[x]

print(f((1<<N)-1,0)%(10**9+7))

Solution (Python):
import sys
input = sys.stdin.readline
sys.setrecursionlimit(10**7)

mod = 10**9+7
n = int(input())
A = tuple(tuple(map(int, input().split())) for _ in range(n))
dp = [[-1]*(1<<n) for _ in range(n+1)]
def f(i, k):
  if i == n:
    return 1
  if dp[i][k] != -1:
    return dp[i][k]
  res = 0
  for j in range(n):
    if A[i][j] and ((k>>j)&1):
      res += f(i+1, k^(1<<j))
      res %= mod
  dp[i][k] = res
  return res
ans = f(0, (1<<n)-1)
print(ans)

Solution (Python):
n=int(input())
a=[list(map(int,input().split())) for _ in range(n)]
dp=[0]*2**n
mod=10**9+7
dp[0]=1
for i in range(2**n-1):
    cnt=-1
    for j in range(n):
        if (i>>j)&1:
            cnt+=1
    for j in range(n):
        if (i>>j)&1==0:
            if a[cnt+1][j]:
                dp[i|pow(2,j)]+=dp[i]
                dp[i|pow(2,j)]%=mod
print(dp[2**n-1])

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)

ra = range
enu = enumerate


def exit(*argv, **kwarg):
    print(*argv, **kwarg)
    sys.exit()


def mina(*argv, sub=1): return list(map(lambda x: x - sub, argv))



def a_int(): return int(readline())


def read_matrix(H, W):
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


from collections import defaultdict, Counter, deque
from operator import itemgetter, xor, add
from itertools import product, permutations, combinations
from bisect import bisect_left, bisect_right  
from functools import reduce

N = a_int()
A = read_matrix(N, N)

'''
一方の集合を管理してDPというアイデア

dp[i][s] ... 男性を:i人目まで見たときに、マッチングした女性の集合がsであるような通りの総数

更新式
dp[i+1][s+1<<j] += dp[i][s] for s のjbit目は0 かつ A[i,j]==1

これでは計算量がO(2^N 2N)なのでTLE
iと|s|は同じなのだから状態をまとめることができる
(二次元dpでいうとi=|s|のセルでしか状態遷移が発生しない)

よって単純sだけでループしても良さそうだし、()内のアイデアで枝刈りしても良さそう
'''

from numba import njit


@njit('(i8,i8[:,:])', cache=True)
def solve_TLE(N, A):
    dp = np.zeros((N + 1, 1 << N), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(N):
        for s in range(1 << N):
            for j in range(N):
                if (s >> j) & 1 == 0 and A[i, j] == 1:
                    dp[i + 1, s + (1 << j)] += dp[i, s]
                    dp[i + 1, s + (1 << j)] %= MOD
    print(dp[N, (1 << N) - 1] % MOD)


@njit('i8(i8)', cache=True)
def popcount(x):
    '''xの立っているビット数をカウントする関数
    (xは64bit整数)'''

    
    x = x - ((x >> 1) & 0x5555555555555555)

    
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)

    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f  
    x = x + (x >> 8)  
    x = x + (x >> 16)  
    x = x + (x >> 32)  
    return x & 0x0000007f


@njit('(i8,i8[:,:])', cache=True)
def solve(N, A):  
    dp = np.zeros((N + 1, 1 << N), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(N):
        for s in range(1 << N):
            if popcount(s) != i:  
                continue
            
            
            for j in range(N):
                if (s >> j) & 1 == 0 and A[i, j] == 1:
                    dp[i + 1, s + (1 << j)] += dp[i, s]
                    dp[i + 1, s + (1 << j)] %= MOD
    print(dp[N, (1 << N) - 1] % MOD)


@njit('(i8,i8[:,:])', cache=True)
def solve2(N, A):  
    dp = np.zeros((1 << N), dtype=np.int64)
    dp[0] = 1
    for s in range(1 << N):
        for j in range(N):
            if (s >> j) & 1 == 0 and A[popcount(s), j] == 1:
                dp[s + (1 << j)] += dp[s]
                dp[s + (1 << j)] %= MOD
    print(dp[(1 << N) - 1] % MOD)



solve(N, A)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)

ra = range
enu = enumerate


def exit(*argv, **kwarg):
    print(*argv, **kwarg)
    sys.exit()


def mina(*argv, sub=1): return list(map(lambda x: x - sub, argv))



def a_int(): return int(readline())


def read_matrix(H, W):
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


from collections import defaultdict, Counter, deque
from operator import itemgetter, xor, add
from itertools import product, permutations, combinations
from bisect import bisect_left, bisect_right  
from functools import reduce

N = a_int()
A = read_matrix(N, N)

'''
一方の集合を管理してDPというアイデア

dp[i][s] ... 男性を:i人目まで見たときに、マッチングした女性の集合がsであるような通りの総数

更新式
dp[i+1][s+1<<j] += dp[i][s] for s のjbit目は0 かつ A[i,j]==1

これでは計算量がO(2^N 2N)なのでTLE
iと|s|は同じなのだから状態をまとめることができる
(二次元dpでいうとi=|s|のセルでしか状態遷移が発生しない)

よって単純sだけでループしても良さそうだし、()内のアイデアで枝刈りしても良さそう
'''

from numba import njit


@njit('(i8,i8[:,:])', cache=True)
def solve_TLE(N, A):
    dp = np.zeros((N + 1, 1 << N), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(N):
        for s in range(1 << N):
            for j in range(N):
                if (s >> j) & 1 == 0 and A[i, j] == 1:
                    dp[i + 1, s + (1 << j)] += dp[i, s]
                    dp[i + 1, s + (1 << j)] %= MOD
    print(dp[N, (1 << N) - 1] % MOD)


@njit('(i8,i8[:,:])', cache=True)
def solve(N, A):  
    dp = np.zeros((N + 1, 1 << N), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(N):
        for s in range(1 << N):
            if dp[i, s] == 0:  
                continue
            for j in range(N):
                if (s >> j) & 1 == 0 and A[i, j] == 1:
                    dp[i + 1, s + (1 << j)] += dp[i, s]
                    dp[i + 1, s + (1 << j)] %= MOD
    print(dp[N, (1 << N) - 1] % MOD)


@njit('i8(i8)', cache=True)
def popcount(x):
    '''xの立っているビット数をカウントする関数
    (xは64bit整数)'''

    
    x = x - ((x >> 1) & 0x5555555555555555)

    
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)

    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f  
    x = x + (x >> 8)  
    x = x + (x >> 16)  
    x = x + (x >> 32)  
    return x & 0x0000007f


@njit('(i8,i8[:,:])', cache=True)
def solve2(N, A):  
    dp = np.zeros((1 << N), dtype=np.int64)
    dp[0] = 1
    for s in range(1 << N):
        for j in range(N):
            if (s >> j) & 1 == 0 and A[popcount(s), j] == 1:
                dp[s + (1 << j)] += dp[s]
                dp[s + (1 << j)] %= MOD
    print(dp[(1 << N) - 1] % MOD)


solve2(N, A)

Solution (Python):
import sys
sys.setrecursionlimit(1 << 25)
readline = sys.stdin.buffer.readline
read = sys.stdin.readline  

import numpy as np
from functools import partial
array = partial(np.array, dtype=np.int64)
zeros = partial(np.zeros, dtype=np.int64)
full = partial(np.full, dtype=np.int64)

ra = range
enu = enumerate


def exit(*argv, **kwarg):
    print(*argv, **kwarg)
    sys.exit()


def mina(*argv, sub=1): return list(map(lambda x: x - sub, argv))



def a_int(): return int(readline())


def read_matrix(H, W):
    '''return np.ndarray shape=(H,W) matrix'''
    lines = []
    for _ in range(H):
        lines.append(read())
    lines = ' '.join(lines)  
    return np.fromstring(lines, sep=' ', dtype=np.int64).reshape(H, W)


MOD = 10**9 + 7
INF = 2**31  


from collections import defaultdict, Counter, deque
from operator import itemgetter, xor, add
from itertools import product, permutations, combinations
from bisect import bisect_left, bisect_right  
from functools import reduce

N = a_int()
A = read_matrix(N, N)

'''
一方の集合を管理してDPというアイデア

dp[i][s] ... 男性を:i人目まで見たときに、マッチングした女性の集合がsであるような通りの総数

更新式
dp[i+1][s+1<<j] += dp[i][s] for s のjbit目は0 かつ A[i,j]==1

これでは計算量がO(2^N 2N)なのでTLE
iと|s|は同じなのだから状態をまとめることができる
(二次元dpでいうとi=|s|のセルでしか状態遷移が発生しない)

よって単純sだけでループしても良さそうだし、()内のアイデアで枝刈りしても良さそう
'''

from numba import njit


@njit('(i8,i8[:,:])', cache=True)
def solve_TLE(N, A):
    dp = np.zeros((N + 1, 1 << N), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(N):
        for s in range(1 << N):
            for j in range(N):
                if (s >> j) & 1 == 0 and A[i, j] == 1:
                    dp[i + 1, s + (1 << j)] += dp[i, s]
                    dp[i + 1, s + (1 << j)] %= MOD
    print(dp[N, (1 << N) - 1] % MOD)


@njit('(i8,i8[:,:])', cache=True)
def solve(N, A):  
    dp = np.zeros((N + 1, 1 << N), dtype=np.int64)
    dp[0, 0] = 1
    for i in range(N):
        for s in range(1 << N):
            if dp[i, s] == 0:  
                continue
            for j in range(N):
                if (s >> j) & 1 == 0 and A[i, j] == 1:
                    dp[i + 1, s + (1 << j)] += dp[i, s]
                    dp[i + 1, s + (1 << j)] %= MOD
    print(dp[N, (1 << N) - 1] % MOD)


solve(N, A)

Solution (Python):

def countSetBits(n):
    count = 0
    while (n):
        n &= (n - 1)
        count += 1

    return count



N=10**9 + 7
n=int(input())
ar = [list(map(int,input().split())) for _ in range(n)]   
total=1<<n
bits=[[]for _ in range(n+1)]
dp=[[0]*(total+5) for j in range(n+5)]

for i in range(n):
    if ar[0][i]==1:
        dp[0][1<<(n-1-i)]=1

for j in range(total):
    i=countSetBits(j)-1
    for k in range(0,n):
        if ar[i][n-k-1]==1 and j&(1<<k):
            dp[i][j]=(dp[i][j]+dp[i-1][j&(~(1<<k))])%N

print(dp[n-1][total-1])

Solution (Python):
import sys
import numpy as np

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

MOD = 1_000_000_007

def permanent(A):
    
    N = A.shape[0]

    
    products = np.ones(1 << N, np.int64)
    for r in range(N):
        
        row_sum = np.zeros(1 << N, np.int64)
        for n in range(N):
            row_sum[1 << n:1 << n + 1] = row_sum[:1 << n] + A[r, n]
        products *= row_sum
        products %= MOD

    
    sgn = np.ones(1 << N, np.int64)
    for n in range(N):
        sgn[1 << n:1 << n + 1] = -sgn[:1 << n]

    products = products * sgn[::-1]

    perm = products.sum() % MOD
    return perm

N = int(readline())
A = np.array(read().split(), np.int64).reshape(N, N)

print(permanent(A))

Solution (Python):
import numpy as np

def main():
    MOD = 10**9 + 7
    
    N = int(input())
    aff = np.array([map(int, input().split()) for i in range(N)])
    dpt = np.zeros((N+1, 1<<N), dtype='int64')
    dpt[0][0] = 1

    for i, row in enumerate(aff):
        for idx, v in enumerate(row):
            if v == 1:
                dpt[i+1][(1<<idx):] += dpt[i][:-(1<<idx)]
        dpt[i+1] %= MOD

    print(dpt[N][(1<<N)-1])

if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np
from numba import njit
MOD = 10**9 + 7

def main():
 
    N = int(input())
    aff = np.array([list(map(int, input().split())) for i in range(N)], dtype='int32')
    dpt = np.zeros((N+1, 1<<N), dtype='int64')
    dpt[0][0] = 1

    print(solve(N, aff, dpt))


def solve(N, aff, dpt):
    for i, row in enumerate(aff):
        for idx, v in enumerate(row):
            if v == 1:
                dpt[i+1][(1<<idx):] += dpt[i][:-(1<<idx)]
        dpt[i+1] %= MOD
    return dpt[N][(1<<N)-1]

if __name__ == '__main__':
    main()

Solution (Python):
import numpy as np

MOD = 10**9 + 7
 
N = int(input())
aff = np.array([map(int, input().split()) for i in range(N)])
dpt = np.zeros((N+1, 1<<N), dtype='int64')
dpt[0][0] = 1

for i, row in enumerate(aff):
    for idx, v in enumerate(row):
        if v == 1:
            dpt[i+1][(1<<idx):] += dpt[i][:-(1<<idx)]
    dpt[i+1] %= MOD

print(dpt[N][(1<<N)-1])

Solution (Python):
N=int(input())
a=[[int(i) for i in input().split()] for _ in range(N)]
MOD=10**9+7

dp=[0 for _ in range(1<<N)]
dp[0]=1

for S in range(1,1<<N):
    i=bin(S).count('1')
    for x in range(N):
        if (S>>x & 1) and a[i-1][x]==1:
            dp[S]=(dp[S]+dp[S^(1<<x)])%MOD

print(dp[-1])

Solution (Python):
n = int(input())
A = [list(map(int,input().split())) for _ in range(n)]
dp = [0]*(1<<n)
dp[0] = 1
mod = 10**9+7
for i in range(1,1<<n):
    cnt = 0
    for j in range(n+1): cnt += (i>>j)&1
    for j in range(n):
        if (i>>j)&1==0: continue
        dp[i] +=dp[i-(1<<j)]*A[cnt-1][j]
        dp[i] %=mod
print(dp[-1]%mod)

Solution (Python):






import sys
input = sys.stdin.readline
sys.setrecursionlimit(500000)

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]


def P(arg): print(arg); return
def Y(): print("Yes"); return
def N(): print("No"); return
def E(): exit()
def PE(arg): print(arg); exit()
def YE(): print("Yes"); exit()
def NE(): print("No"); exit()


def DD(arg): return defaultdict(arg)


def inv(n): return pow(n, MOD-2, MOD)


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count




N = I()
data = ILL(N)


dp = [0 for b in range(1<<N)]

dp[0] = 1

for b in range(1,1<<N):
  i = count_bit(b)
  for j in range(N):
    if (b>>j)&1 and data[i-1][j]:
      dp[b] += dp[b-(1<<j)]
      dp[b] %= MOD
print(dp[-1])

Solution (Python):






import sys
input = sys.stdin.readline
sys.setrecursionlimit(500000)

import math

import bisect
import heapq
from collections import defaultdict
from collections import deque
from collections import Counter
from functools import lru_cache





MOD = 10**9+7
INF = float('inf')
AZ = "abcdefghijklmnopqrstuvwxyz"






def I(): return int(input().strip())
def S(): return input().strip()
def IL(): return list(map(int,input().split()))
def SL(): return list(map(str,input().split()))
def ILs(n): return list(int(input()) for _ in range(n))
def SLs(n): return list(input().strip() for _ in range(n))
def ILL(n): return [list(map(int, input().split())) for _ in range(n)]
def SLL(n): return [list(map(str, input().split())) for _ in range(n)]


def P(arg): print(arg); return
def Y(): print("Yes"); return
def N(): print("No"); return
def E(): exit()
def PE(arg): print(arg); exit()
def YE(): print("Yes"); exit()
def NE(): print("No"); exit()


def DD(arg): return defaultdict(arg)


def inv(n): return pow(n, MOD-2, MOD)


kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n): return kaijo_memo[n]
  if(len(kaijo_memo) == 0): kaijo_memo.append(1)
  while(len(kaijo_memo) <= n): kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]

gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n): return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0): gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n): gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]

def nCr(n,r):
  if n == r: return 1
  if n < r or r < 0: return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret


def factorization(n):
  arr = []
  temp = n
  for i in range(2, int(-(-n**0.5//1))+1):
    if temp%i==0:
      cnt=0
      while temp%i==0: 
        cnt+=1 
        temp //= i
      arr.append([i, cnt])
  if temp!=1: arr.append([temp, 1])
  if arr==[]: arr.append([n, 1])
  return arr


def make_divisors(n):
  divisors = []
  for i in range(1, int(n**0.5)+1):
    if n % i == 0:
      divisors.append(i)
      if i != n // i: 
        divisors.append(n//i)
  return divisors


def make_primes(N):
  max = int(math.sqrt(N))
  seachList = [i for i in range(2,N+1)]
  primeNum = []
  while seachList[0] <= max:
    primeNum.append(seachList[0])
    tmp = seachList[0]
    seachList = [i for i in seachList if i % tmp != 0]
  primeNum.extend(seachList)
  return primeNum


def gcd(a, b):
    while b: a, b = b, a % b
    return a


def lcm(a, b):
    return a * b // gcd (a, b)


def count_bit(n):
  count = 0
  while n:
    n &= n-1
    count += 1
  return count


def base_10_to_n(X, n):
  if X//n: return base_10_to_n(X//n, n)+[X%n]
  return [X%n]

def base_n_to_10(X, n):
  return sum(int(str(X)[-i-1])*n**i for i in range(len(str(X))))

def base_10_to_n_without_0(X, n):
  X -= 1
  if X//n: return base_10_to_n_without_0(X//n, n)+[X%n]
  return [X%n]


def int_log(n, a):
  count = 0
  while n>=a:
    n //= a
    count += 1
  return count




N = I()
data = ILL(N)


dp = [[0 for b in range(1<<N)] for i in range(N+1)]

dp[0][0] = 1

for i in range(N):
  for b in range(1<<N):
    if dp[i][b]:
      for j in range(N):
        if data[i][j] and (b>>j)&1 == 0:
          dp[i+1][b+(1<<j)] += dp[i][b]
          dp[i+1][b+(1<<j)] %= MOD
print(dp[-1][-1])

Solution (Python):
from sys import stdin
import sys
import math
from functools import reduce
import functools
import itertools
from collections import deque,Counter,defaultdict
from operator import mul
import copy


import heapq
sys.setrecursionlimit(10**6)
INF =  float("inf")
import bisect
import statistics
mod = 10**9+7

def popcount(x):
    '''xの立っているビット数をカウントする関数
    (xは64bit整数)'''

    
    x = x - ((x >> 1) & 0x5555555555555555)

    
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)

    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f 
    x = x + (x >> 8) 
    x = x + (x >> 16) 
    x = x + (x >> 32) 
    return x & 0x0000007f

def bit(S, j):
    
    return (S>>j)&1

N = int(input())
a = [[] for i in range(N)]
for i in range(N):
    a[i] = list(map(int, input().split()))

dp = [0 for i in range(1<<N)]
dp[0] = 1

for S in range(1,1<<N):
    i = popcount(S)
    for j in range(N):
        if bit(S,j) == 1 and a[i-1][j] == 1:
            dp[S] = (dp[S] + dp[S^(1<<j)]) % mod

print(dp[(1<<N)-1] % mod)

Solution (Python):
from itertools import combinations
def solve():
  N = int(input())
  A = [list(map(int, input().split())) for _ in range(N)]
  mod = 10**9+7
  dp = [0]*(1<<N)
  dp[0] = 1
  for i in range(N):
    for com in combinations(range(N),i+1):
      S = 0
      for k in com:
        S += 1<<k
      for j in com:
        dp[S] += dp[S-(1<<j)]*A[i][j]
        dp[S] %= mod
  return dp[-1]
print(solve())

Solution (Python):
from itertools import groupby, accumulate, product, permutations, combinations
def solve():
  N = int(input())
  A = [list(map(int, input().split())) for _ in range(N)]
  mod = 10**9+7
  dp = [0]*(1<<N)
  dp[0] = 1
  for i in range(1,N+1):
    for com in combinations(range(N),i):
      S = 0
      for k in com:
        S += 1<<k
      for j in com:
        if A[i-1][j]==1:
          dp[S] += dp[S-(1<<j)]
          dp[S] %= mod
  return dp[-1]
print(solve())

Solution (Python):
from itertools import groupby, accumulate, product, permutations, combinations
def solve():
  N = int(input())
  A = [list(map(int, input().split())) for _ in range(N)]
  mod = 10**9+7
  dp = [[0]*(1<<N) for _ in range(N+1)]
  dp[N][0] = 1
  for i in range(N-1,-1,-1):
    for com in combinations(range(N),N-i):
      S = 0
      for k in com:
        S += 1<<k
      for j in com:
        if A[i][j]==1:
          dp[i][S] += dp[i+1][S-(1<<j)]
          dp[i][S] %= mod
  return dp[0][-1]
print(solve())

Solution (Python):
from itertools import groupby, accumulate, product, permutations, combinations
def solve():
  N = int(input())
  A = [list(map(int, input().split())) for _ in range(N)]
  mod = 10**9+7
  dp = [[0]*(1<<N) for _ in range(N+1)]
  dp[N][-1] = 1
  for i in range(N-1,-1,-1):
    for com in combinations(range(N),i):
      S = 0
      for k in com:
        S += 1<<k
      for j in range(N):
        if not S>>j&1: 
          if A[i][j]==1:
            dp[i][S] += dp[i+1][S|(1<<j)]
            dp[i][S] %= mod
  return dp[0][0]
print(solve())

Solution (Python):
from itertools import groupby, accumulate, product, permutations, combinations
def solve():
  N = int(input())
  A = [list(map(int, input().split())) for _ in range(N)]
  mod = 10**9+7
  dp = [[0]*(1<<N) for _ in range(N+1)]
  dp[N][-1] = 1
  for i in range(N-1,-1,-1):
    for com in combinations(range(N),i):
      S = 0
      for k in com:
        S += 1<<k
      for j in range(N):
        if not S>>j&1: 
          if A[i][j]==1:
            dp[i][S] += dp[i+1][S+(1<<j)]
            dp[i][S] %= mod
  return dp[0][0]
print(solve())

Solution (Python):
def solve():
  N = int(input())
  A = [list(map(int, input().split())) for _ in range(N)]
  mod = 10**9+7
  dp = [[0]*(1<<N) for _ in range(N+1)]
  dp[N][-1] = 1
  for S in range((1<<N)-1,-1,-1):
    i = sum([1&(S>>k) for k in range(N)])
    for j in range(N):
      if not (S>>j)&1: 
        if A[i][j]==1:
          dp[i][S] += dp[i+1][S+(1<<j)]
          dp[i][S] %= mod
  return dp[0][0]
print(solve())

Solution (Python):
from sys import stdin, stdout, setrecursionlimit
from collections import deque, defaultdict, Counter
from heapq import heappush, heappop
from functools import lru_cache
import math

setrecursionlimit(20000)
rl = lambda: stdin.readline()
rll = lambda: stdin.readline().split()
rli = lambda: map(int, stdin.readline().split())
rlf = lambda: map(float, stdin.readline().split())

INF, NINF = float('inf'), float('-inf')

def main():
	MOD = 10**9 + 7
	n = int(rl())
	A = []
	for _ in range(n):
		row = []
		for bit in rli():
			row.append(bit)
		A.append(row)

	popcount = lambda x: bin(x).count("1")
	dp = [0 for _ in range(1 << n)]
	dp[0] = 1
	for mask in range((1 << n)-1):
		new = popcount(mask)
		for cand in range(n):
			if A[new][cand] and not (mask & (1 << cand)):
				nmask = mask ^ (1 << cand)
				dp[nmask] += (dp[mask] % MOD)
				dp[nmask] %= MOD
	print(dp[-1]) 
	stdout.close()

if __name__ == "__main__":
	main()

Solution (Python):
from sys import stdin, stdout, setrecursionlimit
from collections import deque, defaultdict, Counter
from heapq import heappush, heappop
from functools import lru_cache
import math

setrecursionlimit(20000)
rl = lambda: stdin.readline()
rll = lambda: stdin.readline().split()
rli = lambda: map(int, stdin.readline().split())
rlf = lambda: map(float, stdin.readline().split())

INF, NINF = float('inf'), float('-inf')

def popcount(x):
    x = x - ((x >> 1) & 0x5555555555555555)
    x = (x & 0x3333333333333333) + ((x >> 2) & 0x3333333333333333)
    x = (x + (x >> 4)) & 0x0f0f0f0f0f0f0f0f
    x = x + (x >> 8) 
    x = x + (x >> 16) 
    x = x + (x >> 32) 
    return x & 0x0000007f

def main():
	MOD = 10**9 + 7
	n = int(rl())
	A = []
	for _ in range(n):
		row = []
		for bit in rli():
			row.append(bit)
		A.append(row)

	dp = [0 for _ in range(1 << n)]
	dp[0] = 1
	for mask in range((1 << n)-1):
		new = popcount(mask)
		for cand in range(n):
			if A[new][cand] and not (mask & (1 << cand)):
				nmask = mask ^ (1 << cand)
				dp[nmask] += (dp[mask] % MOD)
				dp[nmask] %= MOD
	print(dp[-1]) 
	stdout.close()

if __name__ == "__main__":
	main()

Solution (Python):
import itertools

mod=10**9+7
n=int(input())
arr=[list(map(int,input().split())) for _ in range(n)]
dp=[[0]*(2**n) for _ in range(n+1)]
dp[0][0]=1
for i in range(n):
  for bit in range(2**n):
    if dp[i][bit]!=0:
      for j in range(n):
        if bit&(2**j)==0 and arr[i][j]==1:
          dp[i+1][bit^(2**j)]+=dp[i][bit]
          dp[i+1][bit^(2**j)]%=mod
print(dp[n][2**n-1])