Bingo

time limit per test: 7 seconds

memory limit per test: 512 megabytes

Articles

Getting ready for VK Fest 2021, you prepared a table with $n$ rows and nn columns, and filled each cell of this table with some event related with the festival that could either happen or not: for example, whether you will win a prize on the festival, or whether it will rain.

Forecasting algorithms used in VK have already estimated the probability for each event to happen. Event in row ii and column jj will happen with probability $a_{i,j}⋅10^{−4}$ . All of the events are mutually independent.

Let’s call the table winning if there exists a line such that all $n$ events on it happen. The line could be any horizontal line (cells $(i,1),(i,2),…,(i,n)$ for some ii), any vertical line (cells $(1,j),(2,j),…,(n,j)$ for some $j$ ), the main diagonal (cells $(1,1),(2,2),…,(n,n)$ ), or the antidiagonal (cells $(1,n),(2,n−1),…,(n,1)$ ).

Find the probability of your table to be winning, and output it modulo $31607$ (see Output section).

Input

The first line contains a single integer $n (2≤n≤21)$ — the dimensions of the table.

The $i$ -th of the next $n$ lines contains $n$ integers $ai,1,ai,2,…,ai,n$ ( $0<a_{i,j}<10^4$ ). The probability of event in cell $(i,j)$ to happen is $a_{i,j}⋅10^{−4}$ .

Output

Print the probability that your table will be winning, modulo $31607$ .

Formally, let $M=31607$ . It can be shown that the answer can be expressed as an irreducible fraction $\frac{p}{q}$ , where $p$ and $q$ are integers and $q\equiv̸0(modM)$ . Output the integer equal to $p⋅q^{−1} mod \ M$ . In other words, output such an integer $x$ that $0≤x<M$ and $x⋅q≡p(mod\ M)$ .

Examples

input

1
2
3

2
5000 5000
5000 5000

output

input

1
2
3

2
2500 6000
3000 4000

output

input

3
1000 2000 3000
4000 5000 6000
7000 8000 9000

output

Tutorial

Let $l1,l2,…,l2n+2$ denote the $2n+2$ possible lines that can be formed. Let $Li$ denote the event that line $li$ is formed, and $\overline{Li}$ denote the event that line $li$ is not formed ( $i.e.$ , $P(Li)+P(\overline{Li})=1$ ).

Let’s find the probability that our table is not winning. It is equal to $P(\overline{L_1}∩\overline{L_2}∩…∩\overline{L_{2n+2}})$ .

Note that the following two statements are true:

$P(\overline{L_1}∩\overline{L_2}∩…∩\overline{L_{2n+2}})+P({L_1}∩\overline{L_2}∩…∩\overline{L_{2n+2}})=P(\overline{L_2}∩…∩\overline{L_{2n+2}})$ ;
$P(L_1∩\overline{L_2}∩…∩\overline{L_{2n+2}})=P(\overline{L_2}∩…∩\overline{L_{2n+2}}|L_1)⋅P(L_1).$

The first one follows from the law of total probability, and the second one follows from the definition of conditional probability.

These two statements combined allow us to use the following formula:

$P(\overline{L_1}∩\overline{L_2}∩…∩\overline{L_{2n+2}})=P(\overline{L_2}∩…∩\overline{L_{2n+2}})−P(\overline{L_2}∩…∩\overline{L_{2n+2}}|L_1)⋅P(L_1).$

We can apply this formula recursively. Specifically, we can make a function $f(i,S)$ , where $S={s_1,s_2,…,s_k}$ is a subset of ${1,2,…,i−1}$ , which calculates $P(L_i∩L_{i+1}∩…∩L_{2n+2}|L_{s1}∩…∩L_{sk})$ . For $i=2n+3$ , $f(i,S)=1$ , and for $i≤2n+2$ we can generalize the formula above as follows:

$f(i,S)=f(i+1,S)−f(i+1,S∪{i})⋅P(L_i|L_{s1}∩…∩L_{sk})$ .

Here, $P(L_i|L_{s1}∩…∩L_{sk})$ is the probability that line $li$ is formed given that lines $l_{s1},…,l_{sk}$ are formed. This is equal to the product of probabilities of all cells belonging to $li$ which do not belong to any of $l_{s1},…,l_{sk}$ .

The answer to the problem, $i.e.$ the probability that our table is winning, is $1−f(1,\{\})$ .

This allows us to implement an $O(2^n⋅n)$ solution, which is too slow. In fact, this solution is equivalent to applying inclusion-exclusion principle.

