a too-short sequence related to Latin squares

[hv at crypt.org]

NJAS: N. J. A. Sloane <njas at research.att.com>
HVDS: Hugo van der Sanden <hv at crypt.org>
HP: Hugo Pfoertner <all at abouthugo.de>

NJAS> How many n X n squares are there, filled with
NJAS> the numbers 1...n, such that in row or column k,
NJAS> for all k = 1...n, the number k appears at least once?
NJAS> 
NJAS> For n = 1 there is just 1,
NJAS> 
NJAS> for n = 2 there seem to be 6

HVDS> For n = 3 I make it [...] 2683

HP> I wrote a little program, which confirms a(3)=2683 and gets
HP> a(4)=223041730 after checking the 4^16 possible matrices. Getting the
HP> next term by brute force would require 5^25~=2.98*10^17 checks, which is
HP> a good reason to defer the computation until better ideas become
HP> available.

I can confirm a(4) = 223041730.

Concentrating first on the main diagonal, we can consider a canonically
'correct' example where we put _i_ into the (i, i) square; we can then
count the number of incorrect entries as k from 0 to n, and get:

  a(n) = sum_k=0^n { C(n, k) (n-1)^k n^((n-k)(n-k-1)) b_k(n) }

This breaks down as:
  C(n, k) is the number of ways we can choose which diagonal entries
are wrong;

  (n-1)^k is the number of ways we can fill in those wrong entries;

  n^((n-k)(n-k-1)) is the number of ways we can fill in the non-diagonal
entries in the rows and columns that are satisfied by a correct diagonal
entry;

  b_k(n) is an auxiliary function that counts the number of ways the
rows/columns with incorrect diagonal entries can be satisfied.

b_k(n) looks at the entries in (k^2 - k) primary squares (which are on both
an unsatisfied row and an unsatisfied column) and 2k(n-k) secondary squares
(which have either their row or column satisfied).

Then b_0(n) = 1; b_1(n) = (n^(n-1) - (n-1)^(n-1))^2; and generally:

  b_k(n) = sum_i=0^2k { c_k_i(n) (n^(n-k))^i (n^(n-k) - (n-1)^(n-k))^(2k-i) }

.. for functions c_k_i(n) each of which is a fixed polynomial in n.
Let me illustrate with a picture for k = 3:

  * B C . . . .
  D * F . . . .
  G H * . . . .
  . . .
  . . .
  . . .

The letters represent the (k^2 - k) primary squares, the asterisks the
main diagonal, and the dots the 2k(n-k) secondary squares. c_3_i(n)
represents the number of ways that BCDFGH can be filled in such that
exactly _i_ of these 3 rows and 3 columns are satisfied.

Calculating the polynomials c_k_i(n) for some k involves iterating over
3 values for each of the primary squares: for each square we look at its
row number, its column number (by definition different) and a third value
to represent the (n-2) cases of 'anything else'; for the general case then,
this involves 3^(k^2 - k) iterations for each k; to calculate a(n) we'll
need all the c_k_i(n) for k < n, and c_n_{2n}(n) (note that for the
degenerate case k=n we get b_n(n) = c_n_{2n}(n)).

I don't know how to calculate an individual c_k_i(n) without covering
all values of i, but this approach reduces the number of iterations from
the brute force n^(n^2) to about 3^(n^2 - n).

Here's a short piece of perl code that finds the c_k_i(n), using the very
useful Algorithm::Loops module available on CPAN [1]:

  my $k = 3;   # for example
  my(@primary, @row, @col, @c);
  my $index = 0;
  for my $row (1 .. $k) {
    for my $col (1 .. $k) { 
      next if $row == $col;          # ignore main diagonal
      push @{ $row[$row] }, $index;  # record row indexes
      push @{ $col[$col] }, $index;  # record column indexes
      $primary[$index++] = [ 0, $row, $col ];
      # for each primary, the nested loops will iterate over these 3 values;
      # we use '0' as a marker for the (n-2) 'anything else' values
    }
  }

  use Algorithm::Loops 'NestedLoops';
  NestedLoops(\@primary, sub {
    my $power = grep !$_, @_; # count the zeroes
    my $match = 0;
    for my $v (1 .. $k) { 
      ++$match if grep $_[$_] == $v, @{ $row[$v] };
      ++$match if grep $_[$_] == $v, @{ $col[$v] };
    } 
    ++$c[$match][$power];
  });


It took me 25 seconds to calculate the polynomials for k=4, so I'd expect
k=5 to take about 46 hours, which is at least within the bounds of
reason. (And k=6 to take about 310 years, which isn't; I'm sure that
rewriting in C would offer a speedup between 10- and 100-fold, which
would reduce k=5 to 'easy' but still wouldn't bring k=6 within reach.)

The result of this code is an array of polynomials in (n-2); here are all
the non-zero results for k <= 4, writing 'x' for (n-2):

c_0_0 = 1

c_1_0 = 1

c_2_0 = x^2
c_2_1 = 4x
c_2_2 = 4

c_3_0 = x^6
c_3_1 = 12x^5 + 6x^4
c_3_2 = 54x^4 + 48x^3 + 9x^2
c_3_3 = 112x^3 + 126x^2 + 36x + 2
c_3_4 = 105x^2 + 120x + 30
c_3_5 = 36x + 30
c_3_6 = 2

c_4_0 = x^12
c_4_1 = 24x^11 + 24x^10 + 8x^9
c_4_2 = 240x^10 + 456x^9 + 348x^8 + 120x^7 + 16x^6
c_4_3 = 1296x^9 + 3480x^8 + 3936x^7 + 2328x^6 + 744x^5 + 120x^4 + 8x^3
c_4_4 = 4092x^8 + 13656x^7 + 19480x^6 + 15096x^5 + 6756x^4 + 1736x^3 + 252x^2 + 24x + 2
c_4_5 = 7632x^7 + 29256x^6 + 47880x^5 + 42696x^4 + 22024x^3 + 6432x^2 + 960x + 56
c_4_6 = 8056x^6 + 33384x^5 + 58272x^4 + 54376x^3 + 28344x^2 + 7728x + 844
c_4_7 = 4272x^5 + 18048x^4 + 31184x^3 + 27456x^2 + 12288x + 2232
c_4_8 = 828x^4 + 3312x^3 + 5100x^2 + 3576x + 962

Calculating on these values confirms the results found so far, and can
be used to discover that a(5) = 5788162687776805 + 1024(c_5_10(5));
I hope to have a number (and the c_5_i(n) polynomials) in 2 or 3 days.

I suspect there are further symmetries in the search space that could
bring a(6) within reach, but obviously a formula for these polynomials
would be preferable; no obvious approach (in either direction) springs
out at me though.

Another possibly interesting sequence is b_n(n) = c_n_2n(n) which
represents the same problem for squares missing the main diagonal;
starting at n=0, this goes:
  1, 0, 0, 2, 68258, ?

Hugo van der Sanden

[1] You can find the Algorithm::Loops module at:
  http://search.cpan.org/~tyemq/Algorithm-Loops-1.031/lib/Algorithm/Loops.pm