Number of distinct tilings of an n X n square with free n-polyominoes

The newest "nice" OEIS sequence, A328020, was just published a few minutes ago.

Number of distinct tilings of an n X n square with free n-polyominoes.

This sequence counts tilings up to symmetries of the square. The sequence has six terms, but I'd like to see if folks here can extend it further.

Example

For n=4 there are 22 such grids, as shown in this image from the OEIS. Credit: Jeff Bowermaster, Illustration of A328020(4).

Challenge

Like this previous challenge, the goal of this challenge is to compute as many terms as possible in this sequence, which begins 1, 1, 2, 22, 515, 56734 and where the n-th term is the number of tilings of the n X n grid with n-polyominoes.

Run your code for as long as you'd like. The winner of this challenge will be the user who posts the most terms of the sequence, along with their code to generate it. If two users post the same number of terms, then whoever posts their last term earliest wins.

• So this is modulo symmetries of the square? Oct 8 '19 at 8:59
• @PeterTaylor, that's right. I've disambiguated this in the question. Oct 8 '19 at 17:34
• Naively I would say that the n'th entry would take number_of_fixed_n_polyominoes^(n-1) operations to calculate. So for n = 7, that would take 760^6 ≈ 2^57.4 operations. You can probably cut that down a lot, but it's a big number to start with... Oct 8 '19 at 17:49
• @Sliepen, I expect that you can cut down that by quite a lot just by backtracking. In particular, there are a lot of fixed polynomials that can't be placed in the corner, and once a valid polyomino is placed somewhere, it hugely constrains what can be placed next to it. Oct 8 '19 at 17:58
• @PeterKagey, you're right. I guess it helps if, given that you already have placed m n-polyominoes, you choose the next position to try to place a polyomino in the worst possible position, that you can cut it down a lot. Oct 8 '19 at 18:22

C, 7 terms

The seventh term is 19846102. (The first six are 1, 1, 2, 22, 515, 56734, as stated in the question).

#include <stdio.h>
#include <string.h>
#include <stdint.h>

#define N 7
#define ctz __builtin_ctzl

typedef uint64_t u64;

static u64 board[N*N] = { 0 };
static u64 adjacency_matrix[N*N] = { 0 };
static u64 count = 0;

static u64 check_symmetry()
{
static const u64 symmetries = {
{ 0,     +N, +1 },
{ N-1,   -1, +N },
{ N-1,   +N, -1 },
{ N*N-1, -1, -N },
{ N*N-1, -N, -1 },
{ N*N-N, +1, -N },
{ N*N-N, -N, +1 },
};

int order[N];

for (u64 i = 0; i < 7; ++i) {
u64 start = symmetries[i];
u64 dcol = symmetries[i];
u64 drow = symmetries[i];
memset(order, 0xFF, N*sizeof(int));

for (u64 row = 0, col = 0; col < N || (col = 0, ++row < N); ++col) {
u64 base = board[col + N*row];
u64 symmetry = board[start + dcol*col + drow*row];
u64 lex = 0;

while (order[lex] != symmetry && order[lex] != -1)
++lex;
order[lex] = symmetry;

if (lex < base)
return 0;

if (base < lex)
break;
}
}

return 1;
}

static void recurse(u64 mino, u64 cell, u64 occupied, u64 adjacent, u64 forbidden)
{
if (cell >= N) {
++mino;

if (mino == N) {
count += check_symmetry();
return;
}

u64 next = ctz(~occupied);
board[next] = mino;
recurse(mino, 1, occupied | 1ul << next, adjacency_matrix[next], 0);
return;
}

forbidden |= 1ul << next;
board[next] = mino;
recurse(mino, cell + 1, occupied | 1ul << next, adjacent | adjacency_matrix[next], forbidden);
}
}

int main(void)
{
for (u64 i = 0; i < N*N; ++i) {
if (i % N)
adjacency_matrix[i] |= 1ul << (i - 1);
if (i / N)
adjacency_matrix[i] |= 1ul << (i - N);
if (i % N != N - 1)
adjacency_matrix[i] |= 1ul << (i + 1);
if (i / N != N - 1)
adjacency_matrix[i] |= 1ul << (i + N);
}

recurse(0, 2, 3, 4 | 3 << N, 0);
printf("%ld\n", count);
}

Try it online! (for N=6, since N=7 would time out.)

On my machine, N=6 took 0.171s, and N=7 took 2m23s. N=8 would take a few weeks.

• This is wonderful! Let me know if you'd like to add it to the OEIS—which might result in you doxxing yourself—or if you'd like me to add it. Oct 9 '19 at 20:46
• @PeterKagey Please feel free to add it (: Oct 10 '19 at 8:28
• Fascinating check_symmetry function. Could you please give a brief explanation as I’m not familiar with the approach? Oct 11 '19 at 0:34
• @JohnRees It simply tests that the current board is lexicographically ≤ to all of its symmetries. Thus in any set of symmetrical boards, exactly one is counted: the lexicographical minimum. Oct 11 '19 at 6:34
• To do better than enumerating the solutions one by one, some kind of meet-in-the-middle is required. The problem is that there doesn't seem to be any way of splitting things up which gets significant clustering. E.g. using the same canonical ordering as this answer, placing 3 hexominoes I get an average of about 3.7 sets of hexominos per mask. I conclude that this approach is not likely to be beaten. Oct 19 '19 at 20:11

An extension to @Grimy's code gets N=8

This just underlines that @Grimy deserves the bounty:

I could prune the search tree by extending the code to check, after each finished polyomino, that the remaining free space is not partitioned into components of size not divisible by N.

On a machine where the original code took 2m11s for N=7, this takes 1m4s, and N=8 was computed in 33h46m. The result is 23437350133.

Here is my addition as a diff:

--- tilepoly.c  2019-10-11 12:37:49.676351878 +0200
+++ tilepolyprune.c     2019-10-13 04:28:30.518736188 +0200
@@ -51,6 +51,30 @@
return 1;
}

+static int check_component_sizes(u64 occupied, u64 total){
+    u64 queue[N*N];
+    while (total<N*N){
+        u64 count = 1;
+        u64 start = ctz(~occupied);
+        queue = start;
+        occupied |= 1ul << start;
+        for(u64 current=0; current<count; ++current){
+                free_adjacent &= ~(1ul << next);
+                queue[count++] = next;
+            }
+        }
+        if (count % N){
+            return 0;
+        }
+        total += count;
+    }
+    return 1;
+}
+
static void recurse(u64 mino, u64 cell, u64 occupied, u64 adjacent, u64 forbidden)
{
if (cell >= N) {
@@ -61,6 +85,9 @@
return;
}

+        if(!check_component_sizes(occupied,N*mino))
+            return;
+
u64 next = ctz(~occupied);
board[next] = mino;
recurse(mino, 1, occupied | 1ul << next, adjacency_matrix[next], 0);

Try it online!

• This is very nice.
– user9207
Oct 13 '19 at 15:18
• All we need now is a multithreaded simd version :)
– user9207
Oct 13 '19 at 15:50
• Oh that’s really cool! I actually considered this optimization, but didn’t think it’d be enough to reach N=8 in a reasonable time, so I didn’t bother implementing it. Oct 14 '19 at 21:54