Number of equivalence classes of $w \times h$ matrices under switching rows and columns

Question

If I have a $w \times h$ matrix where each value is an integer $0 \lt n \lt 20$,

how can I count the number of distinct configurations, where $2$ configurations are "distinct" if there is no way to reshuffle the rows and columns that would produce the same matrix?

Can this be counted with the stars and bars method?

---

For example, these are equal (we swapped a row, then a column):

0 0 0    2 0 4
0 2 4    0 0 0

but these are distinct (no way to swap rows or columns to produce the other):

0 0 0    2 0 0
0 2 4    0 4 0

It seems like there ought to be a way to count the rows or columns as "bins" and the values as balls. I realize that in this case there are $18$ different colored balls, but even if the only values possible were $1$ and $0$, (ball or no ball) I can't see how to represent it as stars and bars.

Answer 1

This has a very straightforward answer using the Burnside lemma. With $n$ rows, $m$ columns and $q$ possible values we simply compute the cycle index of the cartesian product group ($S_n \times S_m$, consult Harary and Palmer, Graphical Enumeration, section 4.3) and evaluate it at $a[p]=q$ as we have $q$ possibilities for an assignment that is constant on the cycle. The cycle index is easy too -- for two cycles of length $p_1$ and $p_2$ that originate in a permutation $\alpha$ from $S_n$ and $\beta$ from $S_2$ the contribution is $a[\mathrm{lcm}(p_1, p_2)]^{\gcd(p_1, p_2)}.$

We get for a $3\times3$ the following colorings of at most $q$ colors:

$$1, 36, 738, 8240, 57675, 289716, 1144836, 3780288,\ldots$$

which points us to OEIS A058001 where these values are confirmed.

We get for a $4\times 4$ the following colorings of at most $q$ colors:

$$1, 317, 90492, 7880456, 270656150, 4947097821, \\ 58002778967, 490172624992,\ldots$$

which points us to OEIS A058002 where again these values are confirmed.

We get for a $5\times 5$ the following colorings of at most $q$ colors:

$$1, 5624, 64796982, 79846389608, 20834113243925, 1979525296377132, \\ 93242242505023122, 2625154125717590496,\ldots$$

which points us to OEIS A058003 where here too these values are confirmed.

This was 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_cycleind_symmNM :=
proc(n, m)
local indA, indB, res, termA, termB, varA, varB,
    lenA, lenB, instA, instB, p, lcmv;
option remember;

    if n=1 then
        return pet_cycleind_symm(m);
    else
        indA := pet_cycleind_symm(n);
    fi;

    if m=1 then
        return pet_cycleind_symm(n);
    else
        indB := pet_cycleind_symm(m);
    fi;

    res := 0;

    for termA in indA do
        for termB in indB do
            p := 1;

            for varA in indets(termA) do
                lenA := op(1, varA);
                instA := degree(termA, varA);

                for varB in indets(termB) do
                    lenB := op(1, varB);
                    instB := degree(termB, varB);

                    lcmv := lcm(lenA, lenB);
                    p :=
                    p*a[lcmv]^(instA*instB*lenA*lenB/lcmv);
                od;
            od;

            res := res + lcoeff(termA)*lcoeff(termB)*p;
        od;
    od;

    res;
end;



mat_count :=
proc(n, m, q)

    subs([seq(a[p]=q, p=1..n*m)],
         pet_cycleind_symmNM(n, m));
end;

Addendum Nov 17 2018. Note that a product of powers of variables implements the multiset of cycles concept through indets (distinct elements) and degree (number of occurrences). Here we iterate over pairs of monomials representing a conjugacy class from $Z(S_n)$ and $Z(S_m)$ and compute $a[\mathrm{lcm}(p_1, p_2)]^{\gcd(p_1, p_2)}$ for pairs of cycles $a_{p_1}$ and $a_{p_2}.$ This makes for the highly compact algorithm shown above, which will produce e.g. for a three by four,

$${\frac {{a_{{1}}}^{12}}{144}}+1/24\,{a_{{1}}}^{6}{a_{{2}}}^{3} +1/18\,{a_{{1}}}^{3}{a_{{3}}}^{3}+1/12\,{a_{{2}}}^{6} \\+1/6\,{a_{{4}}}^{3}+1/48\,{a_{{1}}}^{4}{a_{{2}}}^{4} +1/8\,{a_{{2}}}^{5}{a_{{1}}}^{2}+1/6\,a_{{1}}a_{{2}}a_{{3}}a_{{6}} \\+1/8\,{a_{{3}}}^{4}+1/12\,{a_{{3}}}^{2}a_{{6}} +1/24\,{a_{{6}}}^{2}+1/12\,a_{{12}}.$$

Answer 2

After struggling with this problem for a couple weeks and attempting to understand the given code and explanation, I believe I have come up with a somewhat more elegant solution for Python. For those like me who have very little experience with combinatorics, I am also including my explanation of the math behind the code that will hopefully be easy to understand for people new to this stuff. First, the solution in Python (interactive example here):

