How many arrays with crossed cells, order of rows/columns irrelevant

Question

I've been struggling with this simple problem for months though as I am a newbie to… well, maths, there's high chance someone more educated than myself may get it right!

Let's consider an array or a table or a Ferrers diagram or whatever it's called,

• $r$ rows by $c$ columns,
• in which each cell can either be blank or have a cross in it
• but each column must contain exactly $a$ crosses.

For instance, $\begin{array}{|c|c|c|c|c|c|} \hline ×&×&×&×& & \\ \hline ×&×& & &×&× \\ \hline & &×& &×& \\ \hline & & &×& &× \\ \hline \end{array}$ is such an array with $c=6,\ r=4,\ a=2$.

Now two arrays are said to be identical ('isomporphic' is probably the right word?) if we can obtain the second one from the first one by changing the order of columns and/or rows of the first one. So, the following array is identical to the one above:

$\begin{array}{|c|c|c|c|c|c|} \hline & & &×& &× \\ \hline ×&×& & &×&× \\ \hline &×&×& & & \\ \hline ×& &×&×&×& \\ \hline \end{array}$

but the order of columns is changed to $164325$ and the order of rows to $3241$.

Now the question is: how many unidentical (non-isomorphic) arrays are there, such that all the above conditions are satisfied and $c,r,a$ are given?

I know the answer for $c=6,r=4,a=2$ is $32$ because I have painstakingly went through every combination but how to do this more effectively?? :) I've been trying to learn enough group theory but it's all about rotating cubes or necklaces and I can't figure out how to translate this simple problem so that its tools can be deployed. Also, I tried to find a recurrence (by slicing away the last row and column) but it seems to be dependent on the actual placement of the crosses so…

I'd be glad for any help…

Answer 1

Preliminary remark. None of the sequences that appear in this problem have OEIS entries with one exception. Therefore I was not able to verify the following results except in a few cases for small values. The reader is invited to contribute this type of verification. An algorithm to independently verify the results below for cases where $k\ge 5$ say would be most welcome.

Main text. What I am about to contribute is enrichment material to facilitate additional exploration of this problem. We will treat the case of a square grid containing $n$ rows and columns with $k$ marks being placed in each column. The most important observation is that what we have here is an instance of Power Group Enumeration, with the group acting on the slots where a selection from the ${n\choose k}$ possible column configurations are placed being the symmetric group $S_n$ on $n$ elements and the group $Q_{n,k}$ acting on the column vectors being the action induced on the set of columns by permuting rows.

We can compute the number of configurations by Burnside's lemma which says to average the number of assignments fixed by the elements of the power group, which has $n!\times n!$ elements. But this number is easy to compute. Suppose we have a permutation $\alpha$ from $S_n$ (column permutations) and a permutation $\beta$ from $S_n$ (induced action on the columns by row permutations). If we place the appropriate number of complete, directed and consecutive copies of a cycle from $\beta$ on a cycle from $\alpha$ then this assignment is fixed under the power group action, and this is possible iff the length of $\beta$ divides the length of $\alpha$ and there are as many assignments as the length of $\beta.$ We can work with the cycle indices of $S_n$ and $Q_{n,k}$ and do not need to iterate over all $n!^2$ permutations.

To do this we need the cycle index $Z(Q_{n,k})$ of $Q_{n,k}$, which we compute. as follows: for each permutation shape that occurs in the cycle index $Z(S_n)$ of the symmetric group we compute a representative permutation and apply it to the set of column configurations. The result is factored into disjoint cycles and added to $Z(Q_{n,k})$ with the coefficient it had in $Z(S_n).$ Very simple.

