Problem 626: Counting Binary Matrices

Project Euler Problem 626 Solution

EulerSolve provides an optimized solution for Project Euler Problem 626, Counting Binary Matrices, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $E_n$ be the set of binary $n\times n$ matrices whose row sums and column sums are all even: $$E_n=\left\{A=(a_{ij})\in\{0,1\}^{n\times n}:\sum_{j=1}^n a_{ij}\equiv 0\pmod 2,\ \sum_{i=1}^n a_{ij}\equiv 0\pmod 2\right\}.$$ Two matrices are considered equivalent when one can be obtained from the other by permuting rows and permuting columns. The quantity $c_n$ is the number of equivalence classes of $E_n$. The problem asks for $c_{20}\bmod 1001001011$. Mathematical Approach The key idea is to count orbits of matrices with even row and column parity under the group $S_n\times S_n$. Burnside's lemma reduces the problem to fixed matrices for a pair of permutations, and those fixed matrices can be described entirely from the cycle lengths of the two permutations. Step 1: Apply Burnside's Lemma The symmetric group $S_n$ acts on rows, another copy acts on columns, and the combined action is $$A\mapsto P_\sigma A P_\tau^{-1},\qquad (\sigma,\tau)\in S_n\times S_n.$$ Therefore $$c_n=\frac{1}{(n!)^2}\sum_{\sigma\in S_n}\sum_{\tau\in S_n}\operatorname{Fix}(\sigma,\tau),$$ where $\operatorname{Fix}(\sigma,\tau)$ is the number of matrices in $E_n$ unchanged by simultaneously permuting rows by $\sigma$ and columns by $\tau$....

Detailed mathematical approach

Problem Summary

Let $E_n$ be the set of binary $n\times n$ matrices whose row sums and column sums are all even:

$$E_n=\left\{A=(a_{ij})\in\{0,1\}^{n\times n}:\sum_{j=1}^n a_{ij}\equiv 0\pmod 2,\ \sum_{i=1}^n a_{ij}\equiv 0\pmod 2\right\}.$$

Two matrices are considered equivalent when one can be obtained from the other by permuting rows and permuting columns. The quantity $c_n$ is the number of equivalence classes of $E_n$. The problem asks for $c_{20}\bmod 1001001011$.

Mathematical Approach

The key idea is to count orbits of matrices with even row and column parity under the group $S_n\times S_n$. Burnside's lemma reduces the problem to fixed matrices for a pair of permutations, and those fixed matrices can be described entirely from the cycle lengths of the two permutations.

Step 1: Apply Burnside's Lemma

The symmetric group $S_n$ acts on rows, another copy acts on columns, and the combined action is

$$A\mapsto P_\sigma A P_\tau^{-1},\qquad (\sigma,\tau)\in S_n\times S_n.$$

Therefore

$$c_n=\frac{1}{(n!)^2}\sum_{\sigma\in S_n}\sum_{\tau\in S_n}\operatorname{Fix}(\sigma,\tau),$$

where $\operatorname{Fix}(\sigma,\tau)$ is the number of matrices in $E_n$ unchanged by simultaneously permuting rows by $\sigma$ and columns by $\tau$.

Step 2: Replace Permutations by Cycle Types

If $\lambda\vdash n$ is a partition of $n$, write $m_k(\lambda)$ for the number of cycles of length $k$. The number of permutations with cycle type $\lambda$ is

$$N(\lambda)=\frac{n!}{\prod_{k\ge 1} k^{m_k(\lambda)}\,m_k(\lambda)!}.$$

So Burnside's sum depends only on pairs of partitions $(\lambda,\mu)$, not on individual permutations. If the row-cycle lengths are $p_1,\dots,p_u$ and the column-cycle lengths are $q_1,\dots,q_v$, then every quantity in the fixed-count calculation can be expressed in terms of these lengths.

Step 3: Count Cell Orbits for One Type Pair

