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

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