from math import factorial
from fractions import Fraction
import math

def answer(w, h, s):
    total = 0 # initialize return value
    # generate cycle indices for the set of rows and set of columns
    cycidx_cols = cycle_index(w)
    cycidx_rows = cycle_index(h)
    # combine every possible pair of row and column permutations
    for col_coeff, col_cycle in cycidx_cols:
        for row_coeff, row_cycle in cycidx_rows:
            coeff = col_coeff * row_coeff # combine coefficients
            cycle = combine(col_cycle, row_cycle) # combine cycles
            # substitute each variable for s
            value = 1
            for x, power in cycle:
                value *= s ** power
            # multiply by the coefficient and add to the total
            total += coeff * value
    return str(total)

## combines sets of variables with their coefficients to generate a complete cycle index
## returns [ ( Fraction:{coeff}, [ ( int:{length}, int:{frequency} ):{cycle}, ... ]:{cycles} ):{term}, ... ]
def cycle_index(n):
    return [(coeff(term), term) for term in gen_vars(n, n)]

## calculates the coefficient of a term based on values associated with its variable(s)
## this is based off part of the general formula for finding the cycle index of a symmetric group
def coeff(term):
    val = 1
    for x, y in term:
        val *= factorial(y) * x ** y
    return Fraction(1, val)

## generates the solution set to the problem: what are all combinations of numbers <= n that sum to n?
## this corresponds to the set of variables in each term of the cycle index of symmetric group S_n
def gen_vars(n, lim):
    soln_set = [] # store the solution set in a list
    if n > 0: # breaks recursive loop when false and returns an empty list
        for x in range(lim, 0, -1): # work backwards from the limit
            if x == 1: # breaks recursive loop when true and returns a populated list
                soln_set.append([(1, n)])
            else: # otherwise, enter recursion based on how many x go into n
                for y in range(int(n / x), 0, -1):
                    # use recursion on the remainder across all values smaller than x
                    recurse = gen_vars(n - x * y, x - 1)
                    # if recursion comes up empty, add the value by itself to the solution set
                    if len(recurse) == 0:
                        soln_set.append([(x, y)])
                    # otherwise, append the current value to each solution and add that to the solution set
                    for soln in recurse:
                        soln_set.append([(x, y)] + soln)
    return soln_set # return the list of solutions

## combines two terms of a cycle index of the form [ ( int:{length}, int:{frequency} ):{cycle}, ... ]
def combine(term_a, term_b):
    combined = []
    # combine all possible pairs of variables
    for len_a, freq_a in term_a:
        for len_b, freq_b in term_b:
            # new subscript = lcm(len_a, len_b)
            # new superscript = len_a * freq_a * len_b * freq_b / lcm(len_a, len_b)
            lcm = len_a * len_b / math.gcd(len_a, len_b)
            combined.append((lcm, int(len_a * freq_a * len_b * freq_b / lcm)))
    return combined

Now, the explanation: We are asked to find the number of unique matrices given the width $w$, height $h$, and number of possible values $s$. Normally, this would be as simple as counting permutations, which would give us $(w \cdot h)^s$ unique matrices. However, the challenge of this problem comes from the equivalency relationship defined by the ability to shuffle the rows and columns of the matrix. So, we should first consider what happens when we shuffle around rows and columns. We will begin by considering the set of rows separately from the set of columns, but the same methods can be applied to both. Later, we will combine the two results to create a representation of the whole matrix.

We will begin by figuring out the different possible ways to transform one row into another. (In a matrix, this would be equivalent to shuffling the order of the columns.) Let us consider a row of length 4. One possible transformation on would be $\begin{pmatrix}1&2&3&4\\3&1&2&4\end{pmatrix}$, where the top row transforms into the bottom row. If we continually apply this transformation on the same row, you will notice that the value in position 4 stays put while the other three values will follow the cycle $1\rightarrow3\rightarrow2\rightarrow1$. Interestingly, every single possible transformation can be mapped to a unique group of cycles. For example, the above transformation can be mapped to the cycle group $g_8=(132)(4)$. This is one of $4!=24$ unique cycle groups for a row or column of length 4. The complete list is shown here:

$$G=\{(1234), (1243), (1324), (1342), (1423), (1432), (123)(4), (132)(4), (124)(3), (142)(3), (134)(2), (143)(2), (234)(1), (243)(1), (12)(34), (13)(24), (14)(23), (12)(3)(4), (13)(2)(4), (14)(2)(3), (23)(1)(4), (24)(1)(3), (34)(1)(2), (1)(2)(3)(4)\}$$