Consider one row cycle of length $p_i$ and one column cycle of length $q_j$. Their Cartesian block contains $p_iq_j$ cells, and the joint action moves each cell by one step along both cycles. That block splits into

$$d_{ij}=\gcd(p_i,q_j)$$

cell orbits. Hence the total number of binary orbit variables is

$$o(\lambda,\mu)=\sum_{i=1}^{u}\sum_{j=1}^{v}\gcd(p_i,q_j).$$

If there were no parity conditions, the fixed matrices for $(\lambda,\mu)$ would simply number $2^{o(\lambda,\mu)}$.

Step 4: Translate Even Row and Column Sums into $\mathbb{F}_2$ Constraints

Inside one $(p_i,q_j)$-block, let $x_{ij,1},\dots,x_{ij,d_{ij}}\in\mathbb{F}_2$ be the orbit bits. A single row in the row cycle sees each orbit bit exactly $q_j/d_{ij}$ times, so modulo $2$ that block contributes to the row-parity equations if and only if

$$\frac{q_j}{d_{ij}}\equiv 1\pmod 2 \iff \nu_2(q_j)\le \nu_2(p_i).$$

Similarly, a single column in the column cycle sees each orbit bit exactly $p_i/d_{ij}$ times, so the block contributes to the column-parity equations if and only if

$$\frac{p_i}{d_{ij}}\equiv 1\pmod 2 \iff \nu_2(p_i)\le \nu_2(q_j).$$

Only the block sum

$$s_{ij}=x_{ij,1}+\cdots+x_{ij,d_{ij}}\pmod 2$$

enters the parity equations. The remaining $d_{ij}-1$ directions inside that block are always free. This is why the final number of free bits takes the form $o(\lambda,\mu)-\rho(\lambda,\mu)$ for a suitable rank $\rho$.

Step 5: Compute the Parity Rank from 2-Adic Valuations

Now work at the level of cycles. Introduce one row-equation node for each row cycle and one column-equation node for each column cycle. The cancellation rules for the block sums are:

$$\nu_2(p_i)=\nu_2(q_j)\Rightarrow \text{the row node and column node must agree},$$

$$\nu_2(p_i)\gt \nu_2(q_j)\Rightarrow \text{the row node is forced to }0,$$

$$\nu_2(p_i)\lt \nu_2(q_j)\Rightarrow \text{the column node is forced to }0.$$

Merging equal nodes and marking forced-zero nodes gives a small graph problem. Let $f(\lambda,\mu)$ be the number of connected components that are not forced to zero. Since there are $u+v$ equation nodes in total, the rank of the parity system is

$$\rho(\lambda,\mu)=u+v-f(\lambda,\mu).$$

Therefore

$$\operatorname{Fix}(\lambda,\mu)=2^{\,o(\lambda,\mu)-\rho(\lambda,\mu)}.$$

Step 6: Final Burnside Sum

Substituting the type weights and the fixed counts gives the exact formula used by the implementations:

$$\boxed{c_n=\frac{1}{(n!)^2}\sum_{\lambda\vdash n}\sum_{\mu\vdash n}N(\lambda)N(\mu)\,2^{\,o(\lambda,\mu)-\rho(\lambda,\mu)}.}$$

The hard part is no longer matrix enumeration; it is only partition enumeration together with a small rank computation for each pair of cycle types.

Worked Example: $n=3$

The partitions of $3$ are $1+1+1$, $2+1$, and $3$, with permutation counts $1$, $3$, and $2$. Using the orbit formula and the valuation rules above, the fixed-count matrix is

$$\left[\operatorname{Fix}(\lambda,\mu)\right]= \begin{pmatrix} 16 & 4 & 1 \\ 4 & 4 & 1 \\ 1 & 1 & 4 \end{pmatrix}.$$

For instance, for $(2+1,2+1)$ we have $o=2+1+1+1=5$. The four equation nodes split so that three independent constraints remain, hence $\rho=3$ and $\operatorname{Fix}=2^{5-3}=4$.

