On Computing the Number of Latin Rectangles

Abstract

Doyle (circa 1980) found a formula for the number of $k \times n$ Latin rectangles $L_{k,n}$. This formula remained dormant until it was recently used for counting $k \times n$ Latin rectangles, where $k \in \{4,5,6\}$. We give a formal proof of Doyle’s formula for arbitrary k. We also improve a previous implementation of this formula, which we use to find $L_{k,n}$ when $k=4$ and $n \le 150$, when $k=5$ and $n \le 40$ and when $k=6$ and $n \le 15$. Motivated by computational data for $3 \le k \le 6$, some research problems and conjectures about the divisors of $L_{k,n}$ are presented.

1 Introduction

A $k \times n$ Latin rectangle is a $k \times n$ matrix L, with symbols from $\{1,2,\ldots ,n\}$, such that each row and each column contains only distinct symbols. If $k=n$ then L is a Latin square of order n. Let $L_{k,n}$ be the number of $k \times n$ Latin rectangles.

The main purpose of this article is to report the computation of previously unknown values of $L_{k,n}$ (for $4 \le k \le 6$). To the end, we use a formula by Doyle, for which we give a proof for in Sect. 3. We use these values to provide evidentiary support for a question we pose regarding the divisors of $L_{k,n}$ (Question 4). We also present conjectures about divisors of $K_{3,n}$ (Conjecture 5). Questions about divisors of $L_{n,n}$ were previously raised by Alter [1] and about divisors of $L_{k,n}$ in [35].

A Latin rectangle is called normalised if the first row is $(1,2,\dots ,n)$, and reduced if the first row is $(1,2,\dots ,n)$ and the first column is $(1,2,\dots ,k)^T$. Let $K_{k,n}$ denote the number of normalised $k \times n$ Latin rectangles and let $R_{k,n}$ denote the number of reduced $k \times n$ Latin rectangles. In the case of Latin squares, the numbers $L_{n,n}$, $K_{n,n}$ and $R_{n,n}$ will be denoted $L_n$, $K_n$ and $R_n$, respectively. The three numbers $L_{k,n}$, $K_{k,n}$ and $R_{k,n}$ are related by

$$\begin{aligned} L_{k,n}=n! K_{k,n} = \frac{n! (n-1)!}{(n-k)!} R_{k,n}. \end{aligned}$$

In particular

$$\begin{aligned} L_n=n! K_n = n! (n-1)! R_n. \end{aligned}$$

So finding the value of $L_{k,n}$ is essentially the same as finding the value of $R_{k,n}$ or $K_{k,n}$. McKay and Wanless [23] published a table of values for $R_{k,n}$ when $2 \le k \le n \le 11$, which were obtained by lengthy computer enumerations. Some values of $R_{k,n}$ for $k \in \{4,5,6\}$ were reported in [36]. Also see [36] for a survey of the formulae involving the number of Latin rectangles.

Gessel [15] proved that for every $k \ge 1$, there exists a finite $M=M(k)$ such that there exists $M+1$ polynomials $c_i(n)$ such that

$$\begin{aligned} \sum _{0 \le i \le M} c_i(n) L_{k,n+i}=0 \end{aligned}$$

for all $n \ge k$ (other than when each $c_i(n)=0$). So recurrence relations theoretically exist for $L_{k,n}$ for any fixed k. For example, we know that

$$\begin{aligned} (n+1)L_{1,n}-L_{1,n+1}=0 \end{aligned}$$

since $L_{1,n}=n!$ and

$$\begin{aligned} (n+2)(n+1)^2 L_{2,n}+(n+2)(n+1) L_{2,n+1}-L_{2,n+2}=0 \end{aligned}$$

since $L_{2,n}=n! d_n$ where $d_n$ is the number of derangements of n elements, and $d_n$ satisfies the well-known recurrence $d_n=(n-1)(d_{n-1}+d_{n-2})$ where $d_0=1$ and $d_1=0$. Theorem 1 (in the next section) gives such a recurrence for $k=3$.

In the next section we give an historical survey of the formulae for $R_{3,n}$, the case of three-line Latin rectangles. Despite the numerous formulae for $R_{3,n}$, we are still unable to fully explain the growth in the largest power of 2 that divides $R_{3,n}$. Afterwards, we use a formula by Doyle (improving on the implementation used in [35]) to compute values of $R_{k,n}$ when $4 \le k \le 6$. We give a proof of this formula in Sect. 3. We find that the divisors of $R_{k,n}$ with $k \in \{4,5,6\}$ display similar behaviour to that of $R_{3,n}$. Based on this data, in Sect. 5 we discuss the asymptotic behaviour of the divisors of $R_{k,n}$ for fixed k and present some conjectures about the divisors of $R_{3,n}$.

