How many $n\times m$ binary matrices are there, up to row and column permutations?
This is solved here using Pólya enumeration theory. For the square case ($n=m$), see this sequence.
Comment: I found these by searching for $1,2,7$ in the OEIS.
Here is a computational contribution that treats the case of a square matrix. As pointed out this problem can be solved using the Polya Enumeration Theorem. In fact if we are interested only in counting these matrices, then the Burnside lemma will suffice. We just need to compute the cycle index of the group acting on the slots of the matrix.
These cycle indices are easy to compute and we do not need to iterate over all $(n!)^2$ pairs of permutations (acting on rows and columns) but instead it is sufficient to iterate over pairs of terms from the cycle index $Z(S_n)$ of the symmetric group $S_n$ according to their multiplicities to obtain the cycle index $Z(Q_n)$ of the combined action on rows and columns. The number of terms here is the much better count of the number of partitions of $n$ squared (upper bound).
Now for a pair of cycles, one of length $l_1$ from a row permutation $\alpha$ and another of length $l_2$ from a column permutation $\beta$ their contribution to the disjoint cycle decomposition product for $(\alpha,\beta)$ in the cycle index $Z(Q_n)$ is by inspection $$a_{\mathrm{lcm}(l_1, l_2)}^{l_1 l_2 / \mathrm{lcm}(l_1, l_2)} = a_{\mathrm{lcm}(l_1, l_2)}^{\gcd(l_1, l_2)}.$$
The algorithm now becomes very simple -- iterate over pairs of terms as described above, collect the contribution from each pair of cycles and add it to the cycle index being computed.
This gives the following cycle indices (only the first four are shown):
$$Z(Q_2) = 1/4\,{a_{{1}}}^{4}+3/4\,{a_{{2}}}^{2},$$
$$Z(Q_3) = 1/36\,{a_{{1}}}^{9}+1/6\,{a_{{1}}}^{3}{a_{{2}}}^{3}+1/4\,a_{{ 1}}{a_{{2}}}^{4}+2/9\,{a_{{3}}}^{3}+1/3\,a_{{3}}a_{{6}},$$
$$Z(Q_4) = {\frac {{a_{{1}}}^{16}}{576}}+1/48\,{a_{{1}}}^{8}{a_{{2}}}^{4 }+1/16\,{a_{{1}}}^{4}{a_{{2}}}^{6}+1/36\,{a_{{1}}}^{4}{a_{{3} }}^{4}+{\frac {17\,{a_{{2}}}^{8}}{192}}\\+1/6\,{a_{{1}}}^{2}a_{ {2}}{a_{{3}}}^{2}a_{{6}}+1/9\,a_{{1}}{a_{{3}}}^{5}+1/12\,{a_{ {2}}}^{2}{a_{{6}}}^{2}+{\frac {13\,{a_{{4}}}^{4}}{48}}+1/6\,a _{{4}}a_{{12}},$$
and
$$Z(Q_5) = {\frac {{a_{{1}}}^{25}}{14400}}+{\frac {{a_{{1}}}^{15}{a_{{2} }}^{5}}{720}}+{\frac {{a_{{1}}}^{9}{a_{{2}}}^{8}}{144}}+{ \frac {{a_{{1}}}^{10}{a_{{3}}}^{5}}{360}}+{\frac {{a_{{1}}}^{ 5}{a_{{2}}}^{10}}{480}}\\+1/48\,{a_{{1}}}^{3}{a_{{2}}}^{11}+{ \frac {a_{{1}}{a_{{2}}}^{12}}{64}}+1/36\,{a_{{1}}}^{6}{a_{{2} }}^{2}{a_{{3}}}^{3}a_{{6}}+1/36\,{a_{{1}}}^{4}{a_{{3}}}^{7}+{ \frac {{a_{{1}}}^{5}{a_{{4}}}^{5}}{240}}\\+{\frac {{a_{{2}}}^{5 }{a_{{3}}}^{5}}{360}}+1/24\,{a_{{1}}}^{3}a_{{2}}{a_{{4}}}^{5} +1/24\,{a_{{1}}}^{2}{a_{{2}}}^{4}a_{{3}}{a_{{6}}}^{2}+1/36\,{ a_{{2}}}^{5}{a_{{3}}}^{3}a_{{6}}\\+1/16\,a_{{1}}{a_{{2}}}^{2}{a _{{4}}}^{5}+1/24\,{a_{{2}}}^{5}a_{{3}}{a_{{6}}}^{2}+1/18\,{a_ {{2}}}^{2}{a_{{3}}}^{5}a_{{6}}+1/16\,a_{{1}}{a_{{4}}}^{6}\\+1/ 36\,{a_{{2}}}^{2}{a_{{3}}}^{3}{a_{{6}}}^{2}+1/12\,{a_{{1}}}^{ 2}a_{{3}}{a_{{4}}}^{2}a_{{12}}+1/12\,a_{{2}}a_{{3}}{a_{{4}}}^ {2}a_{{12}}+{\frac {13\,{a_{{5}}}^{5}}{300}}\\+1/30\,{a_{{5}}}^ {3}a_{{10}}+1/15\,{a_{{5}}}^{2}a_{{15}}+1/20\,a_{{5}}{a_{{10} }}^{2}+1/10\,a_{{5}}a_{{20}}+1/15\,a_{{10}}a_{{15}}.$$
Evaluating these cycle indices with the variables set to two we quickly obtain the sequence
$$2, 7, 36, 317, 5624, 251610, 33642660, 14685630688,\\ 21467043671008, 105735224248507784,1764356230257807614296,\\ 100455994644460412263071692,19674097197480928600253198363072,\\ 13363679231028322645152300040033513414,\\ 31735555932041230032311939400670284689732948,\ldots$$ which is indeed OEIS A002724.
Note that the cycle indices make it possible to enumerate configurations with more than two possible entries or with entries having different weights. For example, with a 3x3 square and three colors $A,B$ and $C$ we get the generating function
$$Z(Q_3)(A+B+C) = 1/36\, \left( A+B+C \right) ^{9}+1/6\, \left( A+B+C \right) ^{3} \left( {A}^{2}+{B}^{2}+{C}^{2} \right) ^{3}\\+2/9\, \left( {A}^{3}+{B}^{3}+{C}^{3} \right)^{3} +1/4\, \left( A+B+C \right) \left( {A}^{2}+{B}^{2}+{C}^{2 } \right) ^{4}\\+1/3\, \left( {A}^{3}+{B}^{3}+{C}^{3} \right) \left( {A}^{6}+{B}^{6}+{C}^{6} \right)$$
which expands to
$${A}^{9}+{A}^{8}B+{A}^{8}C+3\,{A}^{7}{B}^{2}+3\,{A}^{7}B C+3\,{A}^{7}{C}^{2}+6\,{A}^{6}{B}^{3}+10\,{A}^{6}{B}^{2 }C\\+10\,{A}^{6}B{C}^{2}+6\,{A}^{6}{C}^{3}+7\,{A}^{5}{B}^ {4}+17\,{A}^{5}{B}^{3}C+28\,{A}^{5}{B}^{2}{C}^{2}\\+17\,{ A}^{5}B{C}^{3}+7\,{A}^{5}{C}^{4}+7\,{A}^{4}{B}^{5}+22\, {A}^{4}{B}^{4}C+43\,{A}^{4}{B}^{3}{C}^{2}+43\,{A}^{4}{B }^{2}{C}^{3}\\+22\,{A}^{4}B{C}^{4}+7\,{A}^{4}{C}^{5}+6\,{ A}^{3}{B}^{6}+17\,{A}^{3}{B}^{5}C+43\,{A}^{3}{B}^{4}{C} ^{2}+54\,{A}^{3}{B}^{3}{C}^{3}\\+43\,{A}^{3}{B}^{2}{C}^{4 }+17\,{A}^{3}B{C}^{5}+6\,{A}^{3}{C}^{6}+3\,{A}^{2}{B}^{ 7}+10\,{A}^{2}{B}^{6}C+28\,{A}^{2}{B}^{5}{C}^{2}\\+43\,{A }^{2}{B}^{4}{C}^{3}+43\,{A}^{2}{B}^{3}{C}^{4}+28\,{A}^{ 2}{B}^{2}{C}^{5}+10\,{A}^{2}B{C}^{6}+3\,{A}^{2}{C}^{7}\\+ A{B}^{8}+3\,A{B}^{7}C+10\,A{B}^{6}{C}^{2}+17\,A{B}^{5}{ C}^{3}+22\,A{B}^{4}{C}^{4}+17\,A{B}^{3}{C}^{5}\\+10\,A{B} ^{2}{C}^{6}+3\,AB{C}^{7}+A{C}^{8}+{B}^{9}+{B}^{8}C+3\,{ B}^{7}{C}^{2}+6\,{B}^{6}{C}^{3}+7\,{B}^{5}{C}^{4}\\+7\,{B }^{4}{C}^{5}+6\,{B}^{3}{C}^{6}+3\,{B}^{2}{C}^{7}+B{C}^{ 8}+{C}^{9}.$$
This is the Maple code for this computation. Here we have two slightly different ways of evaluating the count, the first by substituting into the cycle index and the second by skipping the cycle index altogether and evaluating all variables at two during processing. The latter should be used when we are interested ony in the count as opposed to classifying configurations according to the number of each color / value that are present.
pet_cycleind_symm :=
proc(n)
option remember;
if n=0 then return 1; fi;
expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;
pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;
res := ind;
polyvars := indets(poly);
indvars := indets(ind);
for v in indvars do
pot := op(1, v);
subs1 :=
[seq(polyvars[k]=polyvars[k]^pot,
k=1..nops(polyvars))];
subs2 := [v=subs(subs1, poly)];
res := subs(subs2, res);
od;
res;
end;
pet_cycleind_sqmat :=
proc(n)
option remember;
local sind, cind, term_a, term_b, v_a, v_b,
len_a, len_b, inst_a, inst_b, p;
cind := 0;
if n=1 then
sind := [a[1]];
else
sind := pet_cycleind_symm(n);
fi;
for term_a in sind do
for term_b in sind do
p := 1;
for v_a in indets(term_a) do
len_a := op(1, v_a);
inst_a := degree(term_a, v_a);
for v_b in indets(term_b) do
len_b := op(1, v_b);
inst_b := degree(term_b, v_b);
p := p*a[lcm(len_a, len_b)]
^(gcd(len_a, len_b)*inst_a*inst_b);
od;
od;
cind := cind +
lcoeff(term_a)*lcoeff(term_b)*p;
od;
od;
cind;
end;
v :=
proc(n)
option remember;
local cind, vars, sbl;
cind := pet_cycleind_sqmat(n);
vars := indets(cind);
sbl := [seq(v=2, v in vars)];
subs(sbl, cind);
end;
w :=
proc(n)
option remember;
local sind, count, term_a, term_b, v_a, v_b,
len_a, len_b, inst_a, inst_b, p;
count := 0;
if n=1 then
sind := [a[1]];
else
sind := pet_cycleind_symm(n);
fi;
for term_a in sind do
for term_b in sind do
p := 1;
for v_a in indets(term_a) do
len_a := op(1, v_a);
inst_a := degree(term_a, v_a);
for v_b in indets(term_b) do
len_b := op(1, v_b);
inst_b := degree(term_b, v_b);
p := p*
2^(gcd(len_a, len_b)*inst_a*inst_b);
od;
od;
count := count +
lcoeff(term_a)*lcoeff(term_b)*p;
od;
od;
count;
end;
This MSE Meta Link has many more PET computations by various users.