You may notice that the cycle groups can be categorized into five unique types (represented with five unique terms): $a_4=(abcd)$, $a_1a_3=(abc)(d)$, $a_2^2=(ab)(cd)$, $a_1^2a_2=(ab)(c)(d)$, $a_1^4=(a)(b)(c)(d)$, where each variable $a_p^q$ represents a cycle of length $p$ appearing $q$ times in the cycle group. We can generate the complete list of these types for any $n$ by answering the question, "What are all of the different ways for a set of numbers $\{x \in X : 1 \leq x \leq n\}$ to sum to $n$?" For $n=4$, this would be $(4)$, $(3+1)$, $(2+2)$, $(2+1+1)$, and $(1+1+1+1)$. We can rewrite these as a set of vectors $\textbf{j}=(j_1,j_2,j_3,j_4)$, where $j_x$ represents the frequency of $x$ in the sum:

$$J_4=\{(0,0,0,1),(1,0,1,0),(0,2,0,0),(2,1,0,0),(4,0,0,0)\}$$

We will make use of this set later. The function gen_vars(n, lim) recursively generates $J_n$ for any $n$ (initially, lim == n). However, it is returned by the function in the form of a list of lists of pairs of integers [[(p,q),...],...] where each inner list represents a unique sum and each pair represents the value p and its frequency q in the sum. This list representation speeds up calculations later on.

Returning to the notation $a_p^q$ representing cycles, we can form an equation that represents the entire set of possible cycle groups. We do this by adding each of these terms multiplied by their frequency in $G$:

$$6a_4+8a_1a_3+3a_2^2+6a_1^2a_2+a_1^4$$

Furthermore, if we divide the entire polynomial by the total number of cycles, we get each term's contribution to the complete set of cycle groups:

$$\frac{1}{4}a_4+\frac{1}{3}a_1a_3+\frac{1}{8}a_2^2+\frac{1}{4}a_1^2a_2+\frac{1}{24}a_1^4=Z(S_4)$$

This is known as the cycle index $Z(X)$ for the symmetric group $S_4$. This link includes the cycle indices for the first 5 symmetric groups, and you can reverse these steps to verify that each $Z(S_n)$ accurately represents all possible cycle groups for a set of length $n$. Importantly, we are also given a general formula for finding the cycle index of any $S_n$ (cleaned up a little):

$$Z(S_n)=\sum_{\textbf{j} \in J_n} \left(\frac{1}{\prod_{k=0}^n(k^{j_k} \cdot j_k!)}a_1^{j_1}a_2^{j_2}...a_n^{j_n}\right)$$

This is where that set $J_4$ from earlier comes into play. Indeed, if you plug in the associated values, you will come up with the cycle index for the symmetric group $S_4$. The function coeff(term) calculates the $\frac{1}{\prod_{k=0}^n(k^{j_k} \cdot j_k!)}$ portion of the equation. The cycle_index(n) function puts the coefficients with their terms, returning a list that is representative of the appropriate cycle index.

The cycle index will tell us how many different rows are possible such that no row can be transformed into another using any of the transformations that we found. All we have to do is plug in the number of possible values $s$ in for each variable $a_x$ in our equation (regardless of the value of $x$). For example, if we use $s=3$, we find that there should be 15 unique rows. Here is the list of all possible rows for $s=3$ to confirm this result:

$$R=\{(1,1,1,1),(1,1,1,2),(1,1,1,3),(1,1,2,2),(1,1,2,3),(1,1,3,3),(1,2,2,2),(1,2,2,3),(1,2,3,3),(1,3,3,3),(2,2,2,2),(2,2,2,3),(2,2,3,3),(2,3,3,3),(3,3,3,3)\}$$

This same result can be found using the formula for combinations with replacement, however, this equation fails when applied to a matrix, which is why we are using cycle indices. So, once the cycle indices have been calculated for both the set of rows and the set of columns in our matrix, we must combine them to form the cycle index for the entire matrix. This is done term by term, combining each term of the first with each term in the second. Marko Riedel has an excellent step-by-step explanation of how to do this for a $2 \times 3$ matrix in another post linked here. However, I would like to clarify one part that confused me when I first read it. In order to combine two variables $a_p^q$ and $b_x^y$, use the following template (where $\text{lcm}(a,b)$ is the least common multiple of $a$ and $b$):

$$C(a_p^q,b_x^y)=a_{\text{lcm}(p,x)}^{p\cdot q\cdot x\cdot y/\text{lcm}(p,x)}$$

The combining of terms (ignoring the coefficients, which are multiplied in answer(w, h, s)) is done by the function combine(term_a, term_b) which returns the combined term. This entire process is brought together within the function answer(w, h, s). It calls each of the other functions in turn to create the cycle index for the matrix then plugs in $s$ for each variable to give us our final result.

Hope this helps! I will be more than happy to clarify anything in the comments.

Answer 3

Some have asked me about my Python version. It turns out python is missing a lot of what maple provides for symbolic manipulation. Here is my python version. It follows very closely @Marko Riedel's version, and executes on my machine in 0.6 seconds:

from fractions import *
from copy import *


def expand(frac, terml):
    for term in terml:
        term[0] *= frac
    return terml


def multiplyTerm(sub, terml):
    terml = deepcopy(terml)
    for term in terml:
        alreadyIncluded = False
        for a in term[1]:    # term[1] is a list like [[1,1],[2,3]]  where the
            if a[0] == sub:  # first item is subscript and second the exponent
                alreadyIncluded = True
                a[1] += 1
                break
        if not alreadyIncluded:
            term[1].append([sub, 1])

    return terml