2 Three-line Latin Rectangles

The enumeration of three-line Latin rectangles has a long history, which we will now review in detail. We also direct the reader to [36], which gives a survey of formulae for $L_{k,n}$ for general k, but does not go into detail about the three-line case. We list the first few non-zero values of $L_{3,n}$, $K_{3,n}$ and $R_{3,n}$ in Table  1.

Table 1 $L_{3,n}$, $K_{3,n}$ and $R_{3,n}$ for small n

There are many known general formulae for $L_{k,n}$ (see [36]), of which $L_{3,n}$ is a special case. For example, the following result originated with MacMahon [22] in 1898. Let $X=(x_{ij})$ be a $k \times n$ matrix whose symbols are the kn variables $x_{ij}$. We index the rows of X by $[k]:=\{1,2,\ldots ,k\}$ and the columns of X by $[n]:=\{1,2,\ldots ,n\}$. Let ${\mathcal {S}}_{k,n}$ be the set of injections $\sigma : [k] \rightarrow [n]$. We define the permanent of the rectangular matrix X to be

$$\begin{aligned} \mathrm {per}(X)=\sum _{\sigma \in {\mathcal {S}}_{k,n}} \prod _{i=1}^{k} x_{i\sigma (i)} . \end{aligned}$$

Then $L_{k,n}$ is the coefficient of $\prod _{i=1}^{k} \prod _{j=1}^{n} x_{ij}$ in $\mathrm {per}(X)^n$. Let ${\mathcal {B}}_{k,n}$ be the set of $k \times n$ $\,(0,1)$-matrices. Fu [13] and Shao and Wei [31] gave the formula

$$\begin{aligned} L_{k,n}=\sum _{A \in {\mathcal {B}}_{k,n}} (-1)^{\sigma _0(A)} \mathrm {per}(A)^n, \end{aligned}$$

where $\sigma _0(A)$ is the number of 0 elements in A, which was generalised in [36] (see also [23]).

The first published formula specifically for $L_{3,n}$ seems to be by Jacob [19] in 1930, although some errors in [19] were later rectified by Kerawala [20] to give the following theorem.

Theorem 1

For any $n \ge 3$,

$$\begin{aligned} \sum _{0 \le i \le 5} c_i(n) K_{3,n+i}=0 \end{aligned}$$

where

$$\begin{aligned} c_0(n)=&\,-4(n+1)(n+2)(n+3)(n+4)^2,\\ c_1(n)=&\,2(n+2)(n+3)(n+4)(n^2+5n+3),\\ c_2(n)=&\,(n+3)(n+4)(n^2+8n+13),\\ c_3(n)=&\,(n+3)(n+4)(n^2+8n+17),\\ c_4(n)=&\,(n+4)(n^2+8n+17),\\ c_5(n)=&\,-(n+3). \end{aligned}$$

Riordan [26] proved that $K_{3,n} \sim n!^2 \exp (-3)$, which he described as “a result which Kerawala surmised but failed to prove, though his numerical evidence was practically conclusive (agreeing with $\exp (-3)$ to seven decimal places).” Jacob [19] also mentioned the possibility of $n!^2/K_{3,n} \rightarrow \exp (3)$, but it was subsequently abandoned. This asymptotic formula was generalised by Erdős and Kaplansky [11], who found that

$$\begin{aligned} L_{k,n} \sim n!^k \exp \big (-k(k-1)/2\big ) \end{aligned}$$

for $k=O\big ((\log n)^{3/2-\epsilon }\big )$, which has since been extended (see [36] for an history). The most up-to-date asymptotic enumeration is by Godsil and McKay [17], who proved

$$\begin{aligned} L_{k,n} \sim n!^k \left( n(n-1)\cdots (n-k+1)/n^k\right) ^n \left( 1-k/n\right) ^{-n/2}\exp (-k/2) \end{aligned}$$

as $n \rightarrow \infty$ with $k = o(n^{6/7})$.

In [27] (see also [29]), Riordan gave a formula for $K_{3,n}$ that involves the problème de ménages. Specifically

$$\begin{aligned} K_{3,n}=\sum _{i=0}^{\lfloor n/2\rfloor }\left( {\begin{array}{c}n\\ i\end{array}}\right) d_i d_{n-i} u_{n-2i} \end{aligned}$$

(1)

