Counting Total Number of Non-Equivalent Configurations in a 2-D Grid This is a challenging question I've been trying (unsuccessfully) to solve via programming, math or both. 

Suppose you're given a 2D grid, whose width and height, $w$ and $h$, can each range from $1$ thru $12$ inclusive. Each cell in this grid can be in any of $k$ states (where $k$ ranges from $2$ thru $20$).

The problem is to find the total number of distinct (non-equivalent) state configurations of this grid. 

Definition: Two configurations of a grid, $c_1$ and $c_2$, are deemed "equivalent" if it is possible to change $c_1$ to $c_2$ by swapping any pair of rows and/or any pair of columns (you can perform the swap operation as many times as you want).

To take an example, let's consider a grid with $w=2$, $h=2$, and $k=2$. We will represent each of the two states as $0$ and $1$ respectively.

One possible configuration of the grid is one where all the cells are in state $0$:

00
00

You could swap any rows and/or columns without changing anything. So the above counts as 1 configuration.

Another configuration is where only 1 cell is in state $1$. There are 4 possible arrangements:

01 10 00 00
00 00 01 10

But note that all of the 4 arrangements are "equivalent", because you can get any of the arrangements by swapping rows and/or columns. So the above counts as 1 configuration.

Next, exactly 2 cells are in state $1$:

11 00
00 11

Another (non-equivalent) way in which 2 cells are in state $1$:

10 01
10 01

And yet another:

10 01
01 10

So there are 3 non-equivalent configurations in which exactly 2 cells are in state $1$.

Next, exactly 3 cells are in state $1$:

10 01 11 11
11 11 10 01

1 configuration, just as in the case where only 1 cell was in state $1$.

Finally, all 4 cells are in state $1$:

11
11

1 configuration, just as in the case where all 4 cells were in state $0$.

So, counting them up, there are exactly $7$ distinct configurations for the case where $w=2, h=2, k=2$.

For reference, here are some brute-force test cases I ran on some additional examples (where the w,h,k values are fairly small):

w=2,h=1,k=1 Answer: 1
w=2,h=1,k=2 Answer: 3
w=2,h=1,k=3 Answer: 6
w=2,h=1,k=4 Answer: 10
w=2,h=1,k=5 Answer: 15
w=2,h=1,k=6 Answer: 21
...
So at least when w, h are fixed at 1,2 (or at 2,1), the pattern seems pretty obvious as k increases. Here is another suite of test cases for w=2, h=2:

w=2,h=2,k=1 Answer: 1
w=2,h=2,k=2 Answer: 7
w=2,h=2,k=3 Answer: 27
w=2,h=2,k=4 Answer: 76
w=2,h=2,k=5 Answer: 175
w=2,h=2,k=6 Answer: 351
...
This is much less obvious to me. If there is a pattern, I'm unsure how to find it.

One more suite of test cases for good measure:

w=2,h=3,k=1 Answer: 1
w=2,h=3,k=2 Answer: 13
w=2,h=3,k=3 Answer: 92
w=2,h=3,k=4 Answer: 430
w=2,h=3,k=5 Answer: 1505
...
yep I'm lost as to what insight if any can be gained.

So, for some arbitrary $w,h,k$, I'm wondering whether there is some type of closed-form formula or perhaps a way to get the answer via dynamic programming/recursion or perhaps divide-and-conquer? I've tried a few of these approaches and didn't really see any obvious breakthrough. The answers get very large very quickly (for example, for $w=2,h=3,k=4$, the answer is 430). I can only imagine how large the answer would be for, say, $w=12,h=12,k=20$.

I would appreciate any help on this. Thanks!
We will use Burnside's lemma, plus some additional optimizations.

Note that row swaps and column swaps commute, so the group of allowable transformations to the grid is exactly $G = S_w \times S_h$, where $S_n$ is the symmetric group on $n$ symbols.  Note that $G$ embeds in $S_{wh}$ via a canonical embedding $\varphi : G \to S_{wh}$, given by $\varphi(\pi,\sigma) = \tau$ where $\tau$ is defined by $\tau(i,j) = (\pi(i), \sigma(j))$.  Therefore, we will think of $G$ simultaneously as the direct product $S_w \times S_h$ and also as a subgroup of $S_{wh}$.

