[Add-on 2016-10-30]: Case n=4 added.
We calculate the number $a(n,n)$ of equivalent functions $$f:[n]\times[n]\rightarrow[n]$$ for $n=2,3$ and $n=4$ according to the paper of F. Harary and E. Palmer and show
\begin{align*}
a(2,2)&=7\\
a(3,3)&=638\\
a(4,4)&=7643021
\end{align*}
The formula to be applied is stated as formula (14) in the paper. In fact we can use a simplified version of it, which is given in connection with the calculation of $a(2,2,1)=7$ at the end of page 505. The third parameter is not of interest for us and so we instead write $a(2,2)$.
Here I follow the notation of the authors and use $[p,q]:=\operatorname{lcm}(p,q)$ and $\langle p,q\rangle:=\operatorname{gcd}(p,q)$.
[Harary, Palmer]: The following is valid
\begin{align*}
a(n,n)=\frac{1}{\left(n!\right)^2}\sum_{(\alpha,\beta)\in S_n^2}\prod_{p=1}^n\prod_{q=1}^n
\left(\sum_{s|[p,q]}sj_s(\alpha)\right)^{j_p(\alpha)j_q(\beta)\langle p,q\rangle}\tag{1}
\end{align*}
where the sum is taken over all pairs of permutations $(\alpha,\beta)$ of degree $n$ and $j_p(\alpha)$ is denoting the number of cycles of $\alpha$ of length $p$.
Hint: Observe the terms in the sum (1) do not make use of $\alpha$ but instead of $j_p(\alpha)$ only. So, it is not necessary to sum over all $\left(n!\right)^2$ pairs of permutations, as we can conveniently use the cycle index of the permutation group $S_n$ and considerably reduce the number of summands.
Preparatory work: Cycle index
We calculate for $n=2,3,4$ the cycle index based upon the recursion formula
\begin{align*}
Z(S_0)&=1\\
Z(S_n)&=\frac{1}{n}\sum_{j=1}^nz_jZ(S_{n-j})\qquad\qquad n>0
\end{align*}
We obtain
\begin{align*}
Z(S_1)&=z_1Z(S_0)=z_1\\
Z(S_2)&=\frac{1}{2}\left(z_1Z(S_1)+z_2\right)\\
&=\frac{1}{2}\left(z_1^2+z_2\right)\\
Z(S_3)&=\frac{1}{3}\left(z_1\cdot\frac{1}{2}\left(z_1^2+z_2\right)+z_2z_1+z_3\right)\\
&=\frac{1}{6}\left(z_1^3+3z_1z_2+2z_3\right)\\
Z(S_4)&=\frac{1}{4}\left(z_1\cdot\frac{1}{6}\left(z_1^3+3z_1z+2z_3\right)+z_2\cdot\frac{1}{2}\left(z_1^2+z_2\right)
+z_3z_1+z_4\right)\\
&=\frac{1}{24}\left(z_1^4+6z_1^2z_2+8z_1z_3+3z_2^2+6z_4\right)
\end{align*}
$$ $$
Case: $n=2$:
We show the following is valid
\begin{align*}
\color{blue}{a(2,2)=7}
\end{align*}
In order to calculate $a(2,2)$ we consider according to (1)
\begin{align*}
a(2,2)=\frac{1}{4}\sum_{(\alpha,\beta)\in S_2^2}\prod_{p=1}^2\prod_{q=1}^2
\left(\sum_{s|[p,q]}sj_s(\alpha)\right)^{j_p(\alpha)j_q(\beta)\langle p,q\rangle}\tag{2}
\end{align*}
It it convenient to do some bookkeeping by use of tables. We list the permutations of $S_2=\{\operatorname{id},(12)\}$ in cycle notation and write a table with the number of cycles of each length for each permutation. We also write the corresponding monomial from the cycle index.
\begin{array}{l|ccc}
\pi&Z(S_2)&j_1(\pi)&j_2(\pi)\\
\hline
\operatorname{id}& z_1^2& 2& 0\\
(12)& z_2^1& 0 &1\\
\end{array}
Since it is too cumbersome to write each summand from (2) in one long line we use instead a table description as follows:
\begin{array}{cc|cccc|cc|cc|rr}
\alpha&\beta&p&q&[p,q]&<p,q>&s&j_s(\alpha)&j_p(\alpha)&j_q(\beta)&\text{factors}&\text{result}\\
\hline
id&id&1&1&1&1&1&2&2&2&16&16\\
&&1&2&2&1&1&2&2&0&1&\\
&&&&&&2&0&&&&\\
&&2&1&2&1&1&2&0&0&1&\\
&&&&&&2&0&&&&\\
&&2&2&2&2&1&2&0&2&1&\\
&&&&&&2&0&&&&\\
\hline
id&(12)&1&1&1&1&1&2&2&0&1&4\\
&&1&2&2&1&1&2&2&1&4&\\
&&&&&&2&0&&&&\\
&&2&1&2&1&1&2&0&0&1&\\
&&&&&&2&0&&&&\\
&&2&2&2&2&1&2&0&1&1&\\
&&&&&&2&0&&&&\\
\hline
(12)&id&1&1&1&1&1&0&0&2&1&4\\
&&1&2&2&1&1&0&0&0&1&\\
&&&&&&2&1&0&0&&\\
&&2&1&2&1&1&0&1&2&4&\\
&&&&&&2&1&1&2&&\\
&&2&2&2&2&1&0&1&0&1&\\
&&&&&&2&1&1&0&&\\
\hline
(12)&(12)&1&1&1&1&1&0&0&0&1&4\\
&&1&2&2&1&1&0&0&1&1&\\
&&&&&&2&1&0&1&&\\
&&2&1&2&1&1&0&1&0&1&\\
&&&&&&2&1&1&0&&\\
&&2&2&2&2&1&0&1&1&4&\\
&&&&&&2&1&1&1&&\\
\hline
&\color{blue}{\text{Total}}&&&&&&&&&&\color{blue}{28}
\end{array}
Comment:
The table is organised in blocks for pairs of permutation. Although here not eye-catching since we list all $\left(2!\right)^2=4$ pairs, we need in fact only for each cycle type one representative, since we are only interested in the length of cycles of a permutation. The column result gives the summands in (2). Here are the gory details:
Columns: $\alpha,\beta$ correspond to a pair of permutations $(\alpha,\beta)$ which is used as index in the outer sum of (2).
Columns: $p,q$ are the indices of the products in (2)
Column: $s$ gives the divisors of $\operatorname{lcm}(p,q)$
Columns: $j_s(\alpha),j_p(\alpha),j_q(\beta)$ list the cycle lengths
Column: $\text{factor}$ gives
$$\left(\sum_{s|[p,q]}sj_s(\alpha)\right)^{j_p(\alpha)j_q(\beta)\langle p,q\rangle}$$
Column: $\text{result}$ calculates finally the product
$$\prod_{p=1}^2\prod_{q=1}^2
\left(\sum_{s|[p,q]}sj_s(\alpha)\right)^{j_p(\alpha)j_q(\beta)\langle p,q\rangle}$$
Since the total of the table is $28$ we finally conclude according to (2)
\begin{align*}
\color{blue}{a(2,2)=\frac{1}{4}\cdot 28=7}
\end{align*}
and the claim follows.
$$ $$
Case: $n=3$:
We do the calculation similar to above and show the following is valid
\begin{align*}
\color{blue}{a(3,3)=638}
\end{align*}
In order to calculate $a(3,3)$ we consider according to (1)
\begin{align*}
a(3,3)=\frac{1}{\left(3!\right)^2}\sum_{(\alpha,\beta)\in S_3^2}\prod_{p=1}^3\prod_{q=1}^3
\left(\sum_{s|[p,q]}sj_s(\alpha)\right)^{j_p(\alpha)j_q(\beta)\langle p,q\rangle}\tag{3}
\end{align*}
We list the permutations of $S_3=\{\operatorname{id},(12),(13),(23),(123),(132)\}$ in cycle notation and write a table with the number of cycles of each length for each permutation. We also write the corresponding monomial from the cycle index.
\begin{array}{l|cccc}
\pi&Z(S_3)&j_1(\pi)&j_2(\pi)&j_3(\pi)\\
\hline
id&z_1^3&3&0&0\\
(12)&3z_1z_2&1&1&0\\
(123)&2z_3&0&0&1\\
\end{array}
Note: The factors $1,3$ and $2$ in the column $Z(S_3)$ indicate the number of different permutations of the corresponding cycle type. We will use this fact to considerably reduce the calculation of the number of summands in (3).
In the following it is sufficient to calculate tables for the nine pairs
\begin{align*}
\{id,(12),(123)\}\times\{id,(12),(123)\}
\end{align*}
the cycle index provides the supplementary information we need to calculate the complete sum.
Note that in the main table above there is some redundancy to ease traceability. We now use a somewhat more compact notation to ease readability and keep the space small.
Table: $j_s(\pi), j_p(\pi),j_q(\pi)$
\begin{array}{cc|cc|ccc|ccc|ccc}
&&&&&\pi=id&&&\pi=(12)&&&\pi=(123)&\\
p&q&[p,q]&s&j_s&j_p&j_q&j_s&j_p&j_q&j_s&j_p&j_q\\
\hline
1&1&1&1&3&3&3&1&1&1&0&0&0\\
&2&2&1&3&3&0&1&1&1&0&0&0\\
&&&2&0&&&1&&&0&&\\
&3&3&1&3&3&0&1&1&0&0&0&1\\
&&&3&0&&&0&&&1&&\\
2&1&2&1&3&0&3&1&1&1&0&0&0\\
&&&2&0&&&1&&&0&&\\
&2&2&1&3&0&0&1&1&1&0&0&0\\
&&&2&0&&&1&&&0&&\\
&3&6&1&3&0&0&1&1&0&0&0&1\\
&&&2&0&&&1&&&0&&\\
&&&3&0&&&0&&&1&&\\
3&1&3&1&3&0&3&1&0&1&0&1&0\\
&&&3&0&&&0&&&1&&\\
&2&6&1&3&0&0&1&0&1&0&1&0\\
&&&2&0&&&1&&&0&&\\
&&&3&0&&&0&&&1&&\\
&3&3&1&3&0&0&1&0&0&0&1&1\\
&&&3&0&&&0&&&1&&\\
\end{array}
The table above provides all information necessary to calculate the summands in (3) for each of the nine pairs of permutations.
An example of a typical block is given here for $((12),(12))$ as it was done for all four blocks in the case $n=2$ and a summary table follows below.
Table: $\{(12)\}\times\{(12)\}$
\begin{array}{cc|cccc|cc|cc|rr}
\alpha&\beta&p&q&[p,q]&<p,q>&s&j_s(\alpha)&j_p(\alpha)&j_q(\beta)&\text{factors}&\text{result}\\
\hline
(12)&(12)&1&1&1&1&1&1&1&1&1&81\\
&&&2&2&1&1&1&1&1&3&\\
&&&&&&2&1&&&&\\
&&&3&3&1&1&1&1&0&1&\\
&&&&&&3&0&&&&\\
&&2&1&2&1&1&1&1&1&3&\\
&&&&&&2&1&&&&\\
&&&2&2&2&1&1&1&1&9&\\
&&&&&&2&1&&&&\\
&&&3&6&1&1&1&1&0&1&\\
&&&&&&2&1&&&&\\
&&&&&&3&0&&&&\\
&&3&1&3&1&1&1&0&1&1&\\
&&&&&&3&0&&&&\\
&&&2&6&1&1&1&0&1&1&\\
&&&&&&2&1&&&&\\
&&&&&&3&0&&&&\\
&&&3&3&3&1&1&0&0&1&\\
&&&&&&3&0&&&&\\
\end{array}
$$ $$
Summary:
In order to respect all summands of (3) we write the results of the tables above together with the multiplicity of each permutation according to its cycle type. So, e.g. the permutation $(12)$ has cycle type $z_1z_2$ and there are three permutations of this type $\{(12),(13),(23)\}$, we take a factor $3$.
\begin{array}{ll|r|cc|r}
\alpha&\beta&\text{res}&\text{m}_{\alpha}&\text{m}_{\beta}&\text{res}\cdot \text{m}_{\alpha}\cdot \text{m}_{\beta}\\
\hline
id&id&19683&1&1&19683\\
id&(12)&729&1&3&2187\\
id&(123)&27&1&2&54\\
\hline
(12)&id&27&3&1&81\\
(12)&(12)&81&3&3&729\\
(12)&(123)&3&3&2&18\\
\hline
(123)&id&27&2&1&54\\
(123)&(12)&9&2&3&54\\
(123)&(123)&27&2&2&108\\
\hline
\color{blue}{\text{Total}}&&&&&\color{blue}{22968}\\
\end{array}
Since the total of the table is $22968$ we finally conclude according to (3)
\begin{align*}
\color{blue}{a(3,3)=\frac{1}{36}\cdot 22968=638}
\end{align*}
and the claim follows.
$$ $$
Case: $n=4$:
We do the calculation similar to above and show the following is valid
\begin{align*}
\color{blue}{a(4,4)=7643021}
\end{align*}
In order to calculate $a(4,4)$ we consider according to (1)
\begin{align*}
a(4,4)=\frac{1}{\left(4!\right)^2}\sum_{(\alpha,\beta)\in S_4^2}\prod_{p=1}^4\prod_{q=1}^4
\left(\sum_{s|[p,q]}sj_s(\alpha)\right)^{j_p(\alpha)j_q(\beta)\langle p,q\rangle}\tag{4}
\end{align*}
We list the permutations of $S_4$ in cycle notation and write a table with the number of cycles of each length for each permutation. We also write the corresponding monomial from the cycle index.
\begin{array}{l|ccccc}
\pi&Z(S_4)&j_1(\pi)&j_2(\pi)&j_3(\pi)&j_3(\pi)\\
\hline
id&z_1^4&4&0&0&0\\
(12)&6z_1^2z_2&1&1&0&0\\
(123)&8z_1z_3&0&0&1&0\\
(12)(34)&3z_2^2&0&2&0&0\\
(1234)&6z_4&0&0&0&1\\
\end{array}
Note: The factors $1,6,8,3$ and $6$ in the column $Z(S_4)$ indicate the number of different permutations of the corresponding cycle type. We will use this fact to considerably reduce the calculation of the number of summands in (4).
In the following it is sufficient to calculate tables for the $25$ pairs
\begin{align*}
\{id,(12),(123),(12)(34),(1234)\}\times\{id,(12),(123),(12)(34),(1234)\}
\end{align*}
the cycle index provides the supplementary information we need to calculate the complete sum.
Note that in the main table of $n=2$ above there is some redundancy to ease traceability. We now use analogously to $n=3$ above a somewhat more compact notation to ease readability and keep the space small.
Table: $j_s(\pi), j_p(\pi),j_q(\pi)$
\begin{array}{cc|cc|ccc|ccc|ccc}
&&&&&\pi=id&&&\pi=(12)&&&\pi=(123)&\\
p&q&[p,q]&s&j_s&j_p&j_q&j_s&j_p&j_q&j_s&j_p&j_q\\
\hline
1&1&1&1&4&4&4&2&2&2&1&1&1\\
1&2&2&1&4&4&0&2&2&1&1&1&0\\
&&&2&0&&&1&&&0&&\\
1&3&3&1&4&4&0&2&2&0&1&1&1\\
&&&3&0&&&0&&&1&&\\
1&4&4&1&4&4&0&2&2&0&1&1&0\\
&&&2&0&&&1&&&0&&\\
&&&4&0&&&0&&&0&&\\
2&1&2&1&4&0&4&2&1&2&1&0&1\\
&&&2&0&&&1&&&0&&\\
2&2&2&1&4&0&0&2&1&1&1&0&0\\
&&&2&0&&&1&&&0&&\\
2&3&6&1&4&0&0&2&1&0&1&0&1\\
&&&2&0&&&1&&&0&&\\
&&&3&0&&&0&&&1&&\\
2&4&4&1&4&0&0&2&1&0&1&0&0\\
&&&2&0&&&1&&&0&&\\
&&&4&0&&&0&&&0&&\\
3&1&3&1&4&0&4&2&0&2&1&1&1\\
&&&3&0&&&0&&&1&&\\
3&2&6&1&4&0&0&2&0&1&1&1&0\\
&&&2&0&&&1&&&0&&\\
&&&3&0&&&0&&&1&&\\
3&3&3&1&4&0&0&2&0&0&1&1&1\\
&&&3&0&&&0&&&1&&\\
3&4&12&1&4&0&0&2&0&0&1&1&0\\
&&&2&0&&&1&&&0&&\\
&&&3&0&&&0&&&1&&\\
&&&4&0&&&0&&&0&&\\
4&1&4&1&4&0&4&2&0&2&1&0&1\\
&&&2&0&&&1&&&0&&\\
&&&4&0&&&0&&&0&&\\
4&2&4&1&4&0&0&2&0&1&1&0&0\\
&&&2&0&&&1&&&0&&\\
&&&4&0&&&0&&&0&&\\
4&3&12&1&4&0&0&2&0&0&1&0&1\\
&&&2&0&&&1&&&0&&\\
&&&3&0&&&0&&&1&&\\
&&&4&0&&&0&&&0&&\\
4&4&4&1&4&0&0&2&0&0&1&0&0\\
&&&2&0&&&1&&&0&&\\
&&&4&0&&&0&&&0&&\\
&&&&&&&&&&&&\\
&&&&&&&&&&&&\\
\end{array}
$$ $$
Table (cont.): $j_s(\pi), j_p(\pi),j_q(\pi)$
\begin{array}{cc|cc|ccc|cccccc}
&&&&&\pi=(12)(34)&&&\pi=(1234)&\\
p&q&[p,q]&s&j_s&j_p&j_q&j_s&j_p&j_q\\
\hline
1&1&1&1&0&0&0&0&0&0&&&\\
1&2&2&1&0&0&2&0&0&0&&&\\
&&&2&2&&&0&&&&&\\
1&3&3&1&0&0&0&0&0&0&&&\\
&&&3&0&&&0&&&&&\\
1&4&4&1&0&0&0&0&0&1&&&\\
&&&2&2&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
2&1&2&1&0&2&0&0&0&0&&&\\
&&&2&2&&&0&&&&&\\
2&2&2&1&0&2&2&0&0&0&&&\\
&&&2&2&&&0&&&&&\\
2&3&6&1&0&2&0&0&0&0&&&\\
&&&2&2&&&0&&&&&\\
&&&3&0&&&0&&&&&\\
2&4&4&1&0&2&0&0&0&1&&&\\
&&&2&2&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
3&1&3&1&0&0&0&0&0&0&&&\\
&&&3&0&&&0&&&&&\\
3&2&6&1&0&0&2&0&0&0&&&\\
&&&2&2&&&0&&&&&\\
&&&3&0&&&0&&&&&\\
3&3&3&1&0&0&0&0&0&0&&&\\
&&&3&0&&&0&&&&&\\
3&4&12&1&0&0&0&0&0&1&&&\\
&&&2&2&&&0&&&&&\\
&&&3&0&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
4&1&4&1&0&0&0&0&1&0&&&\\
&&&2&2&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
4&2&4&1&0&0&2&0&1&0&&&\\
&&&2&2&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
4&3&12&1&0&0&0&0&1&0&&&\\
&&&2&2&&&0&&&&&\\
&&&3&0&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
4&4&4&1&0&0&0&0&1&1&&&\\
&&&2&2&&&0&&&&&\\
&&&4&0&&&1&&&&&\\
\end{array}
The table above provides all information necessary to calculate the summands in (4) for each of the $25$ pairs of permutations.
Hint: Observe, that we only need to consider $25$ pairs of permutations instead of $\left(4!\right)^2=576$ pairs which are summed up in (4). We will consider all other permutations by respecting multiplicities given by the cycle-index $Z(S_4)$.
An example of a typical block is given here for $((123),(123))$ as it was done for all four blocks in the case $n=2$ and a summary table follows below.
Table: $\{(123)\}\times\{(123)\}$
\begin{array}{cc|cccc|cc|cc|rr}
\alpha&\beta&p&q&[p,q]&<p,q>&s&j_s(\alpha)&j_p(\alpha)&j_q(\beta)&\text{factors}&\text{result}\\
\hline
(123)&(123)&1&1&1&1&1&1&1&1&1&1024\\
&&1&2&2&1&1&1&1&0&1&\\
&&&&&&2&0&&&&\\
&&1&3&3&1&1&1&1&1&4&\\
&&&&&&3&1&&&&\\
&&1&4&4&1&1&1&1&0&1&\\
&&&&&&2&0&&&&\\
&&&&&&4&0&&&&\\
&&2&1&2&1&1&1&0&1&1&\\
&&&&&&2&0&&&&\\
&&2&2&2&2&1&1&0&0&1&\\
&&&&&&2&0&&&&\\
&&2&3&6&1&1&1&0&1&1&\\
&&&&&&2&0&&&&\\
&&&&&&3&1&&&&\\
&&2&4&4&2&1&1&0&0&1&\\
&&&&&&2&0&&&&\\
&&&&&&4&0&&&&\\
&&3&1&3&1&1&1&1&1&4&\\
&&&&&&3&1&&&&\\
&&3&2&6&1&1&1&1&0&1&\\
&&&&&&2&0&&&&\\
&&&&&&3&1&&&&\\
&&3&3&3&3&1&1&1&1&64&\\
&&&&&&3&1&&&&\\
&&3&4&12&1&1&1&1&0&1&\\
&&&&&&2&0&&&&\\
&&&&&&3&1&&&&\\
&&&&&&4&0&&&&\\
&&4&1&4&1&1&1&0&1&1&\\
&&&&&&2&0&&&&\\
&&&&&&4&0&&&&\\
&&4&2&4&2&1&1&0&0&1&\\
&&&&&&2&0&&&&\\
&&&&&&4&0&&&&\\
&&4&3&12&1&1&1&0&1&1&\\
&&&&&&2&0&&&&\\
&&&&&&3&1&&&&\\
&&&&&&4&0&&&&\\
&&4&4&4&4&1&1&0&0&1&\\
&&&&&&2&0&&&&\\
&&&&&&4&0&&&&\\
\end{array}
$$ $$
Summary:
In order to respect all summands of (4) we write the results of the tables above together with the multiplicity of each permutation according to its cycle type. So, e.g. the permutation $(12)$ has cycle type $z_1z_2$ and there are three permutations of this type $\{(12),(13),(23)\}$, we take a factor $3$.
\begin{array}{ll|r|cc|r}
\alpha&\beta&\text{res}&\text{m}_{\alpha}&\text{m}_{\beta}&\text{res}\cdot \text{m}_{\alpha}\cdot \text{m}_{\beta}\\
\hline
id&id&4294967296&1&1&4294967296\\
id&(12)&16777216&1&6&100663296\\
id&(123)&65536&1&8&524288\\
id&(12)(34)&65536&1&3&196608\\
id&(1234)&256&1&6&1536\\
\hline
(12)&id&65536&6&1&393216\\
(12)&(12)&65536&6&6&2359296\\
(12)&(123)&256&6&8&12288\\
(12)&(12)(34)&65536&6&3&1179648\\
(12)&(1234)&256&6&6&9216\\
\hline
(123)&id&256&8&1&2048\\
(123)&(12)&64&8&6&3072\\
(123)&(123)&1024&8&8&65536\\
(123)&(12)(34)&16&8&3&384\\
(123)&(1234)&4&8&6&192\\
\hline
(12)(13)&id&65536&3&1&196608\\
(12)(13)&(12)&65536&3&6&1179648\\
(12)(13)&(123)&256&3&8&6144\\
(12)(13)&(12)(34)&65536&3&3&589824\\
(12)(13)&(1234)&256&3&6&4608\\
\hline
(1234)&id&256&6&1&1536\\
(1234)&(12)&256&6&6&9216\\
(1234)&(123)&16&6&8&768\\
(1234)&(12)(34)&256&6&3&4608\\
(1234)&(1234)&256&6&6&9216\\
\hline
\color{blue}{\text{Total}}&&&&&\color{blue}{4402380096}\\
\end{array}
Since the total of the table is $4402380096$ we finally conclude according to (4)
\begin{align*}
\color{blue}{a(4,4)=\frac{1}{576}\cdot 4402380096=7643021}
\end{align*}
and the claim follows.
Conclusion:
- In order to calculate $a(n,n)$ we do not need $(n!)^2$ summands but can calculate summands corresponding to the square of the number of summands of the cycle index and then use multiplicities of the cycle index for final calculations.
\begin{array}{c|rr}
n&(n!)^2&\left(\text{via }Z(S_n)\right)^2\\
\hline
2&4&4\\
3&36&9\\
4&576&25\\
\end{array}
- It seems feasible to find an efficient implemention based upon this formula. A nice program which coincides with the results of this answer is already given by @ScottBurns.