where $d_n$ is the number of derangements of n elements (Sloane’s [34] A000166) and $u_n$ is given by

$$\begin{aligned} u_n=\sum _{i=0}^n (-1)^i \frac{2n}{2n-i} \left( {\begin{array}{c}2n-i\\ i\end{array}}\right) (n-i)!, \end{aligned}$$

(2)

where $u_0=1$. The numbers $u_n$ for $n \ge 2$ are the ménages numbers (Sloane’s A000179). The ménages numbers can also be defined as the number of permutations $\sigma$ of [n] such that $\sigma (i) \not \equiv i \pmod n$ and $\sigma (i) \not \equiv i+1 \pmod n$ for all $i \in [n]$. However, for (1) to be valid, we require $u_1=-1$, which is inconsistent with the “number of $\sigma$” definition, but is consistent with (2) (this issue was raised by Vladimir Shevelev on Sloane’s A000186). Riordan’s formula (1) was generalised by Moser [24] to count normalised three-line Latin rectangles where the derangement defined by the second row does not have cycle lengths belonging to some set S. A related generalisation was given by Shevelev [32] where instead the cycle lengths belong to S.

Riordan [28] also gave the recurrence equation

$$\begin{aligned} K_{3,n}=n^2 K_{3,n-1}+n(n-1)K_{3,n-2}+2n(n-1)(n-2)K_{3,n-3}+s_n \end{aligned}$$

(3)

where

$$\begin{aligned} s_n=-n s_{n-1}-(n-1)2^n \end{aligned}$$

(4)

and $s_0=1$. We find that using (3)–(4) is the fastest way to compute $K_{3,n}$ in practice. This is unsurprising as it requires only a finite number of arithmetic operations to compute $K_{k,n}$ (from previous terms and the auxiliary function $s_n$). With the exception of Kerawala’s recurrence (Theorem 1), all other formulae for $L_{k,n}$ we survey require a number of arithmetic operations which grows with n.

Dulmage [9] posed an “explicit (though complicated)” formula for $L_{3,n}$ as a problem, which was later refined by Dulmage and McMaster [10] who gave

$$\begin{aligned} K_{3,n}=\sum _{j=0}^{\left\lfloor n/2 \right\rfloor } \sum _{i=0}^{n-2j} (-1)^{j+i} \left( {\begin{array}{c}n\\ j,j,i\end{array}}\right) a_{j+i,i}\,a^2_{n-j-i,n-2j-i}, \end{aligned}$$

where

$$\begin{aligned} a_{n,p}=\sum _{i=0}^p (-1)^i \left( {\begin{array}{c}p\\ i\end{array}}\right) (n-i)!, \end{aligned}$$

which they used to compute $K_{3,n}$ for $n \le 40$ on an IBM 370/158 in approximately five seconds.

Bogart and Longyear [3] gave

$$\begin{aligned} L_{3,n}=n!^2 \sum _{p+q+r+s+d=n} \frac{(-1)^{p+q+r}\,2^s\,(p+d)!\,(q+d)!\,(r+d)!}{d!^2\,p!\,q!\,r!\,s!} \end{aligned}$$

where the sum is over all non-negative integers p, q, r, s and d which sum to n. Bogart and Longyear used their formula to find $L_{3,n}$ for $n \le 11$ (although some typographical errors were pointed out in [39]) and gave an approximation of $L_{3,n}$ for $n \le 20$.

Riordan [28] said Yamamoto [44] found the equation

$$\begin{aligned} R_{3,n}= \sum _{a+b+c=n} n (n-3)! (-1)^b \frac{2^c a!}{c!} \left( {\begin{array}{c}3a+b+2\\ b\end{array}}\right) , \end{aligned}$$

where the sum is over all non-negative integers a, b and c which sum to n.

Goulden and Jackson [18] showed that $K_{3,n}$ is the coefficient of $x^n/n!$ in the series expansion of

$$\begin{aligned} e^{2x} \sum _{i \ge 0} \frac{x^i}{(1+x)^{3i+3}}. \end{aligned}$$

This was generalised by Gessel [14] to count pairs of discordant derangements. Specifically, the number of pairs $(\pi ,\sigma )$ of derangements of [n] such that $\pi \sigma ^{-1}$ is also a derangement, $\pi$ has a cycles and $\sigma$ has b cycles is the coefficient of $\alpha ^a \beta ^b x^n$ in

$$\begin{aligned} e^{2\alpha \beta x} \sum _{i \ge 0} \frac{(\alpha )_i (\beta )_i}{i!} \frac{x^i}{(1+\alpha x)^{i+\beta }(1+\beta x)^{i+\alpha }(1+x)^{n+\alpha \beta }} \end{aligned}$$

