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):
def examA():
    N = I()
    H = LI()
    dp = [inf]*(N)
    dp[0] = 0
    for i in range(N-1):
        dp[i+1] = min(dp[i+1],dp[i]+abs(H[i+1]-H[i]))
        if i<N-2:
            dp[i + 2] = min(dp[i + 2], dp[i]+abs(H[i + 2] - H[i]))
    ans = dp[-1]
    print(ans)
    return

def examB():
    N, K = LI()
    H = LI()
    dp = [inf]*N
    dp[0] = 0
    for i in range(N):
        for k in range(1,K+1):
            if i+k>=N:
                break
            if dp[i+k]>dp[i]+abs(H[i+k]-H[i]):
                dp[i+k] = dp[i]+abs(H[i+k]-H[i])
    ans = dp[-1]
    print(ans)
    return

def examC():
    N = I()
    A = [0]*N; B = [0]*N; C = [0]*N
    for i in range(N):
        A[i],B[i],C[i] = LI()
    dp = [[0]*3 for _ in range(N+1)]
    for i in range(N):
        dp[i+1][0] = max(dp[i][1],dp[i][2]) + A[i]
        dp[i + 1][1] = max(dp[i][0], dp[i][2]) + B[i]
        dp[i + 1][2] = max(dp[i][1], dp[i][0]) + C[i]
    ans = max(dp[N])
    print(ans)
    return

def examD():
    ans = 0
    print(ans)
    return

def examE():
    ans = 0
    print(ans)
    return

def examF():
    
    S = SI(); N = len(S)
    T = SI(); M = len(T)
    dp = [[0] * (M + 1) for _ in range(N + 1)]
    for i in range(N):
        for j in range(M):
            if S[i] == T[j]:
                dp[i + 1][j + 1] = max(dp[i + 1][j + 1], dp[i][j] + 1)
            dp[i+1][j+1] = max(dp[i+1][j+1], dp[i+1][j], dp[i][j+1])

    i = N; j = M
    ans = []
    while i > 0 and j > 0:
        if dp[i][j] == dp[i-1][j]:
            i -= 1
        elif dp[i][j] == dp[i][j-1]:
            j -= 1
        else:
            ans.append(S[i-1])
            i -= 1
            j -= 1
    print("".join(ans[::-1]))
    return

def examG():
    
    
    def topological_sort(n,outs,ins):
        q = deque(v1 for v1 in range(n) if ins[v1] == 0)
        res = []
        while q:
            v1 = q.popleft()
            res.append(v1)
            for v2 in outs[v1]:
                ins[v2] -= 1
                if ins[v2] == 0:
                    q.append(v2)
        return res
    N, M = LI()
    V = [[]for _ in range(N)]
    start = [0]*N
    for _ in range(M):
        x,y = LI()
        x -= 1; y -= 1
        V[x].append(y)
        start[y] += 1
    Order = topological_sort(N, V, start)

    dp = [0]*N
    for i in Order:
        for v in V[i]:
            dp[v] = max(dp[v],dp[i]+1)
    ans = max(dp)
    print(ans)
    return

def examH():
    H, W = LI()
    A = [SI() for _ in range(H)]
    dp = [[0]*W for _ in range(H)]
    dp[0][0] = 1
    for i in range(H):
        for j in range(W):
            if A[i][j]=="
                continue
            dp[i][j] %=mod
            if i+1<H and dp[i+1][j]!="
                dp[i+1][j] +=dp[i][j]
            if j+1<W and dp[i][j+1]!="
                dp[i][j+1] +=dp[i][j]
    ans = dp[-1][-1]
    print(ans)
    return