As an example, here is the cycle index for $Q_{6, 2}:$ $${\frac {{a_{{1}}}^{15}}{720}}+1/48\,{a_{{1}}}^{7}{a_{{2}}}^{4}+1 /18\,{a_{{1}}}^{3}{a_{{3}}}^{4}+1/12\,{a_{{1}}}^{3}{a_{{2}}}^{6} +1/4\,a_{{1}}{a_{{4}}}^{3}a_{{2}}\\+1/6\,a_{{2}}{a_{{3}}}^{2}a_{{1 }}a_{{6}}+1/5\,{a_{{5}}}^{3}+1/18\,{a_{{3}}}^{5}+1/6\,{a_{{6}}}^ {2}a_{{3}}$$ and this is $Z(Q_{7, 4})$ $$1/12\,{a_{{1}}}^{3}{a_{{6}}}^{3}{a_{{3}}}^{4}a_{{2}}+1/6\,{a_{{4 }}}^{7}{a_{{2}}}^{3}a_{{1}}+1/7\,{a_{{7}}}^{5}+{\frac {{a_{{1}}} ^{35}}{5040}}+1/6\,{a_{{6}}}^{5}a_{{2}}a_{{3}}\\+{\frac {{a_{{1}}} ^{15}{a_{{2}}}^{10}}{240}}+{\frac {{a_{{1}}}^{5}{a_{{3}}}^{10}}{ 72}}+1/48\,{a_{{2}}}^{14}{a_{{1}}}^{7}+1/10\,{a_{{5}}}^{7}\\+1/48 \,{a_{{2}}}^{16}{a_{{1}}}^{3}+1/18\,{a_{{1}}}^{2}{a_{{3}}}^{11}+ 1/24\,a_{{1}}{a_{{6}}}^{4}{a_{{3}}}^{2}{a_{{2}}}^{2}\\+1/10\,{a_{{ 5}}}^{3}{a_{{10}}}^{2}+1/12\,a_{{4}}{a_{{12}}}^{2}a_{{6}}a_{{1}}$$ and finally this is $Z(Q_{7,5})$ $$1/12\,{a_{{3}}}^{3}a_{{6}}{a_{{1}}}^{2}{a_{{2}}}^{2}+1/8\,{a_{{4 }}}^{4}{a_{{2}}}^{2}a_{{1}}+1/7\,{a_{{7}}}^{3}+{\frac {{a_{{1}}} ^{21}}{5040}}+1/6\,{a_{{6}}}^{3}a_{{3}}\\+{\frac {{a_{{1}}}^{11}{a _{{2}}}^{5}}{240}}+{\frac {{a_{{3}}}^{5}{a_{{1}}}^{6}}{72}}+1/48 \,{a_{{1}}}^{5}{a_{{2}}}^{8}+1/24\,{a_{{4}}}^{4}a_{{2}}{a_{{1}}} ^{3}+1/10\,{a_{{5}}}^{4}a_{{1}}\\+1/48\,{a_{{1}}}^{3}{a_{{2}}}^{9} +1/18\,{a_{{3}}}^{7}+1/24\,a_{{3}}{a_{{6}}}^{2}{a_{{1}}}^{2}{a_{ {2}}}^{2}\\+1/10\,{a_{{5}}}^{2}a_{{10}}a_{{1}}+1/12\,a_{{4}}a_{{2} }a_{{12}}a_{{3}}.$$

The similarity in the coefficients is because these are inherited from the symmetric group.

Now the Burnside computation is best done with a CAS, here is the Maple code.

with(combinat);

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_flatten_term :=
proc(varp)
        local terml, d, cf, v;

        terml := [];

        cf := varp;
        for v in indets(varp) do
            d := degree(varp, v);
            terml := [op(terml), seq(v, k=1..d)];
            cf := cf/v^d;
        od;

        [cf, terml];
end;


pet_autom2cycles :=
proc(src, aut)
        local numa, numsubs;
        local marks, pos, cycs, cpos, clen;

        numsubs := [seq(src[k]=k, k=1..nops(src))];
        numa := subs(numsubs, aut);

        marks := Array([seq(true, pos=1..nops(aut))]);

        cycs := []; pos := 1;

        while pos <= nops(aut) do
              if marks[pos] then
                 clen := 0; cpos := pos;

                 while marks[cpos] do
                       marks[cpos] := false;
                       cpos := numa[cpos];
                       clen := clen+1;
                 od;

                 cycs := [op(cycs), clen];
              fi;

              pos := pos+1;
        od;

        return mul(a[cycs[k]], k=1..nops(cycs));
end;


pet_flat2rep :=
proc(f)
    local p, q, res, cyc, t, len;

    q := 1; res := [];

    for t in f do
        len := op(1, t);
        cyc := [seq(p, p=q+1..q+len-1), q];
        res := [op(res), seq(cyc[p], p=1..nops(cyc))];
        q := q+len;
    od;

    res;