where $(\alpha )_i=\alpha (\alpha +1)\cdots (\alpha +i-1)$.

Riordan [28] gave the recurrence congruence $R_{3,n+p} \equiv 2 R_{3,n} \pmod p$ for all odd primes p, provided $n \ge 3$, which was generalised by Carlitz [5] to $R_{3,n+t} \equiv 2^t R_{3,n} \pmod t$ for all $t \ge 1$. This was generalised to k-line Latin rectangles in [39], who gave

$$\begin{aligned} R_{k,n+t} \equiv \big ((-1)^{k-1}(k-1)!\big ){}^t\,R_{k,n} \pmod t \end{aligned}$$

(5)

for all $k \ge 1$ and $t \ge 1$ provided $n \ge k$. It was also noted in [39] that some primes p do not divide $R_{3,n}$ for any $n \ge 3$. For example, any prime

$$\begin{aligned} p \in \{3,5,11,29,37,41,43,53,67,79,83,97\} \end{aligned}$$

does not divide $R_{3,n}$ for all $n \ge 3$ (see [36, 39]).

There are also published formulae for four-line Latin rectangles by Light Jr. [21], Pranesachar et al. [2, 25] and Doyle [6]. We will discuss Doyle’s formula in the next section, which is, by far, the best method for finding the exact value of $R_{k,n}$ for $4 \le k \le 6$.

3 Doyle’s Formula

The computation of $L_{k,n}$ for $k \in \{1,2,3\}$ can be considered as effectively “solved”: we have $L_{1,n}=n!$, and $L_{2,n}$ can be found by computing the number of derangements $d_n$. For $n=3$, using Riordan’s recurrence (i.e., (3)–(4)), we computed $L_{3,n}$ for all $n \le 2^{21} \approx 2 \times 10^6$ on a desktop computer in under 38 hours (we use these results in Section 5). As we will see, computing $L_{k,n}$ for $k \ge 4$ is much more difficult practically.

Exact enumeration is difficult for $k>3$. – Skau [33]

In [6], Doyle gave formulae for $K_{k,n}$ for $k \le 4$ and indicated how they could be generalised to arbitrary k. This generalisation was subsequently used in [36] to find values $K_{k,n}$ when $k \le 6$ (although a proof was not given). For the sake of rigour, we will give a proof of Doyle’s generalised formula. The overall idea of the formula is to use Inclusion-Exclusion on ordered n-tuples of columns. We include all such n-tuples, then exclude those that clash. To count the number of n-tuples that clash, Möbius Inversion is used on each column, counting the number of arrangements that “avoid” a certain substructure (e.g. symbol 1 in row 2).

Let ${\mathcal {R}}$ be the set of non-negative integer vectors $\mathbf {s}=(s_i)_{1 \le i \le 2^{k-1}}$ such that $\sum _i s_i=n$. For $1 \le i \le 2^{k-1}$, let $\varDelta _i=(\delta _{ij})_{1 \le j \le 2^{k-1}}$, where $\delta _{ij}$ is the Kronecker $\delta$-function. For any non-negative integer i let $b_j(i)$ be the j-th binary digit of i, for example $\big (b_j(13)\big )_{j \ge 1}=(1,0,1,1,0,0,\ldots )$. Let $||\mathbf {s}||=\sum _{1 \le i \le 2^{k-1}} \sum _{1 \le j \le k-1} s_i b_j(i)$.

Theorem 2

$$\begin{aligned} K_{k,n}=\sum _{\mathbf {s} \in {\mathcal {R}}} (-1)^{||\mathbf {s}||} \left( {\begin{array}{c}n\\ {s_1,s_2,\ldots ,s_{2^{k-1}}}\end{array}}\right) \prod _{i=1}^{2^{k-1}} g\big (\mathbf {s}-\varDelta _i\big )^{s_i} \end{aligned}$$

(6)

where for $\mathbf {a}=(a_1,a_2,\ldots ,a_{2^{k-1}})$,

$$\begin{aligned} g(\mathbf {a})=\sum _{P \in {\mathcal {P}}_{k-1}} \prod _{p \in P} (-1)^{|p|-1}(|p|-1)! f_p(\mathbf {a}) \end{aligned}$$

(7)

where ${\mathcal {P}}_{k-1}$ is the set of partitions of $\{1,2,\ldots ,k-1\}$ and

$$\begin{aligned} f_p(\mathbf {a})=\sum _{i:b_j(i)=0 \forall j \in p} a_i \end{aligned}$$

