Enumerating $4 \times 4$ matrices satisfying parity constraints

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.

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 *)