Enumerating $4 \times 4$ matrices satisfying parity constraints

Question

I've encountered a problem, which requires computer aid, but it seems a little above my Mathematica prowess because it requires counting objects satisfying some simple conditions. It seems doable, and I'm hoping someone might be able to help me out!

The problem is simple. For a fixed $N \in \mathbb{Z}_{\geq 0}$ I'm hoping for a function, say MatCount[N_] which counts how many $4 \times 4$ matrices

$$B_{N}=\begin{bmatrix}b_{00} & b_{01} & b_{02} & b_{03} \\b_{10} & b_{11} & b_{12}& b_{13} \\ b_{20} & b_{21} & b_{22} & b_{23} \\b_{30} & b_{31} & b_{32} & b_{33}\end{bmatrix}$$

there are such that:

1. $b_{ij} \in \mathbb{Z}_{\geq0}$ and $\sum_{i,j} b_{ij} = N$, and

2. If you sum up the elements in the last 3 of the 4 columns, you get either even, even, even or odd, odd, odd. Similarly, if you sum up the elements in the last 3 of the 4 rows, you get either even, even, even or odd, odd, odd.

That last condition is a bit of a mouthful. Basically, I'm interested in the parity of the sums over columns and rows. I don't care about the parity of the first column and first row. But the last 3 columns must be of the same parity and the last 3 rows must be of the same parity (but possibly different than the column's parity).

Can someone help me out with constructing MatCount[N_], or at least an outline of what I should be thinking about? Thanks!

EDIT: If the problem is combinatorially too challenging, I would be content with a solution for "small" $N$. Having even the first 10 or so values would be nice. They quickly become a bear to compute by hand!

Answer 1

You can use SatisfiabilityCount for this. First, represent each matrix element in binary so that:

$$b_{i,j}=\sum _{k=0}^d b(i,j,k) 2^k$$

where the largest number is less than $2^{d+1}$.

Then, the constraint that column $j$ is even can be represented as:

evenColumn[j_] := BooleanCountingFunction[{{0, 2, 4}}, Table[b[i, j, 0], {i, 4}]]

or odd:

oddColumn[j_] := BooleanCountingFunction[{{1, 3}}, Table[b[i, j, 0], {i, 4}]]

Similarly, row $i$ is even:

evenRow[i_] := BooleanCountingFunction[{{0, 2, 4}}, Table[b[i, j, 0], {j, 4}]]

or odd:

oddRow[i_] := BooleanCountingFunction[{{1, 3}}, Table[b[i, j, 0], {j, 4}]]

The constraint that the matrix elements add up to $n$ can also be represented in this way. For example, 5 can be represented with:

1. all the $2^2$ and $2^1$ bits are 0 and five of the $2^0$ bits are 1
2. all the $2^2$ bits are 0 and one of the $2^1$ bits is 1 and three of the $2^0$ bits are 1
3. all the $2^2$ bits are 0 and two of the $2^1$ bits are 1 and one of the $2^0$ bits is 1
4. one of the $2^2$ bits is 1, all of the $2^1$ bits are 0 and one of the $2^0$ bits is 1.

Here is a function that creates this boolean representation (probably could be a bit simpler):

total[n_] := Module[{twos = 2^Range[0, Floor @ Log[2, n]]},
    Or @@ (tobool[Counts[Join[twos, #]]-1]& /@ IntegerPartitions[n, All, twos])
]

tobool[counts_] := And @@ KeyValueMap[
    Function[{k, v}, BooleanCountingFunction[{v}, Tuples[b[Range[4], Range[4], {Log[2,k]}]]]],
    counts
]

A function that computes the counts:

count[n_] := SatisfiabilityCount[
    (oddRow[2] && oddRow[3] && oddRow[4] || evenRow[2] && evenRow[3] && evenRow[4]) &&
    (oddColumn[2] && oddColumn[3] && oddColumn[4] || evenColumn[2] && evenColumn[3] && evenColumn[4]) &&
    total[n]
]

Finally, the results up to $n = 30$:

Print[AbsoluteTiming[# -> count[#]]]& /@ Range[30];

{0.001346,1->1}

{0.00302,2->16}

{0.01568,3->51}

{0.028007,4->276}

{0.072154,5->969}

{0.140554,6->3504}

{0.240101,7->10659}

{0.38415,8->30954}

{0.46881,9->81719}

{0.680237,10->205040}

{0.75767,11->482885}

{1.06645,12->1088100}

{1.24805,13->2340135}

{1.87281,14->4850640}

{2.02903,15->9694845}

{2.61581,16->18789795}

{2.88116,17->35357670}

{3.59111,18->64833120}

{4.05722,19->115997970}

{5.05127,20->203014680}

{5.21446,21->347993910}

{6.39616,22->585292320}

{7.19198,23->966955410}

{8.62268,24->1571349780}

{8.91777,25->2514084066}

{11.093,26->3964589856}

{11.6027,27->6167026726}

{14.3516,28->9470900056}

{14.6145,29->14369476066}

{17.1553,30->21554373984}

The results agree with the available results obtained by @MikeY's brute force method.

Answer 2

Complete update...

OK, a method that solves it directly and is extremely fast for all sizes.

Since we are only concerned with parity relations, the only thing that matters is whether an element of $B_{N}$ is odd or even. Let an odd element be marked by 1 and an even element be marked by 0. Call the matrix $P_{N}$ to denote a parity matrix. There are only $2^{16}=65536$ unique parity matrices. Of those, we can identify in advance the ones that satisfy the parity constraint.

tups = Tuples[{0, 1}, 16];

cmb[vec_] := Module[{mat = Partition[vec, 4], row, col},

  row = Equal @@ (Boole /@ OddQ /@ Total /@ Rest@mat);
  col = Equal @@ (Boole /@ OddQ /@ Total /@ Rest@Transpose@mat);

  row && col
  ]

goodTups = Select[tups, cmb[#] &];

There are 4096 of them. but we don't actually care about the details of those matrices, just the number of odd elements that they have. So tally them up, counting the number with 0 odd elements, 1, etc.

qt = Tally[Total /@ goodTups] // Sort
(* {{0, 1}, {1, 1}, {3, 35}, {4, 140}, {5, 273}, {6, 448}, {7, 715}, 
   {8,870}, {9, 715}, {10, 448}, {11, 273}, {12, 140}, {13, 35}, 
   {15, 1}, {16, 1}} *)

Now we need a function that calculates how many permutations of integers satisfy the constraint that they sum to $N$, and that the number of odd elements = $N_{o}$ (and therefore number of even elements is $N_{e}=N-N_{o}$). Note that all odd numbers can be represented as $2i+1$ and even numbers as $2k$, so write

$$ (2i_{1}+1)+...+(2i_{N_{o}}+1)+2k_{1}+...+2k_{N_{e}}=N $$ $$ 2i_{1}+...+2i_{N_{o}}+2k_{1}+...+2k_{N_{e}}=N-N_{o} $$ $$ i_{1}+...+i_{N_{o}}+k_{1}+...k_{N_{e}}=\frac{N-N_{o}}{2} $$

We can count those using the property that the number of perms of integers with 16 elements summing to $m$ is

$$ \binom{m+15}{15} $$

Note that if $m$ is fractional, let the answer be zero.

numPPO[n_, no_] := If[EvenQ@(n - no), Binomial[(n - no)/2 + 15, 15], 0]

Wrap it up into a function

numParityPerms[n_] := Total@Map[numPPO[n, #[[1]]]*#[[2]] &, qt]

Running it for $N$ from 1 to 15 gives this table. The third column is the ratio of perms satisfying the constraint to total perms. Generating the table is instantaneous.

$$ \left( \begin{array}{ccc} 0 & 1 & 1 \\ 1 & 1 & \frac{1}{16} \\ 2 & 16 & \frac{2}{17} \\ 3 & 51 & \frac{1}{16} \\ 4 & 276 & \frac{23}{323} \\ 5 & 969 & \frac{1}{16} \\ 6 & 3504 & \frac{146}{2261} \\ 7 & 10659 & \frac{1}{16} \\ 8 & 30954 & \frac{469}{7429} \\ 9 & 81719 & \frac{1}{16} \\ 10 & 205040 & \frac{466}{7429} \\ 11 & 482885 & \frac{1}{16} \\ 12 & 1088100 & \frac{465}{7429} \\ 13 & 2340135 & \frac{1}{16} \\ 14 & 4850640 & \frac{13474}{215441} \\ 15 & 9694845 & \frac{1}{16} \\ \end{array} \right) $$

Note the regularity in the third column. All odd values of $N$ result in the ratio being $1/16$. By fiddling around with the even number ratio terms I stumbled onto an orderly relationship, which leads to a function r[ ]

Even terms only for the ratio from the table above:

je = {1, 2/17, 23/323, 146/2261, 469/7429, 466/7429, 465/7429};

Ratios@(je - 1/16);

(* {1/17, 3/19, 5/21, 7/23, 9/25, 11/27}  *)

From this we can derive

r[0] = 1;
r[n_?OddQ] = 1/16;
r[n_?EvenQ] := (r[n - 2] - 1/16) (n - 1)/(n + 15) + 1/16;

parityPermsR[n_] := Binomial[n + 15, 15] r[n];

While not formally proven, it checks out to as large a number as will compute.

parityPermsR[2001] == numParityPerms[2001] // Timing    
(* 0., True *)