(8)

for all $p \subseteq \{1,2,\ldots ,k-1\}$.

Proof

Let ${\mathcal {S}}$ be the set of $n \times (k-1)$ $\,(0,1)$-matrices and let $D=(d_{ij}) \in {\mathcal {S}}$. We say a matrix $D'=(d'_{ij}) \in {\mathcal {S}}$ avoids D if $d'_{ij}=0$ whenever $d_{ij}=1$. Let

• $Z=Z(D)=\{D' \in {\mathcal {S}}:D'$ avoids $D\}$,

• $X=X(D)=\{D' \in Z:$ every column of $D'$ contains a unique 1 and every row of $D'$ contains at most one $1\}$ and

• $Y_P=Y_P(D)=\{D' \in Z:$ every column of $D'$ contains a unique 1 and columns c and $c'$ of $D'$ are identical if there exists a part $p \in P$ for which both $c,c' \in p\}$ for any partition $P \in {\mathcal {P}}_{k-1}$.

Let $L=(l_{ij})$ be a normalised $k \times n$ Latin rectangle. From L we can construct a $k \times n \times n$ $\,(0,1)$-array $P=(p_{ijr})$ where $p_{ijr}=1$ whenever $l_{ij}=r$. We will find it helpful to think of P as an ordered n-tuple of $n \times (k-1)$ $\,(0,1)$-matrices $(M_1,M_2,\ldots ,M_n) \in X^n$ where the (r, i)-th coordinate of $M_j$ is $p_{(i+1)jr}$. An example of this equivalence is given in Fig. 1; the shaded entries must be zero since L is normalised.

For $1 \le j \le n$, and $D \in {\mathcal {S}}$, let $D_j$ be the matrix formed from D after the j-th row has been converted to a row of ones. By Inclusion-Exclusion

$$\begin{aligned} K_{k,n}=\sum _{D \in {\mathcal {S}}} (-1)^{\sigma (D)} |X(D_1)|\cdot |X(D_2)|\cdots |X(D_n)| \end{aligned}$$

(9)

where $\sigma (D)$ is the number of ones in D. The Inclusion-Exclusion works as follows. We include all the n-tuples $(M_1,M_2,\ldots ,M_n) \in X^n$ for which the j-th row of each $M_j$ does not contain a 1. Suppose $D \in {\mathcal {S}}$ where the (r, i)-th element of D is 1. Suppose also that we attempt to construct a normalised Latin rectangle L from $(M_1,M_2,\ldots ,M_n) \in X^n$ in which every $M_j$ avoids D. Then this would imply that symbol r does not appear in row $i+1$ of L, giving rise to a clash. Therefore, we then exclude all such n-tuples for which there exists some $D \in {\mathcal {S}}$ containing a 1 for which every $M_j$ avoids D.

From any $D \in {\mathcal {S}}$ we construct a vector $\mathbf {s}=(s_1,s_2,\ldots ,s_{2^{k-1}})$ in the following way. Each row of D can be considered the binary digits of some number between 0 and $2^{k-1}-1$ (inclusive). Let $s_t$ be the number of rows in D that are the binary digits of t. To be consistent with [36], we use $s_{2^{k-1}}=s_0$ in the statement of the theorem. This does not affect the proof since the first $k-1$ binary digits of both 0 and $2^{k-1}$ are zero. Hence $\sigma (D)=||\mathbf {s}||$. This process is depicted in an example in Fig. 2.

The number of $D \in {\mathcal {S}}$ that give rise to a given $\mathbf {s} \in {\mathcal {R}}$ is given by the multinomial coefficient

$$\begin{aligned} \left( {\begin{array}{c}n\\ s_1,s_2,\ldots ,s_{2^{k-1}}\end{array}}\right) . \end{aligned}$$

We can partition ${\mathcal {S}}$ according to $\mathbf {s} \in {\mathcal {R}}$, whence two matrices D and $D'$ in the same part have $|X(D)|=|X(D')|$. For each $\mathbf {s} \in {\mathcal {R}}$, choose a representative $D^*=D^*(\mathbf {s})$. Hence (9) becomes

$$\begin{aligned} K_{k,n}=\sum _{\mathbf {s} \in {\mathcal {R}}} (-1)^{||\mathbf {s}||} \left( {\begin{array}{c}n\\ s_1,s_2,\ldots ,s_{2^{k-1}}\end{array}}\right) |X(D^*_1)| \cdot |X(D^*_2)| \cdots |X(D^*_n)|. \end{aligned}$$

(10)

If we define $g(\mathbf {a})=|X(D^*(\mathbf {a}))|$, then (10) becomes

$$\begin{aligned} K_{k,n}=\sum _{\mathbf {s} \in {\mathcal {R}}} (-1)^{||\mathbf {s}||} \left( {\begin{array}{c}n\\ s_1,s_2,\ldots ,s_{2^{k-1}}\end{array}}\right) \prod _{i=1}^{2^{k-1}} g\big (\mathbf {s}-\varDelta _i+\varDelta _{2^{k-1}-1}\big )^{s_i} \end{aligned}$$

(11)

since if two rows i and $i'$ of $D^*$ are the same, then $|X(D^*_i)|=|X(D^*_{i'})|$. The first $k-1$ binary digits of $2^{k-1}-1$ are all ones, so we will later delete $\varDelta _{2^{k-1}-1}$ from (11) since it makes no difference to (7) and (8). It is now sufficient to prove (7).