end;

pet_nchoosek_cind :=
proc(n, k)
    option remember;
    local idx_slots, cind, src, aut, q, rep, flat, term;

    cind := 0;

    src := choose(n, k);

    if n=1 then
       idx_slots := [a[1]]
    else
       idx_slots := pet_cycleind_symm(n);
    fi;

    for term in idx_slots do
        flat := pet_flatten_term(term);
        rep := pet_flat2rep(flat[2]);

        aut :=
        map(sel -> sort([seq(rep[sel[q]], q=1..k)]), src);

        cind := cind + flat[1]*pet_autom2cycles(src, aut);
    od;

    cind;
end;


matrix_marks :=
proc(n, k)
    option remember;
    local idx_cols, idx_marks, res, a, b,
    flat_a, flat_b, cyc_a, cyc_b, len_a, len_b, p, q;

    if n=1 then
       idx_cols := [a[1]]
    else
       idx_cols := pet_cycleind_symm(n);
    fi;

    idx_marks := pet_nchoosek_cind(n, k);
    if not type(idx_marks, `+`) then
        idx_marks := [idx_marks];
    fi;

    res := 0;

    for a in idx_cols do
        flat_a := pet_flatten_term(a);

        for b in idx_marks do
            flat_b := pet_flatten_term(b);

            p := 1;
            for cyc_a in flat_a[2] do
                len_a := op(1, cyc_a);
                q := 0;

                for cyc_b in flat_b[2] do
                    len_b := op(1, cyc_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b;
                    fi;
                od;

                p := p*q;
            od;

            res := res + p*flat_a[1]*flat_b[1];
        od;
    od;

    res;
end;

This will produce the following sequence of values for matrices with a single mark per column: $$1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, \ldots$$ which is OEIS A000041.

For two marks per column we get starting at $n=2$ $$1, 3, 11, 35, 132, 471, 1806, 7042, 28494, 118662, 510517, 2262738, \\ 10337474, 48625631,\ldots$$

For three marks per column we get starting at $n=3$ $$1, 5, 35, 410, 6178, 122038, 2921607, 81609320, 2588949454, \\91699869557, 3582942335285, 153048366545566, \\7096576775166579, 355120233277118103,\ldots$$

Finally for four marks per column we get starting at $n=4,$ $$1, 7, 132, 6178, 594203, 85820809, 16341829155, 3875736708590, \\ 1112175913348040, 378860991866916370, 151006214911844288232, \\ 69600017255860985666964, 36729204987785981237238642,\ldots$$ The shared values in these three lists represent the fact that $${n\choose k} = {n\choose n-k}.$$

There is another example of Power Group Enumeration at this MSE link.

Addendum 2014-09-23. By way of an incentive to investigate and develop algorithms to verify the above results for non-trivial values of $n$ and $k$ I present a Perl script (admittedly very simple: warning, exponential algorithm) that can be used to verify small values. It gave the following results.

$ time ./mmpg.pl 5 3
cases: 100000
35

real    0m5.377s
user    0m5.179s
sys     0m0.030s

$ time ./mmpg.pl 6 2
cases: 11390625
132

real    12m24.638s
user    11m35.514s
sys     0m0.390s

$ time ./mmpg.pl 6 3
cases: 64000000
410

real    76m39.028s
user    71m55.408s
sys     0m2.261s

$ time ./mmpg.pl 6 4
cases: 11390625
132

real    16m23.239s
user    15m31.747s
sys     0m0.389s

$ time ./mmpg.pl 7 1
cases: 823543
15

real    2m19.743s
user    2m16.937s
sys     0m0.077s

$ time ./mmpg.pl 8 1
cases: 16777216
22

real    99m47.460s
user    93m49.700s
sys     0m2.667s

This is the Perl script.

#! /usr/bin/perl -w
#

sub permute {
    my ($n) = @_;

    return [[1]] if $n == 1;

    my ($res, $perm) = ([]);

    foreach my $perm (@{ permute($n-1) }){
        my $nxt;

        for(my $pos = 0; $pos < $n; $pos++){
            $nxt = 
                [@$perm[0..($pos-1)], $n, 
                 @$perm[$pos..($n-2)]];
            push @$res, $nxt;
        }
    }

    return $res;
}

sub choose {
    my ($data, $pos, $n) = @_;

    my $size = scalar(@$data);

    return [] if $pos == $size;

    if($n == 1){
        my @res = map { [$_] } @$data[$pos..$size-1];

        return \@res;
    }

    my $rec = choose($data, $pos+1, $n);
    foreach my $sel (@{ choose($data, $pos+1, $n-1) }){
        my $nxt = [@$sel];
        unshift @$nxt, $data->[$pos];

        push @$rec, $nxt;
    }

    return $rec;
}


MAIN: {
    my $n = shift || 4;
    my $k = shift || 3;

    my $range = [0..($n-1)];
    my $cols = choose($range, 0, $k);

    my $opts = scalar(@$cols); my $cases = $opts ** $n;
    print STDERR "cases: $cases\n";

    my $col2ind = {};
    for(my $colind = 0; $colind < $opts; $colind++){
        $col2ind->{join('-', @{ $cols->[$colind] })} 
        = $colind;
    }

    my $seen = {}; my $nperms = permute($n);

    for(my $matind = 0; $matind < $cases; $matind++){
        my $matvec = []; my ($pos, $ind);

        for(($pos, $ind)= (0, $matind); 
            $pos < $n; $pos++){
            my $d = $ind % $opts;

            push @$matvec, $d;
            $ind = ($ind-$d)/$opts;
        }

        for($pos = 0; $pos < $n-1; $pos++){
            last if $matvec->[$pos] > $matvec->[$pos+1];
        }
        next if $pos < $n-1;

        my $admit = 1; my $pid;
        foreach my $perm (@{ $nperms }){
            my @permcols = ();

            foreach(my $colind = 0; $colind < $n; 
                    $colind++){
                my @permcol = 
                    map {
                        $perm->[$_] - 1
                } @{$cols->[$matvec->[$colind]]};


                @permcol = sort { $a <=> $b } @permcol;

                my $pcolind = $col2ind->{join('-', @permcol)};
                push @permcols, $pcolind;
            }

            @permcols = sort { $a <=> $b } @permcols;
            $pid = join('-', @permcols);

            if(exists($seen->{$pid})){
                $admit = undef;
                last;
            }
        }

        $seen->{$pid} = 1 if defined($admit);
    }

    print scalar(keys(%$seen)) . "\n";

    1;
}

The above sequences now have OEIS entries as of today: OEIS A247417, OEIS A247596, OEIS A247597, OEIS A247598.

Addendum 2019-05-01. Consulting the OEIS links from above five years after I first solved this problem by PGE we see that a better, more efficient and more sophisticated algorithm has appeared, placing the above in the category of a historical artefact. Note however that the Maple code admits some improvements which are shown below (duplicate routines have been omitted).

pet_prod2rep :=
proc(varp)
    local v, d, q, res, len, cyc;

    q := 1; res := [];

    for v in indets(varp) do
        d := degree(varp, v);
        len := op(1, v);

        for cyc to d do
            res :=
            [op(res),
             seq(p, p=q+1..q+len-1), q];
            q := q+len;
        od;
    od;

    res;
end;


pet_nchoosek_cind :=
proc(n, k)
    option remember;
    local idx_slots, cind, src, aut, q, rep, term;

    cind := 0;

    src := choose(n, k);

    if n=1 then
       idx_slots := [a[1]]
    else
       idx_slots := pet_cycleind_symm(n);
    fi;

    for term in idx_slots do
        rep := pet_prod2rep(term);

        aut :=
        map(sel -> sort([seq(rep[sel[q]], q=1..k)]), src);

        cind := cind + lcoeff(term)*pet_autom2cycles(src, aut);
    od;

    cind;
end;


matrix_marks :=
proc(n, k)
option remember;
local idx_cols, idx_marks, res, term_a, term_b,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;

    if n=1 then
        idx_cols := [a[1]]
    else
        idx_cols := pet_cycleind_symm(n);
    fi;

    idx_marks := pet_nchoosek_cind(n, k);
    if not type(idx_marks, `+`) then
        idx_marks := [idx_marks];
    fi;

    res := 0;

    for term_a in idx_cols do
        for term_b in idx_marks do
            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in indets(term_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, v_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b*inst_b;
                    fi;
                od;

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;