Burnside now gives

$$c_3=\frac{1}{36}(1,3,2) \begin{pmatrix} 16 & 4 & 1 \\ 4 & 4 & 1 \\ 1 & 1 & 4 \end{pmatrix} \begin{pmatrix} 1\\ 3\\ 2 \end{pmatrix} =\frac{108}{36}=3.$$

This small case already shows the full mechanism: cycle types determine orbit counts, 2-adic valuations determine the parity rank, and Burnside averages the resulting fixed counts.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they enumerate all partitions of $n$ and store, for each cycle type, its cycle-length multiplicities, its list of 2-adic valuations, and the number of permutations with that type. Next they precompute $\gcd(p,q)$ for every $1\le p,q\le n$.

For each pair of cycle types, the implementation computes $o(\lambda,\mu)$ by summing the relevant gcd values. It then builds the small cycle-level equality/forced-zero structure described above, counts the number of unforced connected components, obtains the parity rank $\rho(\lambda,\mu)$, and adds

$$N(\lambda)N(\mu)\,2^{\,o(\lambda,\mu)-\rho(\lambda,\mu)}$$

to an exact integer accumulator. After all type pairs are processed, the total is divided by $(n!)^2$ as required by Burnside's lemma, and only then is the result reduced modulo $1001001011$. The C++ implementation parallelizes the outer traversal, while the Python and Java implementations use the same mathematics in a direct single-process sum.

Complexity Analysis

Let $p(n)$ be the partition number. Enumerating all cycle types costs $O(p(n)\,n)$ space and roughly the same order of bookkeeping work. The main Burnside sum runs over $p(n)^2$ type pairs. For each pair, the orbit count and the parity-rank computation each need at most $O(n^2)$ elementary operations, because there are at most $n$ row cycles and $n$ column cycles.

Thus the overall time complexity is $O(p(n)^2 n^2)$, and the memory usage is $O(p(n)\,n)$ plus the exact-integer accumulator. For $n=20$, this is entirely practical because the number of cycle types is small compared with the astronomical number of matrices themselves.

Footnotes and References

Problem page: https://projecteuler.net/problem=626
Burnside's lemma: Wikipedia - Burnside's lemma
Integer partition: Wikipedia - Partition (number theory)
$p$-adic valuation: Wikipedia - p-adic valuation
Disjoint-set data structure: Wikipedia - Disjoint-set data structure

Problem 626 source code

C++

#include <algorithm>
#include <array>
#include <atomic>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

#include <boost/multiprecision/cpp_int.hpp>

using boost::multiprecision::cpp_int;