def add(termla, termlb):
    terml = termla + termlb

    # now combine any terms with same a's
    if len(terml) <= 1:
        return terml
    #print "t", terml
    for i in range(len(terml) - 1):
        for j in range(i + 1, len(terml)):
            #print "ij", i, j
            if set([(a[0], a[1]) for a in terml[i][1]]) == set([(b[0], b[1]) for b in terml[j][1]]):
                terml[i][0] = terml[i][0] + terml[j][0]
                terml[j][0] = Fraction(0, 1)

    return [term for term in terml if term[0] != Fraction(0, 1)]


def lcm(a, b):
    return abs(a * b) / gcd(a, b) if a and b else 0

pet_cycnn_cache = {}
def pet_cycleind_symm(n):
    global pet_cycnn_cache
    if n == 0:
        return [ [Fraction(1.0), []] ]

    if n in pet_cycnn_cache:
        #print "hit", n
        return pet_cycnn_cache[n]

    terml = []
    for l in range(1, n + 1):
        terml = add(terml, multiplyTerm(l,  pet_cycleind_symm(n - l)) )

    pet_cycnn_cache[n] = expand(Fraction(1, n), terml)
    return pet_cycnn_cache[n]


def pet_cycles_prodA(cyca, cycb):
    alist = []
    for ca in cyca:
        lena = ca[0]
        insta = ca[1]

        for cb in cycb:
            lenb = cb[0]
            instb = cb[1]

            vlcm = lcm(lena, lenb)
            alist.append([vlcm, (insta * instb * lena * lenb) / vlcm])

    #combine terms (this actually ends up being faster than if you don't)
    if len(alist) <= 1:
        return alist

    for i in range(len(alist) - 1):
        for j in range(i + 1, len(alist)):
            if alist[i][0] == alist[j][0] and alist[i][1] != -1:
                alist[i][1] += alist[j][1]
                alist[j][1] = -1

    return [a for a in alist if a[1] != -1]


def pet_cycleind_symmNM(n, m):
    indA = pet_cycleind_symm(n)
    indB = pet_cycleind_symm(m)
    #print "got ind", len(indA), len(indB), len(indA) * len(indB)
    terml = []

    for flatA in indA:
        for flatB in indB:
            newterml = [
                [flatA[0] * flatB[0], pet_cycles_prodA(flatA[1], flatB[1])]
            ]
            #print "b",len(terml)
            #terml = add(terml, newterml)
            terml.extend(newterml)

    #print "got nm"
    return terml


def substitute(term, v):
    total = 1
    for a in term[1]:
    #need to cast the v and a[1] to int or 
    #they will be silently converted to double in python 3 
    #causing answers to be wrong with larger inputs
    total *= int(v)**int(a[1])
    return (term[0] * total)


def answer(w, h, s):
    terml = pet_cycleind_symmNM(w, h)
    #print terml
    total = 0
    for term in terml:
        total += substitute(term, s)

    return int(total)

print answer(12, 12, 20)

Answer 4

Idea

I interpret the Burnside's lemma as: number of orbits = arithmetic average of the fixed points of acts in G on set X. With this in mind, I used the calculation of PET and get something more pythonic (mostly simpler):

Code

I didn't try that hard to shorten the code so I kept things like Polynomial, __str__ to make it easier to play with. For possible a shorter version, replacing Monomial with list and Polynomial with dict can save 20~30 lines.
This is python2 version but with few dependencies. It's not that hard to update it to python3.

from fractions import gcd
from math import factorial

Polynomial = dict

def solution(w, h, s):
    matrix_cycle_index = calc_matrix_cycle_index_no_factorial(w, h)
    exp_coe_map = {}
    for monomial, coefficient in matrix_cycle_index.items():
        exp = monomial.total_exp()
        exp_coe_map[exp] = exp_coe_map.get(exp, 0) + coefficient
    ret = sum([exp_coe_map[exp] * s ** exp for exp in exp_coe_map])
    ret = ret / factorial(w) / factorial(h)
    return str(ret)

class Monomial(list):
    def __init__(self, vars = list()):
        super(Monomial, self).__init__(vars)

    def mul_var(self, var_ind, exp=1):
        while var_ind >= len(self):
            self.append(0)
        self[var_ind] += exp
        return self

    def mul_mono(self, mono):
        m = Monomial()
        for i in range(len(mono)):
            for j in range(len(self)):
                exp_inc = gcd(i + 1, j + 1)
                k = (i + 1) * (j + 1) // exp_inc
                m.mul_var(k - 1, exp_inc * mono[i] * self[j])
        return m

    def total_exp(self):
        return sum(self)

    def __str__(self):
        return '*'.join(['a%d^%d' % (i + 1, self[i]) for i in range(len(self)) if self[i] > 0])

    def __hash__(self):
        # regard self as an integer of base len(self)
        return sum([self[i] * len(self) ** i for i in range(len(self))])