To optimize this solution, note that once it becomes easy to calculate $f(i,S)$ , we don’t have to make any more recursive calls. Why would it become easy to calculate $f(i,S)$ though?

Let’s order the lines in such a way that $l_{n+3},l_{n+4},…,l_{2n+2}$ are the horizontal lines of the table. Consider a call of the form $f(n+3,S)$ . This call is basically asking: “what is the probability that none of lines $l_{n+3},l_{n+4},…,l_{2n+2}$ are formed, given that lines $l_{s1},…,L_{sk}$ are formed?”.

Note that the horizontal lines are independent, and we can actually answer this question in $O(n^2)$ . Specifically, for any horizontal line, the probability that it is not formed is $1$ minus the product of probabilities of all its cells not belonging to any of $l_{s1},…,l_{sk}$ . The overall value of $f(n+3,S)$ is the product of probabilities for individual horizontal lines.

This way, we have built an $O(2^n⋅n^2)$ solution. This might be fast enough depending on your implementation, but there are at least two ways to optimize it to $O(2^n⋅n)$ :

The first way is to maintain the products of probabilities of untouched cells for horizontal lines on the fly. For simplicity, assume that $l_1$ and $l_2$ are the diagonal lines. For each $i=3,4,…,n+2$ , after we process (vertical) line $l_i$ , we can update the products for all horizontal lines with the cells of $l_i$ (in $O(n)$ , make a recursive call, and roll back the updates (in $O(n)$ again). Once we get to $i=n+3$ , instead of going through every cell of the table in $O(n^2)$ , we can just multiply $n$ horizontal line products in $O(n)$ .
The second way is to define $g(i,mask)$ $(1≤i≤n; mask$ is a bitmask of size $n$ ) to be the product of $a_{i,j}$ over all $j$ belonging to $mask$ . All values of $g$ can be calculated in $O(2^n⋅n)$ using dynamic programming: $g(i,0)=1$ , and $g(i,mask)=g(i,mask⊕2^j)⋅a_{i,j}$ , where jj is any bit set in $mask$ . When we arrive at a $f(n+3,S)$ call, for each of the $n$ horizontal lines, instead of going through its cells, we can construct the mask of cells the values in which we want to multiply, and use the corresponding value of $g$ in $O(1)$ .

code

#include<bits/stdc++.h>
using namespace std;

const int mod = 31607;
const int inv = 3973;

int n, a[22][22], g[22][1<<21], dp[1<<23];

void D(int x, int y, int S) {
    if (S >> y & 1) g[x][S] = g[x][S^(1<<y)] * a[x][y+1] % mod;
    if (y+1 < n) D(x, y+1, S|(1<<y+1)), D(x, y+1, S);
}

int dfs(int x, int S) {
    if (x == n+3) {
        return dp[S];
    }
    int k = 1;
    if (x <= n) {
        for(int i = 1; i <= n; ++ i) {
            k = k * a[i][x] % mod;
        }
    }
    else {
        if (x == n+1) {
            for(int i = 0; i < n; ++ i) {
                if (S>>i&1) continue;
                k = k * a[i+1][i+1] % mod;
            }
        }
        else {
            for(int i = 0; i < n; ++ i) {
                if (S>>(n-1-i)&1) continue;
                if ((S>>n&1) && (n&1) && i == n/2) continue;
                k = k * a[i+1][n-i] % mod;
            }
        }
    }
    return (mod + (dfs(x+1, S) - k * dfs(x+1, S|(1ll<<(x-1))) % mod) % mod) % mod;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cin >> n;
    for(int i = 1; i <= n; ++ i) {
        for(int j = 1; j <= n; ++ j) {
            cin >> a[i][j];
            a[i][j] = a[i][j] * inv % mod;
        }
    }
    for(int i = 1; i <= n; ++ i) {
        g[i][0] = 1;
        D(i, 0, 1), D(i, 0, 0);
    }
    for(int S = 0; S < (1<<(n+2)); ++ S) {
        int res = 1;
        for(int i = 1; i <= n; ++ i) {
            int k = S;
            if (k >> n & 1) {
                k ^= 1<<n;
                k |= 1<<(i-1);
            }
            if (k >> (n+1) & 1) {
                k ^= 1<<(n+1);
                k |= 1<<(n-i);
            }
            res = res * (mod + 1 - g[i][((1<<n)-1)^k]) % mod;
        }
        dp[S] = res;
    }
    cout << (mod+1- dfs(1, 0)) % mod << "\n";
    return 0;
}