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):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools

sys.setrecursionlimit(10**7)
inf = 10**20
eps = 1.0 / 10**10
mod = 10**9+7
dd = [(-1,0),(0,1),(1,0),(0,-1)]
ddn = [(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1)]

def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI_(): return [int(x)-1 for x in sys.stdin.readline().split()]
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline())
def F(): return float(sys.stdin.readline())
def S(): return input()
def pf(s): return print(s, flush=True)


def main():
    n = I()
    a = [LI() for _ in range(n)]
    m = 2**n
    s = [[0] * m for _ in range(n+1)]
    s[0][0] = 1
    ii = [2**i for i in range(n)]
    for i in range(n):
        t = s[i+1]
        u = s[i]
        b = [ii[c] for c in range(n) if a[i][c]]
        for k in range(m):
            v = u[k]
            v = v % mod
            if v < 1:
                continue
            for j in b:
                if not (j & k):
                    t[k | j] += v

    return s[-1][-1] % mod


print(main())

Solution (Python):
n = int(input())
dp = [0] * (1 << n)
dp[0] = 1
curr = {0}
MOD = 10 ** 9 + 7
for i in range(n):
    compatibilities = [1 << i for i, v in enumerate(input().split()) if v == '1']
    nxt = set()
    for k in curr:
        p = dp[k] % MOD
        for c in compatibilities:
            if k & c:
                continue
            dp[k | c] += p
            nxt.add(k | c)
    curr = nxt
print(dp[(1 << n) - 1] % MOD)

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

I = [[j for j in range(N) if Ai[j]] for Ai in A]

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

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

cs = [0]*(1 << N)
used = [0]*(1 << N)
s = (1 << N) - 1
que = deque([(s, 0)])
cs[s] = 1
used[s] = 1

while que:
    state, c = que.popleft()
    v = cs[state] % MOD

    for i in range(N):
        if state & (1 << i) and A[c][i]:
            n_state = state ^ (1 << i)
            cs[n_state] += v
            if not used[n_state] and c+1 < N:
                que.append((n_state, c+1))
                used[n_state] = 1
print(cs[0] % MOD)

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):
N = int(input())
A = [list(map(int, input().split())) for i in range(N)]

ALL = 1 << N
S = [0]*ALL
MOD = 10**9 + 7
S[ALL-1] = 1
for k in range(N, 0, -1):
    I = [i for i in range(N) if A[k-1][i]]
    v = (1 << k) - 1
    while v < ALL:
        w = S[v] % MOD
        for i in I:
            if v & (1 << i):
                S[v ^ (1 << i)] += w
        x = v & -v; y = v + x
        v = ((v & ~y) // x >> 1) | y
print(S[0] % MOD)

Solution (Python):
n=int(input())
A=[int(input().replace(' ',''),2)for _ in[0]*n]
p=[0]*2**n
p[0]=1
for I in range(2**n-1):
 a=A[bin(I).count('1')]&~I
 for j in range(n):
  if a&1<<j:p[I|1<<j]+=p[I]
print(p[-1]%(10**9+7))

Solution (Python):
n=int(input())
A=[int(input().replace(' ',''),2)for _ in[0]*n]
p=[0]*2**n
p[0]=1
for I in range(2**n-1):
 a=A[bin(I).count('1')]&~I
 for j in range(n):
  if a&1<<j:p[I|1<<j]+=p[I];p[I|1<<j]%=10**9+7
print(p[-1])

Solution (Python):
n=int(input());
A=[int(input().replace(' ',''),2)for _ in[0]*n]
p=[0]*2**n
p[0]=1
for I in range(2**n-1):
 a=A[bin(I).count('1')]&~I
 for j in range(n):
  if a&1<<j:p[I|1<<j]+=p[I];p[I|1<<j]%=10**9+7
print(p[-1])

Solution (Python):
n=int(input());
A=[int(input().replace(' ',''),2)for _ in[0]*n]
dp=[0]*2**n
dp[0]=1
for I in range(2**n-1):
    a=A[bin(I).count('1')]&~I
    for j in range(n):
        if a&1<<j:dp[I|1<<j]=(dp[I|1<<j]+dp[I])%(10**9+7)
print(dp[-1])

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

N = int(input())
As = [int(input().replace(' ', ''), 2) for _ in range(N)]

patWsList = [[] for _ in range(N+1)]
for patW in range(2**N):
    num1 = sum([(patW >> i) & 1 for i in range(N)])
    patWsList[num1].append(patW)

dp = [0] * (2**N)
dp[-1] = 1
for iM, (A, patWs) in enumerate(zip(As, reversed(patWsList))):
    for patW in patWs:
        A2 = A & patW
        dpNow = dp[patW]
        for iW in range(N):
            if A2 & (1 << iW):
                dp[patW ^ (1 << iW)] += dpNow
                dp[patW ^ (1 << iW)] %= MOD

print(dp[0])