def calc_matrix_cycle_index_no_factorial(w, h):
    polynomial = Polynomial()
    cycleIndexW = calc_cycle_index_z_sn_no_factorial(w)
    cycleIndexH = calc_cycle_index_z_sn_no_factorial(h)
    for entryW in cycleIndexW.items():
        for entryH in cycleIndexH.items():
            m = entryW[0].mul_mono(entryH[0])
            coe = entryW[1] * entryH[1]
            polynomial[m] = polynomial.get(m, 0) + coe
    return polynomial

cycle_indices = [Polynomial({Monomial(): 1})]

def calc_cycle_index_z_sn_no_factorial(setSize):
    if len(cycle_indices) > setSize:
        return cycle_indices[setSize]
    factorial = 1
    ret = Polynomial()
    for k in range(setSize):
        sub_cycle_index_no_factorial = calc_cycle_index_z_sn_no_factorial(setSize - (k + 1))
        for mono, coe in sub_cycle_index_no_factorial.items():
            mono = Monomial(mono).mul_var(k) # copy(mono)
            ret[mono] = ret.get(mono,0) + factorial * coe
        factorial *= setSize - (k + 1)
    if len(cycle_indices) == setSize:
        cycle_indices.append(ret)
    return ret

Algorithm

Understanding the example

I first verified my algorithm with manual calculation:

• Take (w,h,s)=(2,3,4) as an example
• We know that the number of orbits is n=430, |G| = 322 =12. So n|G| = 12 * 430 = 5160. Let's verify $\Sigma_{g \in G} |S_g| = 5160$.
• Write out all elements in G: G = {e, c, 12, 12c, 23, 23c, 13, 13c, 123, 123c, 132, 132c}. |G| = 12. The c switches the two columns and numbers represents permutations of the lines of the matrix.

g	why	calc	count	cycle index
e	all x in S are fixed points	4^6	4096	a1^6
c	losing 3 degrees of freedoms	4^3	64	a2^3
12	losing 2 degrees of freedoms	4^4	256	a1^2 a2^2
12c	losing 3 degrees of freedoms	4^3	64	a3^2
23	losing 2 degrees of freedoms	4^4	256	a1^2 a2^2
23c	losing 3 degrees of freedoms	4^3	64	a3^2
13	losing 2 degrees of freedoms	4^4	256	a1^2 a2^2
13c	losing 3 degrees of freedoms	4^3	64	a3^2
123	only 2 degrees of freedoms left	4^2	16	a2^3
123c	only 1 degree of freedoms left	4^1	4	a6
132	only 2 degrees of freedoms left	4^2	16	a2^3
132c	only 1 degree of freedoms left	4^1	4	a6

"cycle index" of the table: see Burnside's Lemma applied to grids with interchanging rows and columns

The lost degrees of freedoms for a_m^n is n*(m-1) The degrees of freedoms left is sum of the exponential $6-(3-1)*2=2 = 6/3, 6-(2-1)*3=3, 6-(1-1)*6=6, 6-(6-1)*1=1, k-(k/n-1)*n=n, k- (k - n) = n$

And 4096 + 64 + 256 + 64 + 256 + 64 + 256 + 64 + 16 + 4 + 16 + 4 = 5160.

After this I found the "cycle index" way is much easier than my "counting lost degrees of freedom", although these above gave me a solid understanding of how to use PET & Burnside's Lemma.

Cycle Index

To calculate the "cycle index", I modified the last equation/formula of this wiki, and implemented it with the last function in my code, calc_cycle_index_z_sn_no_factorial. But if you are allowed to use sympy, you can do something like below to get complete exponential Bell polynomial.

def ZNF(n):
return sum([
    bell(n, k, [factorial(m-1) * Symbol("a"+str(m)) for m in range(1,n+1) ] )
    for k in range(n+1)]) # bell from package `sympy`

And this would save the trouble of writing calc_cycle_index_z_sn_no_factorial by yourself.

The purpose was to avoid Fraction() but this also saved the trouble of gen_vars and coeff. In @ Koky Puebla 's code.

Conclusion

This is a combinatorics / enumeration problem. I spent quite some time reviewing my algebra. And now I feel more confident about the Burnside's lemma.

Appendix (TL;DR)