We can decompose any $N \in Y_P$ into |P| submatrices $N_p$ for each $p \in P$. An example of this decomposition is given in Fig. 3; we shade the entries that N and the $N_p$ avoid. For a given $D \in {\mathcal {S}}$, the number of possible submatrices $N_p$ is equal to the number of rows of D that map to a number t for which the j-th binary digit of t is zero for all $j \in p$ (informally, we can choose a single row of ones in $N_p$ that avoids the ones in D); this number is given by (8). Therefore $|Y_P|=\prod _{p \in P} f_p(\mathbf {a})$. It is now sufficient to show that $g(\mathbf {a})=\sum _{P \in {\mathcal {P}}_{k-1}} |Y_P| \prod _{p \in P} (-1)^{|p|-1}(|p|-1)!$.

If $P,Q \in {\mathcal {P}}_{k-1}$ and P is a refinement of Q, then we write $P \unlhd Q$ and $Q \unrhd P$. If additionally $P \ne Q$, then we write $P \lhd Q$ and $Q \rhd P$. For the remainder of this proof, $P \in {\mathcal {P}}_{k-1}$ will be the partition of cardinality $k-1$, i.e., we will henceforth assume $P=\{\{1\},\{2\},\ldots ,\{k-1\}\}$.

By Inclusion-Exclusion,

$$\begin{aligned} |X|=|Y_P| - \Big |\bigcup _{Q \rhd P} Y_Q\Big |=\sum _{Q \unrhd P} \mu (P,Q) |Y_Q| \end{aligned}$$

for some integer coefficients $\mu (P,Q)$, where $\mu (P,P)=1$. To ensure the matrices $M \in Y_P$ with some duplicated columns are counted 0 times overall, the coefficients $\mu (P,Q)$ must satisfy $\sum _{P \unlhd R \unlhd Q} \mu (P,R)=0$, whenever $P \lhd Q$. Therefore $\mu (P,Q)=-\sum _{P \unlhd R \lhd Q} \mu (P,R)$. Hence $\mu$ is the Möbius Function for the lattice of partitions ${\mathcal {P}}_{k-1}$ (see [30, pp. 359–360] or [43, pp. 333–336] for example). Since P is the partition of cardinality $k-1$, we find $\mu (P,Q)=(-1)^{k-1-|Q|} \prod _{q \in Q}(|q|-1)!=\prod _{q \in Q}(-1)^{|q|-1}(|q|-1)!$ $\square$

Fig. 1

Converting between L and $M_1,M_2,\ldots ,M_n$

Fig. 2

Finding $\mathbf {s}$ from D

Fig. 3

An example of an $N \in Y_P(D)$ where $P=\{\{1,2\},\{3\}\}$

Doyle’s equation for $K_{k,n}$ might seem intimidating at first, but if we fix a value of k, then $g(\mathbf {a})$ and the $f_p(\mathbf {a})$ are fixed multivariate polynomials. Afterwards, for each vector $\mathbf {s} \in {\mathcal {R}}$, computing its contribution to (6) is a straightforward task.

Doyle’s formula was used in [36] to find $R_{4,n}$ for $n \le 80$, $\,R_{5,n}$ for $n \le 28$ and $R_{6,n}$ for $n \le 13$. Motivation for finding more of these numbers comes from the following question in [35], which was posed after considering the data for $R_{4,n}$ and $R_{5,n}$ in [36].

Question 3

How do the prime power divisors $p^a$ of $R_{k,n}$ behave asymptotically for a fixed $k>p$ as $n \rightarrow \infty$ or as both $k \rightarrow \infty$ and $n \rightarrow \infty$?

