Counting on a matrix
Copied from here.
The mathematics
Consider the group  $G=S_w\times S_h$, where $S_w$ and $S_h$ are the symmetric group of the sets $W=\{1,2,3,\ldots,w\}$ and $H=\{1,2,3,\ldots,h\}$ with $w$ and $h$ elements, respectively. The group $G$ acts on the set $X=W\times H$, which we view as the set of indexes of the entries of the matrices.
Each matrix is a function $f:X\to S$, and element $f\in S^X$.
So, the problem consists of computing the cardinality
$$\left|S^X/G\right|$$
Cycle index polynomial of a group of permutations
For a permutation group $G$ acting on a set $X=\{1,2,3,\ldots,n\}$ define the polynomial
$$Z(G,s_1,s_2,s_3,\ldots,s_n)=|G|^{-1}\sum_{g\in G}\prod_{k=1}^{n}s_k^{c_k(g)}$$
where $c_k(g)$ is the number of cycles of length $k$ in the cycle decomposition of the permutation $g$.
Polya's enumeration theorem tell us that
$$\left|S^X/G\right|=Z(G,s,s,\ldots,s)=|G|^{-1}\sum_{g\in G}s^{c_k(g)}$$
Cycle index of the Cartesian product
In our case, $G$ is the Cartesian product $S_w\times S_h$. In this paper they prove the following formula
$$Z(G_1\times G_2,s_1,s_2,\ldots,s_{n_1\cdot n_2}) = Z(G_1,s_1,\ldots,s_{n_1})\otimes Z(G_1,s_1,s_2,\ldots,s_{n_2})$$
where the product $\otimes$ (In the paper they use the symbol text reference mark) of polynomials is defined by:
If
$$f(x_1,x_2,...,x_m)=\sum a_{i_1i_2\ldots i_m}x_1^{i_1}x_2^{i_2}\dotsm x_m^{i_m}\text{ and }g(x_1,x_2,...,x_n)=\sum b_{j_1j_2\ldots j_n}x_1^{j_1}x_2^{j_2}\dotsm x_n^{j_n}$$
then
$$f(x_1,x_2,\ldots,x_m)\otimes g(x_1,x_2,\ldots,x_n)=\sum a_{i_1i_2\ldots i_m}b_{j_1j_2\ldots j_n}\times\prod_{\substack{1\leq r\leq m\\1\leq s\leq n}}(x_{r}^{i_l}\otimes x_{s}^{j_s})$$
where
$$x_{r}^{i_l}\otimes x_{s}^{j_s} = x_{\operatorname{lcm}(r,s)}^{\gcd(r,s)i_rj_s}$$
Cycle index of the symmetric group
For the case of the symmetric group $S_n$ acting on $\\{1,2,3,\ldots,n\\}$ we have
$$Z(S_n)=\frac{1}{n}\sum_{k=1}^{n}s_kZ(S_{n-k})$$
or explicitly
$$Z(S_n,s_1,s_2,\ldots,s_n)=\sum_{r_1+2r_2+3r_3+\ldots nr_n=n}\frac{s_1^{r_1}s_2^{r_2}\dotsm s_n^{r_n}}{1^{r_1}r_1!2^{r_2}r_2!\dotsm n^{r_n}r_n!}$$
Solution of the problem
Combining the results above we obtain that the number of classes of $w\times h$ matrices on $s$ symbols up to permutations of the rows and columns is
$$N=\frac{1}{w!h!}\sum_{\substack{i\in P(W)\\j\in P(H)}}\frac{w!}{1^{i_1}i_1!2^{i_2}i_2!\dotsm w^{i_w}i_w!}\frac{h!}{1^{j_1}j_1!2^{j_2}j_2!\dotsm h^{j_h}j_h!}s^{\sum_{\substack{a\in i\\b\in j}}\gcd(a,b)}$$
where the outer sums runs over the partitions $i=(i_1,i_2,\ldots,i_w)$ and $j=(j_1,j_2,\ldots,j_h)$ of $w$ and $h$, respectively. This is, $i_1+2i_2+\ldots+wi_w=w$ and $j_1+2j_2+\ldots+hj_h=h$. The inner sums run over the elements $a$ and $b$ of the partitions $i$ and $j$. This is, $a$ takes the value $1$, $i_1$ times, the value $2$ $i_2$ times, and so on and likewise $b$ takes the value $1$ $j_1$ times, the value $2$ $j_2$ times, etc.
Implementation
As seen in the final formula, we need to compute:

Partitions of a number
Greatest common divisors (gcd) of many numbers.
Factorials of many numbers.

For these, we can do:

To compute all partitions one can use the iterative algorithms here.
An efficient way to compute gcd one could use Euclidean algorithm. However, since we are going to need the gcd of all pairs of numbers in a range and each one many times. It is better to pre-compute the full table of gcd all at once by dynamic programming. If a&gt;b, then gcd(a,b)=gcd(a-b,b). This recurrence equation allows to compute gcd of larger pairs in terms of that of smaller pairs. In the table, one has the initial values gcd(1,a)=gcd(a,1)=1 and gcd(a,a)=a, for all a.
The same happens for factorials. The formula will require the factorials of all numbers in a range many times each. So, it is better to compute them all from the bottom up using that n! = n(n-1)! and 0!=1!=1.
The powers of $s$, which are also needed repeatedly in the sum could also be pre-computed.

An implementation in Python could look like this. Feel free to improve it.