Now let's count the number of grids that are fixed by a group element $g=(\pi,\sigma) \in G$.  Define the cycle type of a permutation $\pi$ to be the multiset of lengths of the cycles in the cycle decomposition of $\pi$.  If $\pi \in S_n$, the cycle type of $\pi$ is a partition of $n$.  With this in mind, we will let $\ell(\pi)$ denote the cycle type of $\pi \in S_w$ and $\ell(\sigma)$ the cycle type of $\sigma \in S_h$.  Finally, let $\ell((\pi,\sigma))$ denote the cycle type of $(\pi,\sigma)$, thinking of $(\pi,\sigma)$ as an element of $S_{wh}$.  We can express $\ell((\pi,\sigma))$ in terms of $\ell(\pi)$ and $\ell(\sigma)$: let $x,y$ range over $x \in \ell(\pi), y \in \ell(\sigma)$, and for each $x,y$, $\ell((\pi,\sigma))$ has $\gcd(x,y)$ copies of $\operatorname{lcm}(x,y)$.

Now the number of grids that are fixed by $(\pi,\sigma)$ is

$$f(\pi,\sigma) = k^{|\ell(\pi,\sigma)|},$$

where $|\ell(\pi,\sigma)|$ counts the number of cycles in $(\pi,\sigma) \in G$ (multiplicities are counted).  By the above characterization of $\ell(\pi,\sigma)$, we see that

$$f(\pi,\sigma) = k^{\sum_{x \in \ell(\pi)} \sum_{y \in \ell(\sigma)} \gcd(x,y)}.$$

To apply Burnside's lemma, it remains to sum $f(g)$ over all elements $g \in G$, as Burnside's lemma tells us that the total number of non-equivalent grids is

$$N = {1 \over |G|} \sum_{g \in G} f(g) = {1 \over w! h!} \sum_{g \in G} f(g).$$

In principle, we could compute this sum by enumerating all $w! h!$ elements in $G$, then evaluating $f$ on each, and summing the result.

It turns out this can be optimized further, as $f((\pi,\sigma))$ depends only on $\ell(\pi)$ and $\ell(\sigma)$, but not on any other aspect of $\pi$ or $\sigma$.  Thus, we can more efficiently evaluate the sum using the formula

$$\sum_{g \in G} f(g) = \sum_{P,Q} n(P,Q) \times k^{\sum_{x \in P} \sum_{y \in Q} \gcd(x,y)},$$

where the inner sum ranges over all partitions $P$ of $w$ and all partitions $Q$ of $h$, and where we define

$$n(P,Q) = |\{(\pi,\sigma) \in G : \ell(\pi)=P, \ell(\sigma)=Q\}|$$

to count the number of group elements of $G$ whose cycle types are $P,Q$.  Now we have

$$n(P,Q) = n_w(P) \times n_h(Q),$$

where $n_w(P) = |\{\pi \in S_w : \ell(\pi)=P\}|$ is the number of permutations in $S_w$ with cycle type $P$, and similarly $n_h(Q)$ is the number of permutations in $S_h$ with cycle type $Q$.  Finally, it is easy to compute $n_w(P)$ for any desired partition $P$; it is

$$n_w(P) = {w! \over \prod_{x \in P} x}.$$

Plugging in, we find that the total number of non-equivalent grids is given by

$$N = \sum_{P,Q} {k^{\sum_{x \in P} \sum_{y \in Q} \gcd(x,y)} \over \prod_{x \in P} x \prod_{y \in Q} y}.$$

This can be computed by two nested for-loops that iterate over all partitions $P$ of $w$ and all partitions $Q$ of $h$.

What is the running time of this algorithm?  The  asymptotic running time is approximately $O(p(w) p(h))$, where $p(n)$ is the number of partitions of the integer $n$.  It is known that asymptotically $p(n) = O(\exp\{\pi \sqrt{2n/3} \})$, so it follows that the asymptotic running time of this algorithm is approximately $O(\exp\{\pi \sqrt{2/3} (\sqrt{w}+\sqrt{h})\})$, or equivalently, approximately $O(2^{3.7(\sqrt{w}+\sqrt{h})})$.  This should allow you to handle the full range of values $1 \le w,h \le 12$, as for that range we have $\sqrt{w}+\sqrt{h} &lt; 7$ and $2^{3.7 \times 7} \approx 2^{26}$ is not too large.