namespace {

constexpr std::uint64_t MOD = 1001001011ULL;
constexpr int MAX_N = 20;
constexpr int MAX_VARS = 2 * MAX_N;

struct CycleType {
    std::array<int, MAX_N + 1> mult{};
    std::vector<int> valuations;  // v2 of each cycle length
    int cycles = 0;
    cpp_int perm_count = 0;       // number of permutations with this cycle type
};

int v2(int x) {
    int c = 0;
    while ((x & 1) == 0) {
        x >>= 1;
        ++c;
    }
    return c;
}

std::vector<cpp_int> factorials_up_to(int n) {
    std::vector<cpp_int> fact(static_cast<std::size_t>(n + 1), 1);
    for (int i = 1; i <= n; ++i) fact[static_cast<std::size_t>(i)] = fact[static_cast<std::size_t>(i - 1)] * i;
    return fact;
}

cpp_int ipow_int(int base, int exp) {
    cpp_int result = 1;
    cpp_int b = base;
    int e = exp;
    while (e > 0) {
        if (e & 1) result *= b;
        b *= b;
        e >>= 1;
    }
    return result;
}

void build_cycle_types_rec(int n,
                           int remaining,
                           int min_part,
                           std::vector<int>& current,
                           const std::vector<cpp_int>& fact,
                           std::vector<CycleType>& out) {
    if (remaining == 0) {
        CycleType t;
        t.cycles = static_cast<int>(current.size());
        for (int len : current) {
            ++t.mult[static_cast<std::size_t>(len)];
            t.valuations.push_back(v2(len));
        }

        cpp_int denom = 1;
        for (int len = 1; len <= n; ++len) {
            int m = t.mult[static_cast<std::size_t>(len)];
            if (m == 0) continue;
            denom *= ipow_int(len, m);
            denom *= fact[static_cast<std::size_t>(m)];
        }
        t.perm_count = fact[static_cast<std::size_t>(n)] / denom;
        out.push_back(std::move(t));
        return;
    }

    for (int part = min_part; part <= remaining; ++part) {
        current.push_back(part);
        build_cycle_types_rec(n, remaining - part, part, current, fact, out);
        current.pop_back();
    }
}

std::vector<CycleType> build_cycle_types(int n, const std::vector<cpp_int>& fact) {
    std::vector<CycleType> types;
    std::vector<int> cur;
    build_cycle_types_rec(n, n, 1, cur, fact, types);
    return types;
}

struct DSU {
    std::array<int, MAX_VARS> parent{};
    std::array<unsigned char, MAX_VARS> rank{};
    std::array<unsigned char, MAX_VARS> forced{};

    void init(int n) {
        for (int i = 0; i < n; ++i) {
            parent[static_cast<std::size_t>(i)] = i;
            rank[static_cast<std::size_t>(i)] = 0;
            forced[static_cast<std::size_t>(i)] = 0;
        }
    }

    int find(int x) {
        int r = x;
        while (parent[static_cast<std::size_t>(r)] != r) r = parent[static_cast<std::size_t>(r)];
        while (parent[static_cast<std::size_t>(x)] != x) {
            int p = parent[static_cast<std::size_t>(x)];
            parent[static_cast<std::size_t>(x)] = r;
            x = p;
        }
        return r;
    }

    void unite(int a, int b) {
        a = find(a);
        b = find(b);
        if (a == b) return;
        if (rank[static_cast<std::size_t>(a)] < rank[static_cast<std::size_t>(b)]) std::swap(a, b);
        parent[static_cast<std::size_t>(b)] = a;
        forced[static_cast<std::size_t>(a)] |= forced[static_cast<std::size_t>(b)];
        if (rank[static_cast<std::size_t>(a)] == rank[static_cast<std::size_t>(b)]) ++rank[static_cast<std::size_t>(a)];
    }