Test code for python

  def test():
  def assertEq(actual, expected):
  assert expected == actual, (
  "Expected " + str(expected) + " but got " + str(actual)
  )

    assertEq(solution(2, 2, 2), "7")
    assertEq(solution(2, 3, 4), "430")

    threeSquare = ["1", "36", "738", "8240", "57675", "289716", "1144836", "3780288"]
    for i in range(len(threeSquare)):
        assertEq(solution(3, 3, i + 1), threeSquare[i])  # https://oeis.org/A058001

    fourSquare = [
        "1",
        "317",
        "90492",
        "7880456",
        "270656150",
        "4947097821",
        "58002778967",
        "490172624992",
    ]
    for i in range(len(fourSquare)):
        assertEq(solution(4, 4, i + 1), fourSquare[i])  # https://oeis.org/A058002

    fiveSquare = ["1", "5624", "64796982", "79846389608", "20834113243925"]
    for i in range(len(fiveSquare)):
        assertEq(solution(5, 5, i + 1), fiveSquare[i])  # from https://oeis.org/A058003

    sixSquare = ["1", "251610", "302752867740"]
    for i in range(len(sixSquare)):
        assertEq(solution(6, 6, i + 1), sixSquare[i])  # https://oeis.org/A058004

    binaryMatrices = ["1", "2", "7", "36", "317", "5624", "251610", "33642660"]
    for i in range(len(binaryMatrices)):
        assertEq(solution(i, i, 2), binaryMatrices[i])  # https://oeis.org/A058005

    assertEq(solution(2, 5, 4), "7882")
    assertEq(solution(3, 5, 4), "1757384")

    print("All tests passed")def test():
    def assertEq(actual, expected):
        assert expected == actual,  (
            "Expected " + str(expected) + " but got " + str(actual)
        )

    assertEq(solution(2, 2, 2), "7")
    assertEq(solution(2, 3, 4), "430")

    threeSquare = ["1", "36", "738", "8240", "57675", "289716", "1144836", "3780288"]
    for i in range(len(threeSquare)):
        assertEq(solution(3, 3, i + 1), threeSquare[i])  # https://oeis.org/A058001

    fourSquare = [
        "1",
        "317",
        "90492",
        "7880456",
        "270656150",
        "4947097821",
        "58002778967",
        "490172624992",
    ]
    for i in range(len(fourSquare)):
        assertEq(solution(4, 4, i + 1), fourSquare[i])  # https://oeis.org/A058002

    fiveSquare = ["1", "5624", "64796982", "79846389608", "20834113243925"]
    for i in range(len(fiveSquare)):
        assertEq(solution(5, 5, i + 1), fiveSquare[i])  # from https://oeis.org/A058003

    sixSquare = ["1", "251610", "302752867740"]
    for i in range(len(sixSquare)):
        assertEq(solution(6, 6, i + 1), sixSquare[i])  # https://oeis.org/A058004

    binaryMatrices = ["1", "2", "7", "36", "317", "5624", "251610", "33642660"]
    for i in range(len(binaryMatrices)):
        assertEq(solution(i, i, 2), binaryMatrices[i])  # https://oeis.org/A058005

    assertEq(solution(2, 5, 4), "7882")
    assertEq(solution(3, 5, 4), "1757384")

    print("All tests passed")

Java code (same algorithm)

I wrote Java8 first, but not being able to use Math.BigInteger for some reason (with google), same as for python2.7 and no sympy. So I had to switch to Python. I also adhere the full commented code in Java8, it passes all tests after changing long in solution to (java.Math.)BigInteger and it uses the same algorithm as Python above.

public class Solution {
    public static void main(String[] args) {
        long startTime = System.nanoTime();

        // assertEq(solution(5, 5, 6), "1979525296377132");
        assertEq(solution(2, 2, 2), "7");
        assertEq(solution(2, 3, 4), "430");

        String[] threeSquare = { "1", "36", "738", "8240", "57675", "289716", "1144836", "3780288" };
        for (int i = 0; i < threeSquare.length; i++) {
            assertEq(solution(3, 3, i + 1), threeSquare[i]);
        } // https://oeis.org/A058001

        String[] fourSquare = { "1", "317", "90492", "7880456", "270656150", "4947097821", "58002778967",
                "490172624992" };
        for (int i = 0; i < fourSquare.length; i++) {
            assertEq(solution(4, 4, i + 1), fourSquare[i]);
        } // https://oeis.org/A058002

        String[] fiveSquare = { "1", "5624", "64796982",
                "79846389608", "20834113243925" };
        // cannot test with long , "1979525296377132", "93242242505023122"
        for (int i = 0; i < fiveSquare.length; i++) {
            assertEq(solution(5, 5, i + 1), fiveSquare[i]);
        } // from https://oeis.org/A058003

        String[] sixSquare = { "1", "251610", "302752867740" }; // , "9178323524804624", "28125393244553141210" };
        for (int i = 0; i < sixSquare.length; i++) {
            assertEq(solution(6, 6, i + 1), sixSquare[i]);
        } // https://oeis.org/A058004

        String[] binaryMatrices = {
                "1", "2", "7", "36", "317", "5624", "251610", "33642660" };
        // , "14685630688", "21467043671008","105735224248507784",};
        for (int i = 0; i < binaryMatrices.length; i++) {
            assertEq(solution(i, i, 2), binaryMatrices[i]);
        } // https://oeis.org/A058005

        // I made these up
        assertEq(solution(2, 5, 4), "7882");
        assertEq(solution(3, 5, 4), "1757384");
        long endTime = System.nanoTime();
        long duration = (endTime - startTime);
        System.out.println("Time: " + duration);
    }