def examI():
    N = I()
    P = LFI()
    dp = [[0]*(N+1) for _ in range(N+1)]
    dp[0][0] = 1
    for i in range(N):
        for j in range(i+1):
            dp[i+1][j+1] += dp[i][j]*P[i]
            dp[i+1][j] += dp[i][j]*(1-P[i])
    ans = 0
    for v in dp[N][(1+N//2):]:
        ans += v
    print(ans)
    return

def examJ():
    N = I()
    A = LI()
    dp = [[[0]*(N+1) for _ in range(N+1)]for _ in range(N+1)]
    S = Counter(A)
    print(S)
    dp[S[1]][S[2]][S[3]] = 0
    que = deque()
    que.append((S[1],S[2],S[3]))
    used = defaultdict(set)
    while(que):
        n1,n2,n3 = que.popleft()
        if used[(n1,n2,n3)]:
            continue
        used[(n1,n2,n3)] = True
        now = n1+n2+n3
        if n1>0:
            dp[n1-1][n2][n3] += dp[n1][n2][n3] + N/now * n1/now
            que.append((n1-1, n2, n3))
        if n2>0:
            dp[n1+1][n2-1][n3] += dp[n1][n2][n3] + N/now * n2/now
            que.append((n1+1,n2-1,n3))
        if n3>0:
            dp[n1][n2+1][n3-1] += dp[n1][n2][n3] + N/now * n3/now
            que.append((n1,n2+1,n3-1))

    ans = dp[0][0][0]
    print(ans)
    return

def examK():
    N, K = LI()
    A = LI()
    dp = [-1]*(K+1)
    for i in range(1,K+1):
        flag = False
        for a in A:
            if a>i:
                break
            if dp[i-a]==-1:
                flag = True
                break
        if flag:
            dp[i] = 1

    ans = dp[K]
    if ans==1:
        print("First")
    else:
        print("Second")
    return

def examM():
    N, K = LI()
    A = LI()
    dp = [[0]*(K+1) for _ in range(N+1)]
    dp[0][0] = 1
    for i in range(N):
        cum = [0]*(K+1)
        cum[0] = dp[i][0]
        for c in range(K):
            cum[c+1] = cum[c] + dp[i][c+1]
        for j in range(K+1):
            if j<=A[i]:
                dp[i+1][j] = cum[j]
            else:
                dp[i+1][j] = cum[j] - cum[j-A[i]-1]
            dp[i+1][j] %=mod

    ans = dp[N][K]%mod
    print(ans)
    return

def examO():
    N = I()
    A = [LI() for _ in range(N)]
    loop = 2**N
    dp = [[0]*loop for _ in range(N+1)]
    dp[0][0] = 1
    for mask in range(loop):
        cur = 0
        for j in range(N):
            if mask&(1<<j)==(1<<j):
                cur +=1
        for i in range(N):
            if mask&(1<<i)==0 and A[i][cur]==1:
                dp[cur+1][mask|(1<<i)] += dp[cur][mask]
                dp[cur + 1][mask | (1 << i)] %=mod

    ans = dp[N][loop-1]
    print(ans)
    return

def examP():
    def dfs_dp(n,s,fords,visited):
        W = 1; B = 1
        visited[s] = True
        for v in fords[s]:
            if visited[v]:
                continue
            cur = dfs_dp(n, v, fords,visited)
            W *= sum(cur)
            B *= cur[0]
            W %= mod
            B %= mod

        return (W,B)
    N = I()
    V = [[]for _ in range(N)]
    for _ in range(N-1):
        x,y = LI()
        x -=1; y -=1
        V[x].append(y)
        V[y].append(x)
    visited = [False]*N
    ans = dfs_dp(N,0,V,visited)
    print(sum(ans)%mod)
    return

import sys,copy,bisect,itertools,heapq,math
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))
def LFI(): return list(map(float,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
alphabet = [chr(ord('a') + i) for i in range(26)]

sys.setrecursionlimit(10**6)

if __name__ == '__main__':
    examO()

"""

"""

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

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

dp = [[0]*(2**N) for i in range(N)]

for k in range(N):
  if a[0][k]==1:
    dp[0][(1<<k)] = 1
    
for i in range(1, N):
  for b in range(2**N):
    if dp[i-1][b]==0:
      continue
    for k in range(N):
      if a[i][k]==1 and b&(1<<k)==0:
        target = (b|(1<<k))
        dp[i][target] += dp[i-1][b]
        if dp[i][target]>DIV:
          dp[i][target] %= DIV
          
s = 0
for i in range(2**N):
  s += dp[N-1][i]
  if s>DIV:
    s %= DIV
    
print(s)

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

def popcnt(x):
  ret = 0
  for i in range(n):
    if x & (1<<i):
      ret += 1
  return ret

for bit in range(1<<n):
  if ans[bit]:
    i = popcnt(bit)
    v = 1
    for j in range(n):
      if not bit & v and a[i][j]:
        ans[bit|v] = (ans[bit|v] + ans[bit])%mod
      v <<= 1
print(ans[-1])

Solution (Python):
def main():

    n = int(input())
    a = [list(map(int, input().split())) for _ in range(n)]

    MOD = 10**9 +7
    dp = [0] * (1<<n)

    pair_dic = {}
    for i in range(n+1):
        pair_dic[i] = []

    for i in range(1<<n):
        pair_dic[ bin(i).count('1') ].append(i)

    dp[0] = 1

    for i in range(n):
        for k in range(n):
            if(a[i-1][k] == 0):
                continue
            for j in pair_dic[i]:
                if(  ((j&(1<<k)) >>k == 0 ) ):
                    dp[j + (1<<k)] += dp[j]
                    dp[j + (1<<k)] = dp[j + (1<<k)] % MOD

    print(dp[-1])

if __name__ == '__main__':
    main()

Solution (Python):
N, = map(int, input().split())
M=10**9+7
n=N
_ALL=1<<n
def itercomb(n,k):
    if k==0:
        yield 0
    else:
        _a = (1<<k)-1
        yield _a
        while 1:
            _x = _a & -_a
            _y = _a + _x
            _a = (((_a & ~_y) // _x) >> 1) | _y
            if _a >= _ALL: break
            yield _a


xs = []
for _ in range(N):
    xs.append([1<<i for i, x in enumerate(input().split()) if int(x)])
d = [0]*(1<<N)
d[0] = 1
for i, y in enumerate(xs):
    for x in itercomb(N, i):
        for ai in y:
            if (not x&ai):
                d[x|ai] = (d[x] + d[x|ai])%M
print(d[-1])

Solution (Python):
N, = map(int, input().split())
M=10**9+7
n=N
_ALL=1<<n
def itercomb(n,k):
    if k==0:
        yield 0
    else:
        _a = (1<<k)-1
        yield _a
        while 1:
            _x = _a & -_a
            _y = _a + _x
            _a = (((_a & ~_y) // _x) >> 1) | _y
            if _a >= _ALL: break
            yield _a


xs = []
for _ in range(N):
    xs.append([1<<i for i, x in enumerate(input().split()) if int(x)])
d = [0]*(1<<N)
d[0] = 1
for i, y in enumerate(xs):
    for x in itercomb(N, i):
        for ai in y:
            if (not x&ai):
                d[x|ai] += d[x]%M
                if d[x|ai] >M:
                    d[x|ai]-=M
print(d[-1])

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

n = int(input())
a = [[int(x) for x in input().split()] for _ in range(n)]

dp = [-1 for __ in range(1 << n)]
dp[0] = 1


def dfs(bit, i):
    if dp[bit] != -1:
        return dp[bit]
    else:
        tmp = 0
        for j in range(n):
            if bit & (1 << j) and a[i][j]:
                tmp += dfs(bit ^ (1 << j), i+1)
    dp[bit] = tmp % mod
    return dp[bit]


print(dfs(2**n-1, 0))

Solution (Python):

import sys
INF = float("inf")

MOD = 1000000007  


def bit_sum(n):
    wa = 0
    while n > 0:
        wa += n & 1
        n >>= 1
    return wa


def bit(n, k):
    return (n >> k) & 1


def solve(N: int, a: "List[List[int]]"):
    dp = [0]*(1 << N)
    dp[0] = 1
    for S in range(1, 1 << N):
        i = bit_sum(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])
    return


def main():

    def iterate_tokens():
        for line in sys.stdin:
            for word in line.split():
                yield word
    tokens = iterate_tokens()
    N = int(next(tokens))  
    a = [[int(next(tokens)) for _ in range(N)]
         for _ in range(N)]  
    solve(N, a)


if __name__ == '__main__':
    main()

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

dp = [[0] * (1 << n) for i in range(n+1)]
dp[0][0] = 1

for i in range(n):
    bit_state = (1 << i) - 1
    if bit_state == 0:
        for j in range(n):
            if not(bit_state & (1 << j)) and a[i][j] == 1:
                dp[i+1][bit_state|(1 << j)] += dp[i][bit_state]
                dp[i+1][bit_state|(1 << j)] %= MOD
        continue
    while bit_state < (1 << n):
        for j in range(n):
            if not(bit_state & (1 << j)) and a[i][j] == 1:
                dp[i+1][bit_state|(1 << j)] += dp[i][bit_state]
                dp[i+1][bit_state|(1 << j)] %= MOD
        x = bit_state & -bit_state
        y = bit_state + x
        bit_state = ((bit_state & ~y) // x >> 1) | y
print(dp[n][(1 << n) - 1])

Solution (Python):
N = int(input())
A = [[int(a) for a in input().split()] for _ in range(N)]

mod = 10**9+7

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

POPCOUNT_TABLE16 = [0] * 2**16
for index in range(len(POPCOUNT_TABLE16)):
    POPCOUNT_TABLE16[index] = (index & 1) + POPCOUNT_TABLE16[index >> 1]

def popcnt(v):
    return (POPCOUNT_TABLE16[ v        & 0xffff] +
            POPCOUNT_TABLE16[(v >> 16) & 0xffff] +
            POPCOUNT_TABLE16[(v >> 32) & 0xffff] +
            POPCOUNT_TABLE16[(v >> 48)         ])


for s in range(1, 2**N):
    i = popcnt(s)
    cnt = 0
    for j in range(N):
        if A[i-1][j] == 1 and s&(1<<j):
            cnt += dp[s^(1<<j)]%mod
    dp[s] = cnt%mod
            
print(dp[-1]%mod)

Solution (Python):
def main():
    N = int(input())
    adj = []
    for i in range(N):
        adj.append(list(map(int, input().split())))
    mod = 10 ** 9 + 7

    state_list = [[] for _ in range(N + 1)]
    for state in range(2 ** N):
        one_num = 0
        for j in range(N):
            if state >> j & 1:
                one_num += 1
        state_list[one_num].append(state)

    dp = [0] * (2**N)
    dp[0] = 1
    for i in range(N):
        for j in range(N):
            if adj[i][j]:
                for state in state_list[i]:
                    state_new = state | (2**j)
                    if state != state_new:
                        dp[state_new] += dp[state]
                        if dp[state_new] >= mod:
                            dp[state_new] -= mod

    print(dp[-1])


if __name__ == '__main__':
    main()

Solution (Python):
def main():
    N = int(input())
    adj = []
    for i in range(N):
        adj.append(list(map(int, input().split())))
    mod = 10 ** 9 + 7

    state_list = [[] for _ in range(N + 1)]
    for state in range(2 ** N):
        one_num = 0
        for j in range(N):
            if state >> j & 1:
                one_num += 1
        state_list[one_num].append(state)

    dp = [0] * (2**N)
    dp[0] = 1
    for i in range(N):
        for j in range(N):
            if adj[i][j]:
                for state in state_list[i]:
                    state_new = state | (2**j)
                    if state != state_new:
                        dp[state_new] = (dp[state_new] + dp[state]) % mod

    print(dp[-1])


if __name__ == '__main__':
    main()

Solution (Python):
MOD = int(1e9 + 7)

n = int(input())
can = []
for _ in range(n):
    can.append([int(c) for c in input().split()])

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

for mask in range(2 ** n - 1):
    a = bin(mask).count('1')
    for b in range(n):
        if can[a][b] and not(mask & (1 << b)):
            m2 = mask ^ (1 << b)
            dp[m2] = dp[m2] + dp[mask]
            if dp[m2] > MOD:
                dp[m2] -= MOD


print(dp[(2**n) - 1])

Solution (Python):
def bitcount(i): 
    i = i - ((i >> 1) & 0x55555555)
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
    return (((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) & 0xffffffff) >> 24

MOD = 10**9+7
n = int(input())
can = [[int(i) for i in input().split()] for j in range(n)]
dp = [0] * (1 << n)
dp[0] = 1
for mask in range((1 << n) - 1):
    a = bin(mask).count('1')
    
    for b in range(n):
        if can[a][b] and not (mask & (1 << b)):
            mask2 = mask ^ (1 << b)
            dp[mask2] += dp[mask]
            if dp[mask2] >= MOD:
                dp[mask2] -= MOD
print(dp[(1 << n) - 1])

Solution (Python):
def bitcount(i): 
    i = i - ((i >> 1) & 0x55555555)
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
    return (((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) & 0xffffffff) >> 24

MOD = 10**9+7
n = int(input())
can = [[int(i) for i in input().split()] for j in range(n)]
dp = [0] * (1 << n)
dp[0] = 1
for mask in range((1 << n) - 1):
    
    a = bitcount(mask)
    for b in range(n):
        if can[a][b] and not (mask & (1 << b)):
            mask2 = mask ^ (1 << b)
            dp[mask2] += dp[mask]
            if dp[mask2] >= MOD:
                dp[mask2] -= MOD
print(dp[(1 << n) - 1])

Solution (Python):
import sys

sys.setrecursionlimit(10 ** 6)
int1 = lambda x: int(x) - 1
p2D = lambda x: print(*x, sep="\n")
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(rows_number): return [list(map(int, sys.stdin.readline().split())) for _ in range(rows_number)]

def main():
    md=10**9+7
    n=int(input())
    aa=LLI(n)
    lov=[[] for _ in range(n)]
    for i,ar in enumerate(aa):
        for j,a in enumerate(ar):
            if a:lov[i].append(j)
    dp=[0]*(1<<n)
    dp[0]=1
    for s in range((1<<n)-1):
        pre=dp[s]
        if pre==0:continue
        i=bin(s).count("1")
        for j in lov[i]:
            if s>>j&1:continue
            ns=s|1<<j
            dp[ns]=(dp[ns]+pre)%md
    print(dp[-1])
    

main()

Solution (Python):
def main():
    mod = 10**9+7
    N = int(input())
    a = [list(map(int, input().split())) for _ in range(N)]
    dp = [0]*(1<<N)
    dp[0] = 1
    def bitcount(i): 
        i = i - ((i >> 1) & 0x55555555)
        i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
        return (((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) & 0xffffffff) >> 24
    for s in range(1, 1<<N):
        women = []
        m = bitcount(s)-1
        for w in range(N):
            if a[m][w] == 1 and s&1<<w:
                dp[s]+=dp[s^(1<<w)]
                dp[s]%=mod
    print(dp[-1])
main()

Solution (Python):
mod = 10**9+7
N = int(input())
a = [list(map(int, input().split())) for _ in range(N)]
dp = [0]*(1<<N)
dp[0] = 1
def bitcount(i): 
    i = i - ((i >> 1) & 0x55555555)
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333)
    return (((i + (i >> 4) & 0xF0F0F0F) * 0x1010101) & 0xffffffff) >> 24
for s in range(1, 1<<N):
    women = []
    m = bitcount(s)-1
    for w in range(N):
        if a[m][w] == 1 and s&1<<w:
            dp[s]+=dp[s^(1<<w)]
            dp[s]%=mod
print(dp[-1])

Solution (Python):
mod = 10**9+7
N = int(input())
a = [list(map(int, input().split())) for _ in range(N)]
dp = [0]*(1<<N)
dp[0] = 1
for s in range(1, 1<<N):
    women = []
    for i in range(N):
        if (1<<i)&s:
            women.append(i)
    m = len(women)-1
    for w in women:
        if a[m][w] == 1:
            dp[s]+=dp[s^(1<<w)]
            dp[s]%=mod
print(dp[-1])

Solution (Python):
N = int(input())
A = [[int(_) for _ in input().split()] for i in range(N)]

dp = [[0] * ((1 << 21) + 10) for _ in range(N + 1)]
dp[0][0] = 1
MOD = 10**9+7
for i in range(N):
    
    for s in range(1 << N):  
        if dp[i][s] == 0: continue

        for j in range(N):
            jb = 1 << j
            if not (s & jb) and A[i][j]:
                
                dp[i + 1][s | jb] += dp[i][s]
print(dp[N][(1 << N) - 1] % MOD)

Solution (Python):
N = int(input())
A = [[int(_) for _ in input().split()] for i in range(N)]

dp = [[0] * ((1 << 21) + 10) for _ in range(N + 1)]
dp[0][0] = 1
MOD = 10**9+7
for i in range(N):
    
    for s in range(1 << N):  
        if dp[i][s] == 0: continue

        for j in range(N):
            jb = 1 << j
            if not (s & jb) and A[i][j]:
                
                dp[i + 1][s | jb] += dp[i][s] % MOD
print(dp[N][(1 << N) - 1] % MOD)

Solution (Python):
def solve(n, a):
    mod = 10 ** 9 + 7

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

    for mask in range(1, 1 << n):
        man = bin(mask).count('1') - 1
        for woman in range(n):
            if (mask >> woman) & 1 == 1 and a[man][woman] == 1:
                dp[mask] = (dp[mask] + dp[mask & ~(1 << woman)]) % mod

    return dp[(1 << n) - 1]


if __name__ == '__main__':
    n = int(input())
    a = tuple(tuple(map(int, input().split())) for _ in range(n))
    print(solve(n, a))

Solution (Python):
import sys
MAX_INT = int(10e15)
MIN_INT = -MAX_INT
mod = 1000000007
sys.setrecursionlimit(1000000)
def IL(): return list(map(int,input().split()))
def SL(): return input().split()
def I(): return int(sys.stdin.readline())
def S(): return input()

N = I()
a = [IL() for _ in range(N)]

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


for i in range(2, N+1): 
  n = (1 << i) - 1
  while n < (1 << N):
    x = n
    
    for j in range(N):
      if (x >> j) & 1:
        if a[i - 1][j] == 1:
          dp[x] += dp[x ^ (1 << j)] %mod

    """
    while x != 0:
      if x & 1 == 1:
        if a[i - 1][cnt] == 1:
          DP[n] = (DP[n] + DP[n & ~(1 << (cnt))]) % MOD
      x >>= 1
      cnt += 1
    """
    
    x = n & -n
    y = n + x
    n = ((n & ~y) // x >> 1) | y

print(dp[-1] %mod)

Solution (Python):
def popcount(x):
    c = 0
    while x > 0:
        if x & 1:
            c += 1
        x = x >> 1
    return c

def f(A):
    MOD = (10 ** 9) + 7
    n = len(A)

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

    for mask in range((1 << n) - 1):
        a = popcount(mask)
        for b in range(n):
            if A[a][b] and not (mask & (1 << b)):
                m2 = mask ^ (1 << b)
                dp[m2] += dp[mask]
                if dp[m2] > MOD:
                    dp[m2] -= MOD

    return dp[(1 << n) - 1]

assert f([[0, 1, 1], [1, 0, 1], [1, 1, 1]]) == 3

n = int(input())
A = []
for _ in range(n):
    R = list(map(int, input().split()))
    A.append(R)

ans = f(A)
print(ans)

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

def main():
  N = int(input())
  A = [tuple(map(int,input().split())) for _ in range(N)]
  MOD =10**9+7
  N2 = 1<<N

  dp = [1] + [0 for _ in range(N2-1)]

  for i in range(N):
    nextdp = [0 for _ in range(N2)]
    for j in range(N2):
      if dp[j] != 0:
        for k in range(N):
          if (j>>k)&1 == 0 and A[i][k] == 1:
            idx = j|(1<<k)
            nextdp[idx] += dp[j]
            
    dp = nextdp

  print(dp[N2-1]%MOD)

if __name__ == '__main__':
    main()

Solution (Python):
N = int(input())

dp = [0] * (1 << N)
dp[0] = 1
curr = {0}
MOD = 10**9+7

for _ in range(N):
    a = [1 << i for i, v in enumerate(input().split()) if v == '1']
    nxt = set()
    for i in curr:
        dpi = dp[i] % MOD
        for j in a:
            if i & j:
                continue
            dp[i | j] += dpi
            nxt.add(i | j)
    curr = nxt
print(dp[(1 << N) - 1] % MOD)

Solution (Python):



import sys
import math
from bisect import bisect_right as br
from bisect import bisect_left as bl
sys.setrecursionlimit(1000000000)
from heapq import heappush, heappop,heappushpop
from collections import defaultdict
from itertools import accumulate
from collections import Counter
from collections import deque
from operator import itemgetter
from itertools import permutations
mod = 10**9 + 7
inf = float('inf')
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))

n = I()
a = [LI() for _ in range(n)]
dp = [[0]*(2**n) for _ in range(n+1)]
dp[0][0] = 1
for i in range(n):
    for j in range(2**n):
        if dp[i][j]:
            for k in range(n):
                if (j>>k)&1 == 0 and a[i][k] == 1:
                    dp[i+1][j|(1<<k)] += dp[i][j]
                    dp[i+1][j|(1<<k)] %= mod
print(dp[-1][-1])

Solution (Python):
import subprocess as sp
import sys
 
code = r"""

using namespace std;
const int mod=1e9+7;
int add(int x,int y){return(x+=y)<mod?x:x-mod;}
int sub(int x,int y){return(x-=y)>= 0?x:x+mod;}
int dp[1<<21];
int n,a[21][21];
int main(){
  scanf("%d",&n);
  for(int i=0;i<n;++i){
    for(int j=0;j<n;++j){
      scanf("%d",&a[i][j]);
    }
  }
  dp[0]=1;
  for(int b=1;b<(1<<n);++b){
    int bc=__builtin_popcount(b)-1;
    for(int i=0;i<n;++i){
      if((b>>i&1)&&a[bc][i]){
        dp[b]=add(dp[b],dp[b^(1<<i)]);
      }
    }
  }
  printf("%d\n",dp[(1<<n)-1]);
}
"""
 
with open("A.cpp","w") as f:
	f.write(code)
    
sp.Popen(['g++','-std=c++11','-O2',"A.cpp"]).communicate()
sp.Popen(['./a.out'],stdin=sys.stdin,stdout=sys.stdout).communicate()

Solution (Python):
mod=1000000007


n=int(input())
dp=[0] * (1<<n)
dp[0]=1
curr={0}
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):
import numpy as np

MOD = 10 ** 9 + 7
N=int(input())
A = [[int(x) for x in input().split()] for _ in range(N)]

dp = np.zeros(1<<N, dtype=np.uint64)
dp[0] = 1

for A_row in A:
    prev = dp.copy()
    dp = np.zeros_like(prev)
    for n,a in enumerate(A_row):
        if not a:
            continue
        x = 1 << n
        dp[x:] += prev[:-x]
    dp %= MOD

answer = dp[-1]
print(answer)

Solution (Python):

from collections import defaultdict
from collections import deque
from heapq import heappush, heappop
import sys
import math
import bisect
import random
def LI(): return list(map(int, sys.stdin.readline().split()))
def I(): return int(sys.stdin.readline())
def LS():return list(map(list, sys.stdin.readline().split()))
def S(): return list(sys.stdin.readline())[:-1]
def IR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = I()
    return l
def LIR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = LI()
    return l
def SR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = S()
    return l
def LSR(n):
    l = [None for i in range(n)]
    for i in range(n):l[i] = SR()
    return l
sys.setrecursionlimit(1000000)
mod = 1000000007



def A():
    n,k = LI()
    a = LI()
    dp = [0 for i in range(k+1)]
    f = [0 for i in range(k+1)]
    for i in range(k+1):
        if not f[i]:
            for j in a:
                if i-j >= 0:
                    dp[i] = 1-dp[i-j]
                    if dp[i] == 0:
                        break
            else:
                dp[i] = 1
    print(["First","Second"][dp[k]])

def B():
    n = I()
    a = LI()
    ma = n*(n+1)//2
    su = sum(a)
    s = [1 for i in range(n+1)]
    for i in range(1,n):
        s[i+1] = s[i]+n-i+1
    d = [None for i in range(ma+1)]
    d[0] = 0
    k = 1
    for i in range(1,n+1):
        for j in range(n+1-i):
            d[k+j] = i
        k += n+1-i
    k = n%2
    dp = [su if (d[i]+k)%2 else -su for i in range(ma+1)]
    dp[0] = 0
    if k:
        for i in range(n):
            dp[i+1] = a[n-1-i]
    else:
        for i in range(n):
            dp[i+1] = -a[n-1-i]
    for i in range(1,ma):
        r = n-1-(i-s[d[i]])
        l = r+1-d[i]
        j = i+(n-d[i])
        j2 = j+1
        if (d[i]+k)%2:
            if r < n-1:
                if a[r+1]+dp[i] > dp[j]:
                    dp[j] = a[r+1]+dp[i]
            if l > 0:
                if a[l-1]+dp[i] > dp[j2]:
                    dp[j2] = a[l-1]+dp[i]
        else:
            if r < n-1:
                if -a[r+1]+dp[i] < dp[j]:
                    dp[j] = -a[r+1]+dp[i]
            if l > 0:
                if -a[l-1]+dp[i] < dp[j2]:
                    dp[j2] = -a[l-1]+dp[i]
    print(dp[ma])
    return


def C():
    def comb(a,b):
        return fact[a]*inv[b]*inv[a-b]%mod

    n,k = LI()
    a = LI()
    """
    fact = [1]
    for i in range(n+k):
        fact.append(fact[-1]*(i+1)%mod)
    inv = [1]*(n+k+1)
    inv[n+k] = pow(fact[n+k],mod-2,mod)
    for i in range(n+k)[::-1]:
        inv[i] = inv[i+1]*(i+1)%mod
    f = [comb(i+n-1,n-1) for i in range(k+1)]
    for i in a:
        for j in range(k-i)[::-1]:
            f[j+i+1] -= f[j]
            f[j+i+1] %= mod
    print(f[k])
    """
    dp = [[0]*(k+1) for i in range(n+1)]
    for i in range(k+1):
        dp[0][i] = 1
    for i in range(n):
        ni = i+1
        ai = a[i]
        for j in range(k+1):
            if j >= ai+1:
                dp[ni][j] = (dp[i][j]-dp[i][j-ai-1])%mod
            else:
                dp[ni][j] = dp[i][j]
        if i < n-1:
            for j in range(k):
                dp[ni][j+1] += dp[ni][j]
                dp[ni][j+1] %= mod
    print(dp[n][k]%mod)
    return


def D():
    n = I()
    a = LI()
    a.insert(0,0)
    for i in range(n):
        a[i+1] += a[i]
    dp = [[0]*(n+1) for i in range(n)]
    for le in range(1,n):
        for l in range(n-le):
            r = l+le+1
            res = float("inf")
            for i in range(l+1,r):
                m = dp[l][i]+dp[i][r]
                if m < res:
                    res = m
            res += a[r]-a[l]
            dp[l][r] = res
    print(dp[0][n])
    return


def E():

    n = I()
    a = LIR(n)
    m = 1<<n
    dp = [0]*m
    dp[0] = 1
    a = [[j for j in range(n) if a[i][j]] for i in range(n)]
    bit_count = [-1]
    for i in range(n):
        for j in range(len(bit_count)):
            bit_count.append(bit_count[j]+1)
    for k in range(1,m):
        i = bit_count[k]
        for j in a[i]:
            bj = 1<<j
            if bj > k:
                break
            if not k&bj:
                continue
            nk = k^bj
            dp[k] += dp[nk]
    print(dp[m-1]%mod)

def F():
    return


def G():
    return


def H():
    return


if __name__ == "__main__":
    E()