    void set_forced(int a) {
        int r = find(a);
        forced[static_cast<std::size_t>(r)] = 1;
    }
};

int parity_rank(const std::vector<int>& rows_v2, const std::vector<int>& cols_v2) {
    const int r = static_cast<int>(rows_v2.size());
    const int s = static_cast<int>(cols_v2.size());
    const int vars = r + s;
    DSU dsu;
    dsu.init(vars);

    for (int i = 0; i < r; ++i) {
        for (int j = 0; j < s; ++j) {
            if (rows_v2[static_cast<std::size_t>(i)] == cols_v2[static_cast<std::size_t>(j)]) {
                dsu.unite(i, r + j);                 // A_i = B_j
            } else if (rows_v2[static_cast<std::size_t>(i)] > cols_v2[static_cast<std::size_t>(j)]) {
                dsu.set_forced(i);                   // A_i = 0
            } else {
                dsu.set_forced(r + j);               // B_j = 0
            }
        }
    }

    std::array<unsigned char, MAX_VARS> seen{};
    int free_components = 0;
    for (int v = 0; v < vars; ++v) {
        int root = dsu.find(v);
        if (seen[static_cast<std::size_t>(root)]) continue;
        seen[static_cast<std::size_t>(root)] = 1;
        if (!dsu.forced[static_cast<std::size_t>(root)]) ++free_components;
    }
    return vars - free_components;
}

int orbit_count(const CycleType& a, const CycleType& b, const int gcd_tbl[MAX_N + 1][MAX_N + 1], int n) {
    int o = 0;
    for (int p = 1; p <= n; ++p) {
        int mp = a.mult[static_cast<std::size_t>(p)];
        if (mp == 0) continue;
        for (int q = 1; q <= n; ++q) {
            int mq = b.mult[static_cast<std::size_t>(q)];
            if (mq == 0) continue;
            o += mp * mq * gcd_tbl[p][q];
        }
    }
    return o;
}

cpp_int compute_c_exact(int n, unsigned threads) {
    std::vector<cpp_int> fact = factorials_up_to(n);
    std::vector<CycleType> types = build_cycle_types(n, fact);

    int gcd_tbl[MAX_N + 1][MAX_N + 1]{};
    for (int i = 1; i <= n; ++i) {
        for (int j = 1; j <= n; ++j) gcd_tbl[i][j] = std::gcd(i, j);
    }

    const std::size_t tcount = types.size();
    if (threads == 0) threads = 1;
    threads = std::min<unsigned>(threads, static_cast<unsigned>(tcount));
    if (threads == 0) threads = 1;

    std::atomic<std::size_t> next_i{0};
    std::vector<cpp_int> partial(static_cast<std::size_t>(threads), 0);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (unsigned tid = 0; tid < threads; ++tid) {
        pool.emplace_back([&, tid]() {
            cpp_int local = 0;
            while (true) {
                std::size_t i = next_i.fetch_add(1, std::memory_order_relaxed);
                if (i >= tcount) break;
                const CycleType& a = types[i];
                for (std::size_t j = 0; j < tcount; ++j) {
                    const CycleType& b = types[j];
                    int o = orbit_count(a, b, gcd_tbl, n);
                    int rk = parity_rank(a.valuations, b.valuations);
                    int exp = o - rk;
                    assert(exp >= 0);
                    local += a.perm_count * b.perm_count * (cpp_int(1) << exp);
                }
            }
            partial[static_cast<std::size_t>(tid)] = std::move(local);
        });
    }

    for (auto& th : pool) th.join();

    cpp_int total = 0;
    for (const auto& p : partial) total += p;

    cpp_int denom = fact[static_cast<std::size_t>(n)] * fact[static_cast<std::size_t>(n)];
    assert(total % denom == 0);
    return total / denom;
}

std::uint64_t to_u64(const cpp_int& x) {
    return x.convert_to<std::uint64_t>();
}

void run_checks(unsigned threads) {
    assert(to_u64(compute_c_exact(1, 1)) == 1ULL);
    assert(to_u64(compute_c_exact(3, 1)) == 3ULL);
    assert(to_u64(compute_c_exact(5, threads)) == 39ULL);
    assert(to_u64(compute_c_exact(8, threads)) == 656108ULL);
}

}  // namespace

int main() {
    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) threads = 1;

    run_checks(threads);

    cpp_int c20 = compute_c_exact(20, threads);
    std::uint64_t answer = to_u64(c20 % MOD);
    std::cout << answer << "\n";
    return 0;
}

Python

import math

MOD = 1001001011

def v2(x):
    c = 0
    while (x & 1) == 0:
        x >>= 1
        c += 1
    return c

def build_cycle_types(n, fact):
    types = []
    
    def build_rec(remaining, min_part, current):
        if remaining == 0:
            mult = [0] * (n + 1)
            valuations = []
            for length in current:
                mult[length] += 1
                valuations.append(v2(length))
                
            denom = 1
            for length in range(1, n + 1):
                m = mult[length]
                if m == 0: continue
                denom *= (length ** m)
                denom *= fact[m]
                
            perm_count = fact[n] // denom
            types.append({
                'mult': mult,
                'valuations': valuations,
                'perm_count': perm_count
            })
            return
            
        for part in range(min_part, remaining + 1):
            current.append(part)
            build_rec(remaining - part, part, current)
            current.pop()
            
    build_rec(n, 1, [])
    return types