Some general results concerning the divisors of $R_{k,n}$ were given by [39]. For example, we know that $p^{\left\lfloor n/p \right\rfloor }$ divides $R_{k,n+d}$ and $K_{k,n}$ when p is a prime and $d\ge k>p$.

4 Implementation

The first author gave a basic implementation of Doyle’s formula in C available from

$$\begin{aligned} \text {{http://code.google.com/p/latinrectangles/downloads/list}} \end{aligned}$$

(12)

which was used in [36] to find $R_{4,n}$, $R_{5,n}$ and $R_{6,n}$ for $n \le 80$, $n \le 28$ and $n \le 13$, respectively. Some basic, but significant improvements to the code were made, most notably:

• 1. Parallelisation.   In order to implement Doyle’s formula in parallel, we partition ${\mathcal {R}}$ into x parts, where x is the number of processes to be run in parallel. Using lexicographic ordering on ${\mathcal {R}}$, we find $x-1$ vectors $\mathbf {z}_1,\mathbf {z}_2,\ldots ,\mathbf {z}_{x-1}$, such that the x sets

  $$\begin{aligned}&\{\mathbf {s} \in {\mathcal {R}}:\mathbf {s}<\mathbf {z}_1\},\\&\{\mathbf {s} \in {\mathcal {R}}:\mathbf {z}_i \le \mathbf {s}<\mathbf {z}_{i+1}\}&\text {for } 1 \le i \le x-2,\\&\{\mathbf {s} \in {\mathcal {R}}:\mathbf {z}_{x-1} \le \mathbf {s}\} \end{aligned}$$

  have roughly the same cardinality. The Maple code used to find the vectors $\mathbf {z}_i$ is also available from (12).

• 2. Symmetry. We identify a group G that acts on ${\mathcal {R}}$ such that the contribution to (6) is invariant. For $k \in \{4,5,6\}$, the order of this group is $|G|=(k-1)!$. If $G(\mathbf {s})$ denotes the orbit of $\mathbf {s} \in {\mathcal {R}}$, then we only include the contribution of $\mathbf {s}$ to (6) whenever $\mathbf {s}=\min G(\mathbf {s})$, but multiplied by $|G(\mathbf {s})|$. While useful, exploiting this symmetry does not reduce run-time by a factor of |G| since (a) there is additional overhead and (b) sometimes $|G(\mathbf {s})|<|G|$.

Using the aforementioned improvements, we were able to find $R_{k,n}$ when

• $k=4$ and $n \le 150$, as at Sloane’s [34] A000573,

• $k=5$ and $n \le 40$, as in Fig. 4 (Fig. 5 gives the prime factorisations),

• $k=6$ and $n \le 15$, as in Fig. 6.

The improved version of the code can also be downloaded from (12).

Theoretically, if we found enough values of $R_{4,n}$, we could use Sister Celine’s technique [45] to propose a recurrence for $R_{4,n}$. However, it seems that very many values of $R_{4,n}$ would be required to find the recurrence, well beyond the values for $n \le 150$ we currently have.

Fig. 4

$R_{5,n}$ for small n

Fig. 5

Factorisation of $R_{5,n}$ for small n

Fig. 6

$R_{6,n}$ for small n

We can use (5) to check that the values are indeed correct. It is highly unlikely that an erroneous implementation would still return numbers that satisfy all of the recurrences implied by (5). We can also check that the number of digits is roughly what we expect based on the asymptotic formula and neighbouring data points.

Doyle’s formula iterates through each $\mathbf {s} \in {\mathcal {R}}={\mathcal {R}}(k,n)$; a “stars and bars” argument gives

$$\begin{aligned} |{\mathcal {R}}(k,n)|=\left( {\begin{array}{c}{n+2^{k-1}-1}\\ {2^{k-1}-1}\end{array}}\right) . \end{aligned}$$

So the number of iterations required to compute $R_{4,150}$, $R_{5,40}$ and $R_{6,15}$ is as given below:

$$\begin{aligned} \begin{array}{rcccr} |{\mathcal {R}}(4,150)| &{} = &{} \left( {\begin{array}{c}157\\ 7\end{array}}\right) &{} = &{} 407,340,975,756, \\ |{\mathcal {R}}(5,40)| &{} = &{} \left( {\begin{array}{c}55\\ 15\end{array}}\right) &{} = &{} 11,899,700,525,790, \\ |{\mathcal {R}}(6,15)| &{} = &{} \left( {\begin{array}{c}46\\ 31\end{array}}\right) &{} = &{} 511,738,760,544,\\ \end{array} \end{aligned}$$

although we do reduce the number of required iterations by a factor of $(k-1)!$ through symmetry. For the purpose of inspecting the divisors of $R_{k,n}$, we are particularly interested in finding more values of $R_{5,n}$.

We can estimate how many times longer it would take to compute $R_{k,n+1}$ than $R_{k,n}$ with the ratio

$$\begin{aligned} r_{k,n}:=\frac{|{\mathcal {R}}(k,n+1)|}{|{\mathcal {R}}(k,n)|}=\frac{n+2^{k-1}}{n+1}, \end{aligned}$$

which is close to 1 when n is much larger than $2^{k-1}$. For example, $r_{4,150} \approx 1.1$, $\,r_{5,40} \approx 1.4$ and $r_{6,15} \approx 2.9$. So it would be possible to find more values of $R_{4,n}$, but we had to stop somewhere.

Fig. 7

$\omega _2(R_{k,n})$ for $k \in \{3,4,5,6\}$ and small n a $\omega _2(R_{3,n})$ vs. n. b $\omega _2(R_{4,n})$ vs. n. c $\omega _2(R_{5,n})$ vs. n. d $\omega _2(R_{6,n})$ vs. n

Fig. 8

$\omega _3(R_{k,n})$ for $k \in \{4,5\}$ and small n. a $\omega _3(R_{4,n})$ vs. n. b $\omega _3(R_{5,n})$ vs. n

Fig. 9

Approximations of the gradients in Figs. 7 and 8

5 Divisors

The question of divisors of $R_n$, and in particular powers of 2 and 3, was originally raised by Alter [1] about Latin squares. This question, along with analogous questions for Latin rectangles, has been a hot research topic for Stones [35], and is linked to the autotopisms and automorphisms of Latin rectangles [4, 39] (see also [12, 38]), and orthomorphisms and partial orthomorphisms of finite cyclic groups [40, 41]. Divisors for the number of even/odd Latin squares have been used in proving special cases of the Alon-Tarsi Conjecture [7] (see also [8, 16, 37, 42]).

For any prime p, let $\omega _p(n)$ be the largest non-negative integer such that $p^{\omega _p(n)}$ divides n, that is the p-adic valuation of n. In Fig. 7, we plot $\omega _2(R_{k,n})$ for $k \in \{3,4,5,6\}$ and small n. In Fig. 8, we plot $\omega _3(R_{k,n})$ for $k \in \{4,5\}$ and small n. This data leads us to the following question.

Question 4

For any fixed k and fixed prime $p<k$, is

$$\begin{aligned} \omega _p(R_{k,n}) \ge \frac{\lfloor (k-1)/p \rfloor }{p-1} n-o(n)? \end{aligned}$$

The gradients implied by Question 4 match those in Fig. 9, which are suggested by the empirical data.

From [39], we already know that $\omega _p(R_{k,n+k}) \ge \left\lfloor n/p \right\rfloor$ for all primes $p<k$. Furthermore, it was shown in [4] that $\omega _p(R_n) \ge n/(p-1)-O(\log ^2 n)$ as $n \rightarrow \infty$. The divisors of $R_{k,n}$ for primes $p \ge k$ do not display the same characteristics as when $p<k$. For example, 3 does not divide $R_{3,n}$ for any $n \ge 3$, whereas a growing power of 3 divides $R_{4,n}$ [39].

In the case of three-line Latin rectangles, Riordan’s recurrence (3) gives us the ability to quickly compute values of $K_{3,n}$. We will now list some conjectures for three-line Latin rectangles that are motivated by these numbers.

Conjecture 5

• 1. $\omega _2(K_{3,n})=\omega _2(n!)=n-1$ whenever $n=2^i$ and $i \ge 2$. (True for $i \le 21$.)

• 2. $\omega _2(K_{3,n}) \ge \omega _2(n!)$ for all $n \ge 3$. (True for $n \le 2^{18}$.)

• 3. $\omega _2(K_{3,n})=n-i+1$ whenever $n=2^i-3$ and $i \ge 3$. (True for $i \le 21$.)

Kerawala’s recurrence (Theorem 1) suggests that $\omega _2(K_{3,n})$ might be less than usual when n has the form $2^i-3$, which is what we see in the numerical data.

Conjecture 5 implies that $\omega _2(K_{3,n})$ is not bounded below by $n-c$ for any constant c. Furthermore, Conjecture 5 implies that the error term for the first order approximation in Question 4 is at least $O(\log n)$.