    /**
     * Solution of this problem: The count of orbits of the M_{w,h} under group
     * action of S_w x S_h
     *
     * @param w - width of the grid, 1 <= w <= 12
     * @param h - height of the grid, 1 <= h <= 12
     * @param s - number of states of each celestial body, 2 <= s <= 20
     * @return
     */
    public static String solution(int w, int h, int s) {
        Polynomial matrixCycleIndex = calcMatrixCycleIndexNoFactorial(w, h);
        long ret = 0;
        HashMap<Integer, Integer> expCoeMap = new HashMap<>();
        long minExp = Long.MAX_VALUE;
        for (Map.Entry<Monomial, Integer> entry : matrixCycleIndex.entrySet()) {
            int exp = entry.getKey().totalExp();
            minExp = Math.min(minExp, exp);
            expCoeMap.put(exp, expCoeMap.getOrDefault(exp, 0) + entry.getValue());
        }
        for (Map.Entry<Integer, Integer> entry : expCoeMap.entrySet()) {
            long fixedPointNum = entry.getValue() * longPower(s, entry.getKey());
            ret = ret + fixedPointNum;
        }
        for (int factorialLimit : Arrays.asList(w, h)) {
            for (int i = 2; i <= factorialLimit; i++) {
                ret /= i;
            }
        }
        return String.valueOf(ret);
    }

    public static class Monomial { // or extends ArrayList<Integer>
        ArrayList<Integer> vars;

        public Monomial(List<Integer> vars) {
            this.vars = new ArrayList<>(vars);
        }

        public Monomial(Monomial m) {
            this.vars = new ArrayList<>(m.vars);
        }

        public Monomial mulVar(int varInd, int exp) {
            while (varInd >= vars.size()) {
                vars.add(0);
            }
            vars.set(varInd, vars.get(varInd) + exp);
            return this;
        }

        public Monomial mulVar(int varInd) {
            return mulVar(varInd, 1);
        }

        /**
         * Multiply monomial by a special rule.
         * Variable a_m * a_n = a_{lcm(m,n)} ^ gcd(m,n)
         *
         * Test:
         * assertEq(new Mono(Arrays.asList(0, 1))
         * .mulMono(new Mono(Arrays.asList(1, 1))),
         * new Mono(Arrays.asList(0, 3)));
         * assertEq(new Mono(Arrays.asList(2))
         * .mulMono(new Mono(Arrays.asList(3))),
         * new Mono(Arrays.asList(6)));
         *
         * @param mono
         * @return
         */
        public Monomial mulMono(Monomial mono) {
            Monomial m = new Monomial(new ArrayList<>());
            for (int i = 0; i < mono.vars.size(); i++) {
                for (int j = 0; j < this.vars.size(); j++) {
                    int k = lcmOf(i + 1, j + 1); // to 1-based
                    int expInc = gcdOf(i + 1, j + 1); // to 1-based
                    m.mulVar(k - 1, expInc * mono.vars.get(i) * this.vars.get(j)); // TODO: ref: use 1 based
                }
            }
            return m;
        }

        public int totalExp() {
            return this.vars.stream().mapToInt(Integer::intValue).sum();
        }

        @Override
        public String toString() {
            StringBuilder sb = new StringBuilder();
            for (int i = 0; i < this.vars.size(); i++) {
                if (vars.get(i) > 0) {
                    sb.append("a");
                    sb.append(i + 1); // using 1-based index
                    if (vars.get(i) > 1) {
                        sb.append("^");
                        sb.append(vars.get(i));
                    }
                }
            }
            return sb.toString();
        }

        @Override
        public boolean equals(Object obj) {
            if (obj instanceof Monomial) {
                Monomial m = (Monomial) obj;
                return this.vars.equals(m.vars); // or use hashCode
            }
            return false;
        }

        @Override
        public int hashCode() {
            int id = 0;
            int len = vars.size();
            for (int i = 0; i < len; i++) {
                id += vars.get(i) * Math.pow(len, i);
            }
            return id;
        }
    }

    public static class Polynomial extends HashMap<Monomial, Integer> {
    }

    public static Polynomial calcMatrixCycleIndexNoFactorial(int w, int h) {
        Polynomial polynomial = new Polynomial();
        Polynomial cycleIndexW = calcCycleIndexZSnNoFactorial(w);
        Polynomial cycleIndexH = calcCycleIndexZSnNoFactorial(h);
        for (Map.Entry<Monomial, Integer> entryW : cycleIndexW.entrySet()) {
            for (Map.Entry<Monomial, Integer> entryH : cycleIndexH.entrySet()) {
                Monomial m = entryW.getKey().mulMono(entryH.getKey());
                int coe = entryW.getValue() * entryH.getValue();
                polynomial.put(m, polynomial.getOrDefault(m, 0) + coe);
            }
        }
        return polynomial;
    }

    static ArrayList<Polynomial> cycleIndices = new ArrayList<>();