class DSU:
    def __init__(self, n):
        self.parent = list(range(n))
        self.rank = [0] * n
        self.forced = [0] * n
        
    def find(self, x):
        r = x
        while self.parent[r] != r:
            r = self.parent[r]
        curr = x
        while self.parent[curr] != curr:
            p = self.parent[curr]
            self.parent[curr] = r
            curr = p
        return r
        
    def unite(self, a, b):
        a = self.find(a)
        b = self.find(b)
        if a == b: return
        if self.rank[a] < self.rank[b]:
            a, b = b, a
        self.parent[b] = a
        self.forced[a] |= self.forced[b]
        if self.rank[a] == self.rank[b]:
            self.rank[a] += 1
            
    def set_forced(self, a):
        r = self.find(a)
        self.forced[r] = 1

def parity_rank(rows_v2, cols_v2):
    r = len(rows_v2)
    s = len(cols_v2)
    vars = r + s
    dsu = DSU(vars)
    
    for i in range(r):
        for j in range(s):
            if rows_v2[i] == cols_v2[j]:
                dsu.unite(i, r + j)
            elif rows_v2[i] > cols_v2[j]:
                dsu.set_forced(i)
            else:
                dsu.set_forced(r + j)
                
    seen = [0] * vars
    free_components = 0
    for v in range(vars):
        root = dsu.find(v)
        if seen[root]: continue
        seen[root] = 1
        if not dsu.forced[root]:
            free_components += 1
            
    return vars - free_components

def orbit_count(a, b, gcd_tbl, n):
    o = 0
    for p in range(1, n + 1):
        mp = a['mult'][p]
        if mp == 0: continue
        for q in range(1, n + 1):
            mq = b['mult'][q]
            if mq == 0: continue
            o += mp * mq * gcd_tbl[p][q]
    return o