    /**
     * Recursively calculates the partitions of Cycle Index without Factorial(n)
     * Impl of https://en.wikipedia.org/wiki/Cycle_index#Symmetric_group_Sn
     * The coefficients add up to Bell Number B_n
     * This is not https://oeis.org/A178867
     *
     * Algorithm: (1-based index)
     * ZNF(S_n) = Z(S_n) * (n!)
     * Z(S_n) = 1/n \Sigma_{k=1}^n {a_k * Z(S_{n-k})}
     * ZNF(S_n) = \Sigma_{k=1}^n { (n-1)! * a_k * ZNF(S_{n-k}) / (n-k)! }
     * ZNF(S_n) = \Sigma_{k=1}^n { (n-1)! / (n-k)! * a_k * ZNF(S_{n-k}) }
     *
     * To verify: (on https://live.sympy.org/ use expand, should all be 0)
     *
     * def ZNF(n):
     * return sum([
     * bell(n, k, [factorial(m-1) * Symbol("a"+str(m)) for m in range(1,n+1) ] )
     * for k in range(n+1)])
     *
     * expand(ZNF(3)-2*Symbol("a3")*ZNF(0)-2*Symbol("a2")*ZNF(1)-1*Symbol("a1")*ZNF(2))
     * ZNF(4)-6*Symbol("a4")*ZNF(0)-6*Symbol("a3")*ZNF(1)-3*Symbol("a2")*ZNF(2)-1*Symbol("a1")*ZNF(3)
     * expand(ZNF(5)-24*Symbol("a5")*ZNF(0)-24*Symbol("a4")*ZNF(1)-12*Symbol("a3")*ZNF(2)
     * -4*Symbol("a2")*ZNF(3)-1*Symbol("a1")*ZNF(4))
     *
     * @param setSize - the number of elements in the set
     * @return Poly, all monomials and their coefficients of the
     *         polynomial
     *
     * @example
     *          * ZNF(S_0) = 1
     *          * ZNF(S_1) = a1
     *          * ZNF(S_2) = a1^2 + a2
     *          * ZNF(S_3) = a1^3 + 3a1a2 + 2a3
     *          * ZNF(S_4) = a1^4 + 6a1^2a2 + 8a1a3 + 3a2^2 + 6a4
     *          * ZNF(S_5) = a1^5 + 10a1^3a2 + 20a1^2a3 + 15a1a2^2 + 30a1a4 + 20a2a3
     *          + 24a5
     *
     * @see https://live.sympy.org/
     *      A easier way is hard coding the 12 results of
     *      ```
     *      [
     *      sum([
     *      bell(n, k, [factorial(m-1) * Symbol("a"+str(m)) for m in range(1,n+1)])
     *      for k in range(n+1)]) # these 3 lines are same as ZNF(n)
     *      for n in range(1,13)]
     *      ```
     *      in SymPy
     */
    public static Polynomial calcCycleIndexZSnNoFactorial(Integer setSize) {
        if (cycleIndices.size() > setSize) {
            return cycleIndices.get(setSize);
        }
        if (setSize == 0) {
            Polynomial ret = new Polynomial();
            ret.put(new Monomial(new ArrayList<>()), 1);
            if (cycleIndices.isEmpty()) {
                cycleIndices.add(ret);
            }
            return ret;
        }
        int factorial = 1;
        Polynomial ret = new Polynomial();
        for (int k = 0; k < setSize; k++) {
            int eleLeft = (setSize - (k + 1));
            Polynomial subCycleIndexNoFactorial = calcCycleIndexZSnNoFactorial(eleLeft);
            for (Map.Entry<Monomial, Integer> entry : subCycleIndexNoFactorial.entrySet()) {
                Monomial mono = new Monomial(entry.getKey());
                mono.mulVar(k);
                Integer oldCoe = ret.get(mono);
                if (oldCoe == null) {
                    oldCoe = 0;
                }
                ret.put(mono, oldCoe + factorial * entry.getValue());
            }
            factorial *= (eleLeft);
        }
        if (cycleIndices.size() == setSize) {
            cycleIndices.add(ret);
        }
        return ret;
    }

    public static int gcdOf(int m, int n) {
        if (m == 0 || n == 0) {
            return 0;
        }
        if (m == n) {
            return m;
        }
        if (m > n) {
            return gcdOf(m - n, n);
        }
        return gcdOf(m, n - m);
    }

    public static int lcmOf(int m, int n) {
        if (m == 0 || n == 0) {
            return 0;
        }
        return m * n / gcdOf(m, n);
    }

    public static long longPower(long base, long power) {
        long ret = 1;
        for (long i = 0; i < power; i++) {
            ret *= base;
        }
        return ret;
    }

    public static <T> void assertEq(T actual, T expected) {
        if (expected == actual || expected.equals(actual)) {
            return;
        }
        throw new AssertionError("Expected " + expected.toString() + " but got " + actual.toString());
    }
}

Python

In python I did some refactor to improve the readability but it slows down the python execution (from 0.015s to 0.024s on my machine running all the tests below). I think the implementation can still be improved.