def compute_c_exact(n):
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i
        
    types = build_cycle_types(n, fact)
    
    gcd_tbl = [[0] * (n + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        for j in range(1, n + 1):
            gcd_tbl[i][j] = math.gcd(i, j)
            
    total = 0
    tcount = len(types)
    for i in range(tcount):
        a = types[i]
        for j in range(tcount):
            b = types[j]
            o = orbit_count(a, b, gcd_tbl, n)
            rk = parity_rank(a['valuations'], b['valuations'])
            exp = o - rk
            total += a['perm_count'] * b['perm_count'] * (1 << exp)
            
    denom = fact[n] * fact[n]
    return total // denom

def solve():
    c20 = compute_c_exact(20)
    return str(c20 % MOD)

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler626 {
    static final long MOD = 1001001011L;

    static int v2(int x) {
        int c = 0;
        while ((x & 1) == 0) {
            x >>= 1;
            c++;
        }
        return c;
    }

    static class CycleType {
        int[] mult;
        List<Integer> valuations;
        BigInteger permCount;

        CycleType(int n) {
            mult = new int[n + 1];
            valuations = new ArrayList<>();
        }
    }

    static void buildCycleTypesRec(int n, int remaining, int minPart, List<Integer> current, BigInteger[] fact,
            List<CycleType> out) {
        if (remaining == 0) {
            CycleType t = new CycleType(n);
            for (int len : current) {
                t.mult[len]++;
                t.valuations.add(v2(len));
            }

            BigInteger denom = BigInteger.ONE;
            for (int len = 1; len <= n; len++) {
                int m = t.mult[len];
                if (m == 0)
                    continue;
                denom = denom.multiply(BigInteger.valueOf(len).pow(m));
                denom = denom.multiply(fact[m]);
            }
            t.permCount = fact[n].divide(denom);
            out.add(t);
            return;
        }

        for (int part = minPart; part <= remaining; part++) {
            current.add(part);
            buildCycleTypesRec(n, remaining - part, part, current, fact, out);
            current.remove(current.size() - 1);
        }
    }

    static List<CycleType> buildCycleTypes(int n, BigInteger[] fact) {
        List<CycleType> types = new ArrayList<>();
        buildCycleTypesRec(n, n, 1, new ArrayList<>(), fact, types);
        return types;
    }

    static class DSU {
        int[] parent;
        int[] rank;
        int[] forced;

        DSU(int n) {
            parent = new int[n];
            rank = new int[n];
            forced = new int[n];
            for (int i = 0; i < n; i++) {
                parent[i] = i;
            }
        }

        int find(int x) {
            int r = x;
            while (parent[r] != r)
                r = parent[r];
            int curr = x;
            while (parent[curr] != curr) {
                int p = parent[curr];
                parent[curr] = r;
                curr = p;
            }
            return r;
        }

        void unite(int a, int b) {
            a = find(a);
            b = find(b);
            if (a == b)
                return;
            if (rank[a] < rank[b]) {
                int tmp = a;
                a = b;
                b = tmp;
            }
            parent[b] = a;
            forced[a] |= forced[b];
            if (rank[a] == rank[b])
                rank[a]++;
        }

        void setForced(int a) {
            int r = find(a);
            forced[r] = 1;
        }
    }

    static int parityRank(List<Integer> rowsV2, List<Integer> colsV2) {
        int r = rowsV2.size();
        int s = colsV2.size();
        int vars = r + s;
        DSU dsu = new DSU(vars);

        for (int i = 0; i < r; i++) {
            for (int j = 0; j < s; j++) {
                if (rowsV2.get(i).equals(colsV2.get(j))) {
                    dsu.unite(i, r + j);
                } else if (rowsV2.get(i) > colsV2.get(j)) {
                    dsu.setForced(i);
                } else {
                    dsu.setForced(r + j);
                }
            }
        }

        int[] seen = new int[vars];
        int freeComponents = 0;
        for (int v = 0; v < vars; v++) {
            int root = dsu.find(v);
            if (seen[root] == 1)
                continue;
            seen[root] = 1;
            if (dsu.forced[root] == 0)
                freeComponents++;
        }
        return vars - freeComponents;
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static int orbitCount(CycleType a, CycleType b, int[][] gcdTbl, int n) {
        int o = 0;
        for (int p = 1; p <= n; p++) {
            int mp = a.mult[p];
            if (mp == 0)
                continue;
            for (int q = 1; q <= n; q++) {
                int mq = b.mult[q];
                if (mq == 0)
                    continue;
                o += mp * mq * gcdTbl[p][q];
            }
        }
        return o;
    }

    static BigInteger computeCExact(int n) {
        BigInteger[] fact = new BigInteger[n + 1];
        fact[0] = BigInteger.ONE;
        for (int i = 1; i <= n; i++) {
            fact[i] = fact[i - 1].multiply(BigInteger.valueOf(i));
        }

        List<CycleType> types = buildCycleTypes(n, fact);

        int[][] gcdTbl = new int[n + 1][n + 1];
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                gcdTbl[i][j] = gcd(i, j);
            }
        }

        BigInteger total = BigInteger.ZERO;
        int tcount = types.size();

        for (int i = 0; i < tcount; i++) {
            CycleType a = types.get(i);
            for (int j = 0; j < tcount; j++) {
                CycleType b = types.get(j);
                int o = orbitCount(a, b, gcdTbl, n);
                int rk = parityRank(a.valuations, b.valuations);
                int exp = o - rk;
                BigInteger term = a.permCount.multiply(b.permCount).multiply(BigInteger.valueOf(2).pow(exp));
                total = total.add(term);
            }
        }

        BigInteger denom = fact[n].multiply(fact[n]);
        return total.divide(denom);
    }

    public static String solve() {
        BigInteger c20 = computeCExact(20);
        return c20.remainder(BigInteger.valueOf(MOD)).toString();
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}