# Counting generalized polyominoes

This challenge will have you count pseudo-polyforms on the snub square tiling.

I think that this sequence does not yet exist on the OEIS, so this challenge exists to compute as many terms as possible for this sequence.

Update: this is now on the OEIS as A309159: Number of generalized polyforms on the snub square tiling with n cells.

# Definitions

The snub square tiling is a semiregular tiling of the plane that consists of equilateral triangles and squares.

A pseudo-polyform on the snub square tiling is a plane figure constructed by joining together these triangles and squares along their shared sides, analogous to a polyomino. Here is an example of a six-cell and an eight-cell pseudo-polyform:

# Examples

For n = 1 there are two 1-cell pseudo-polyforms, namely the square and the triangle:

For n = 2 there are two 2-cell pseudo-polyforms, namely a square with a triangle and two triangles.

For n = 3 there are four 3-cell pseudo-polyforms.

# Challenge

The goal of this challenge is to compute as many terms as possible in this sequence, which begins 2, 2, 4, ... and where the n-th term is the number of n-cell pseudo-polyforms up to rotation and reflection.

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. If two users post the same number of terms, then whoever posts their last term earliest wins.

(Once there are enough known terms to prove that this sequence does not already exist in the OEIS, I'll create an entry in the OEIS and list the contributor as a co-author if he or she desires.)

• In a code-challenge challenge, answers should be code, not terms of a sequence. Also, the number of terms one manages to discover before too much time passes in their opinion (or before the program overflows) isn't an objective winning criterion. Commented Jul 6, 2019 at 17:09
• @EriktheOutgolfer, last time I used a code-challenge for a similar-style problem, and it worked really well. Commented Jul 6, 2019 at 17:11
• Hm... I disagree on the part that you can run the program for as long as you like, since different people have different levels of patience (e.g. person 1 might let it run for 7 days, while person 2 might let it run for 30 days), and so it's subjective, although others may consider it to be a "good subjective". Just saying that it's not really objective. Regarding the code-challenge tag, I just wanted to make sure you're asking for code and not just a list of terms (Project Euler style). :P Commented Jul 6, 2019 at 17:19
• Anyone want to confirm or dispute 2, 2, 4, 10, 28, 79, 235, 720, 2254, 7146, 22927, 74137, 241461, 790838, 2603210, 8604861? Commented Jul 12, 2019 at 10:02
• @PeterTaylor I'm getting the same numbers Commented Jul 12, 2019 at 18:48

Now that not only the comments document that Peter Taylor was the first one to give enough terms to search on OEIS, I can give my results.

( 1 - 10) 2, 2, 4, 10, 28, 79, 235, 720, 2254, 7146,
(11 - 15) 22927, 74137, 241461, 790838, 2603210,
(16 - 18) 8604861, 28549166, 95027832,
(19 - 22) 317229779, 1061764660, 3562113987, 11976146355


Earlier, I counted hexagonal polyominoes. Except for some optimizations, what I'm doing here is very similar.

The elements of the tiling are represented like this: You can go in an almost straight line from left to right (in the first picture), alternating between squares and rectangles. There are almost parallel further lines, wiggling in opposite directions. Together, they miss some triangles. There are similar almost straight parallel lines from bottom to top, containing the missing triangles. Now ignore the wiggling and use a Cartesian coordinate system, but only use odd numbers for the coordinates of the squares. Then the triangles naturally get coordinate pairs with one even and one odd coordinate. Pairs with both coordinates even don't represent elements of the tiling.

(You could as well use even numbers for the coordinates of the squares. I guess I decided this way because I thought about reflection before rotation.)

Save the program in a file with a name like cgp.hs and compile with ghc -O2 -o cgp cgp.hs. It takes either one numeric command line argument and computes the number of polyominoes of that size, or none, in which case it computes values until stopped.

{-# LANGUAGE BangPatterns #-}

import Data.List(sort)
import qualified Data.Set as S
import System.Environment(getArgs)

data Point = P !Int !Int deriving (Eq,Ord)

start :: Point
start = P 1 1

redsq :: Point -> Bool
redsq (P x y) = (x+y) mod 4 == 2

neighs :: Point -> [Point]
neighs (P x y) =
case (even x, even y) of
(False,False) -> [P x (y+1), P (x+1) y, P x (y-1), P (x-1) y]
(True, False) -> (P x (c y (x+y+1))) : opt [P (x-1) y, P (x+1) y]
(False,True ) -> (P (c x (x+y-1)) y) : opt [P x (y-1), P x (y+1)]
where
opt = filter ok
ok p = p>start || not (redsq p)
c z m = if m mod 4 == 0 then z+2 else z-2

count :: S.Set Point -> S.Set Point -> [Point] -> Int -> Int -> Int -> Int -> Int
count use _    _            0 c r y =
if check (S.toAscList use) (y==r)
then c+1
else c
count _   _    []           _ c _ _ = c
count use seen (p:possible) n c r y =
let !c' = count use seen possible n c r y
new = filter (S.notMember seen) $neighs p !r' = if redsq p then r+1 else r !y' = if redsq (mirror p) then y+1 else y !n' = n-1 in if r'+n' < y' then c' else count (S.insert p use) (foldr S.insert seen new) (new++possible) n' c' r' y' class Geom g where translate :: Int -> Int -> g -> g rot :: g -> g mirror :: g -> g instance Geom Point where translate dx dy (P x y) = P (dx+x) (dy+y) rot (P x y) = P (2-y) x -- rotate around (1,1) mirror (P x y) = P x (-y) instance (Geom g, Ord g) => Geom [g] where translate x y = map$ translate x y
rot = sort . map rot
mirror = sort . map mirror

normalize :: [Point] -> [Point]
normalize pol = let (P x y) = head (filter redsq pol)
in translate (1-x) (1-y) pol

check :: [Point] -> Bool -> Bool
check pol !cm = let rotated = take 4 $iterate rot pol mirrored = if cm then map mirror rotated else [] alts = map normalize (tail rotated ++ mirrored) in all (pol<=) alts f :: Int -> Int f 0 = 1; f 1 = 2; f 2 = 2 f n = count S.empty S.empty [start] n 0 0 0 output :: Int -> IO () output n = putStrLn$ show n ++ ": " ++ show (f n)

main = do args <- getArgs
case args of
[]  -> mapM_ output [1..]



Try it online!

• Looks like you have a better tile representation than me. Would you mind explaining how it works? Commented Jul 13, 2019 at 20:02

## 2, 2, 4, 10, 28, 79, 235, 720, 2254, 7146, 22927, 74137, 241461, 790838, 2603210, 8604861, 28549166, 95027832

I'm going to put a stake in the ground before Christian Sievers posts an answer for n=18. This is as far as I can go with the current code and 16GB of RAM. I've already had to sacrifice some speed to reduce the memory usage, and I'm going to have to do so even more. I have some ideas...

This snippet is the SVG from the first comment.

<svg xmlns="http://www.w3.org/2000/svg" width="130" height="130">
<path style="stroke:none; fill:#f22" d="M 72,72 l -14.235,53.1259 -53.1259,-14.235 14.235,-53.1259 z" />  <!-- "Anticlockwise" square -->
<path style="stroke:none; fill:#44f" d="M 72,72 l 53.1259,-14.235 -14.235,-53.1259 -53.1259,14.235 z" />  <!-- "Clockwise" square -->

<path style="stroke:none; fill:#4f4" d="M 72,72 l 38.89,38.89 14.235,-53.1259 z" />  <!-- "NE" triangle -->
<path style="stroke:none; fill:#ff4" d="M 72,72 l 38.89,38.89 -53.1259,14.235 z" />  <!-- "SW" triangle -->
<path style="stroke:none; fill:#4ff" d="M 72,72 m -53.1259,-14.235 l 38.89,-38.89 -53.1259,-14.235 z" />  <!-- "NW" triangle -->

<path style="stroke:#000; fill:none" d="M 72,72 m 38.89,38.89 l 14.235,-53.1259 -14.235,-53.1259 -53.1259,14.235 -53.1259,-14.235 14.235,53.1259 -14.235,53.1259 53.1259,14.235 53.1259,-14.235" />
</svg>

Code is C#. I ran it with .Net Core 2.2.6 under Linux.

#define SUPERLIGHT
using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;

namespace Sandbox
{
// https://codegolf.stackexchange.com/questions/187763/counting-generalized-polyominoes
// Count polyominos on the snub square tiling.

// We index the tiles using the following basic element, which tiles like a square:
/*
<?xml version="1.0" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg xmlns="http://www.w3.org/2000/svg" width="130" height="130">
<path style="stroke:none; fill:#f22" d="M 72,72 l -14.235,53.1259 -53.1259,-14.235 14.235,-53.1259 z" />  <!-- "Anticlockwise" square -->
<path style="stroke:none; fill:#44f" d="M 72,72 l 53.1259,-14.235 -14.235,-53.1259 -53.1259,14.235 z" />  <!-- "Clockwise" square -->

<path style="stroke:none; fill:#4f4" d="M 72,72 l 38.89,38.89 14.235,-53.1259 z" />  <!-- "NE" triangle -->
<path style="stroke:none; fill:#ff4" d="M 72,72 l 38.89,38.89 -53.1259,14.235 z" />  <!-- "SW" triangle -->
<path style="stroke:none; fill:#4ff" d="M 72,72 m -53.1259,-14.235 l 38.89,-38.89 -53.1259,-14.235 z" />  <!-- "NW" triangle -->
<!-- There's a "SE" triangle, but it's unfilled -->

<path style="stroke:#000; fill:none" d="M 72,72 m 38.89,38.89 l 14.235,-53.1259 -14.235,-53.1259 -53.1259,14.235 -53.1259,-14.235 14.235,53.1259 -14.235,53.1259 53.1259,14.235 53.1259,-14.235" />
</svg>
*/
// In terms of symmetries, we have rotation by 90 degrees and reflection, possibly with glide.
// We obviously want a canonical representation.
//   Reflection interchanges "anticlockwise" and "clockwise" squares, so we shall require at least as many anticlockwise as clockwise.
//   Rotation anticlockwise by 90 maps NE -> NW -> SW -> SE -> NE. We rotate to get a standard necklace.
//   Further ties must be broken lexicographically, after translating to give minimum X and Y of 0.
class PPCG187763
{

internal static void Main()
{
SanityChecks();

var polyominos = new HashSet<TileSet>();
polyominos.Add(new TileSet(Enumerable.Repeat(new Tile { X = 0, Y = 0, Shape = TileShape.SE }, 1)));
polyominos.Add(new TileSet(Enumerable.Repeat(new Tile { X = 0, Y = 0, Shape = TileShape.Anticlockwise }, 1)));
Console.WriteLine($"1\t{polyominos.Count}"); for (int tileCount = 2; tileCount < 60; tileCount++) { var sw = new Stopwatch(); sw.Start(); var nextPolyominos = new HashSet<TileSet>(); // TODO This can be greatly optimised by tracking discarded insertion points foreach (var polyomino in polyominos) { foreach (var neighbour in polyomino.SelectMany(tile => tile.Neighbours).Distinct()) { if (!polyomino.Contains(neighbour)) nextPolyominos.Add(new TileSet(polyomino.Concat(Enumerable.Repeat(neighbour, 1)))); } } polyominos = nextPolyominos; Console.WriteLine($"{tileCount}\t{polyominos.Count}\t{sw.ElapsedMilliseconds}ms");
}
}

private static void SanityChecks()
{
var cluster = new HashSet<Tile>();
cluster.Add(new Tile { Shape = TileShape.Anticlockwise });
for (int i = 0; i < 3; i++)
{
foreach (var tile in cluster.SelectMany(tile => tile.Neighbours).ToList()) cluster.Add(tile);
}

foreach (var tile in cluster)
{
foreach (var neighbour in tile.Neighbours)
{
if (!neighbour.Neighbours.Contains(tile))
{
throw new Exception("Assertion failed: adjacency isn't symmetric");
}

if (!tile.Flip().Neighbours.Contains(neighbour.Flip()))
{
throw new Exception("Assertion failed: flip doesn't preserve adjacency");
}

if (!tile.Rot().Neighbours.Contains(neighbour.Rot()))
{
throw new Exception("Assertion failed: rot doesn't preserve adjacency");
}

if (!tile.Equals(tile.Rot().Rot().Rot().Rot()))
{
throw new Exception("Assertion failed: rot^4 should be identity");
}
}
}
}

struct Tile : IComparable<Tile>
{
public TileShape Shape { get; set; }
public sbyte X { get; set; }
public sbyte Y { get; set; }

public IEnumerable<Tile> Neighbours
{
get
{
switch (Shape)
{
case TileShape.Anticlockwise:
yield return new Tile { X = X, Y = Y, Shape = TileShape.SE };
yield return new Tile { X = X, Y = Y, Shape = TileShape.SW };
yield return new Tile { X = X, Y = (sbyte)(Y - 1), Shape = TileShape.NW };
yield return new Tile { X = (sbyte)(X - 1), Y = Y, Shape = TileShape.NE };
break;

case TileShape.Clockwise:
yield return new Tile { X = X, Y = Y, Shape = TileShape.SE };
yield return new Tile { X = X, Y = Y, Shape = TileShape.NE };
yield return new Tile { X = X, Y = (sbyte)(Y + 1), Shape = TileShape.SW };
yield return new Tile { X = (sbyte)(X + 1), Y = Y, Shape = TileShape.NW };
break;

case TileShape.NE:
yield return new Tile { X = X, Y = Y, Shape = TileShape.SW };
yield return new Tile { X = X, Y = Y, Shape = TileShape.Clockwise };
yield return new Tile { X = (sbyte)(X + 1), Y = Y, Shape = TileShape.Anticlockwise };
break;

case TileShape.NW:
yield return new Tile { X = X, Y = Y, Shape = TileShape.SE };
yield return new Tile { X = (sbyte)(X - 1), Y = Y, Shape = TileShape.Clockwise };
yield return new Tile { X = X, Y = (sbyte)(Y + 1), Shape = TileShape.Anticlockwise };
break;

case TileShape.SE:
yield return new Tile { X = X, Y = Y, Shape = TileShape.NW };
yield return new Tile { X = X, Y = Y, Shape = TileShape.Clockwise };
yield return new Tile { X = X, Y = Y, Shape = TileShape.Anticlockwise };
break;

case TileShape.SW:
yield return new Tile { X = X, Y = Y, Shape = TileShape.NE };
yield return new Tile { X = X, Y = (sbyte)(Y - 1), Shape = TileShape.Clockwise };
yield return new Tile { X = X, Y = Y, Shape = TileShape.Anticlockwise };
break;

default:
throw new NotSupportedException();
}
}
}

public Tile Flip()
{
// We'll flip vertically.
switch (Shape)
{
case TileShape.Anticlockwise:
return new Tile { Shape = TileShape.Clockwise, X = X, Y = (sbyte)-Y };
case TileShape.Clockwise:
return new Tile { Shape = TileShape.Anticlockwise, X = (sbyte)(X + 1), Y = (sbyte)-Y };
case TileShape.NE: // G
return new Tile { Shape = TileShape.SE, X = (sbyte)(X + 1), Y = (sbyte)-Y };
case TileShape.NW: // Cy
return new Tile { Shape = TileShape.SW, X = X, Y = (sbyte)-Y };
case TileShape.SE: // W
return new Tile { Shape = TileShape.NE, X = X, Y = (sbyte)-Y };
case TileShape.SW: // Y
return new Tile { Shape = TileShape.NW, X = (sbyte)(X + 1), Y = (sbyte)-Y };
default:
throw new NotSupportedException();
}
}

public Tile Rot()
{
// Anti-clockwise rotation: (x, y) = (-y, x)
// But there will be offsets to account for the positions within the cell
switch (Shape)
{
case TileShape.Anticlockwise:
return new Tile { Shape = TileShape.Anticlockwise, X = (sbyte)-Y, Y = X };
case TileShape.Clockwise:
return new Tile { Shape = TileShape.Clockwise, X = (sbyte)(-Y - 1), Y = X };
case TileShape.NE:
return new Tile { Shape = TileShape.NW, X = (sbyte)-Y, Y = X };
case TileShape.NW:
return new Tile { Shape = TileShape.SW, X = (sbyte)(-Y - 1), Y = X };
case TileShape.SE:
return new Tile { Shape = TileShape.NE, X = (sbyte)(-Y - 1), Y = X };
case TileShape.SW:
return new Tile { Shape = TileShape.SE, X = (sbyte)-Y, Y = X };
default:
throw new NotSupportedException();
}
}

public override int GetHashCode() => (X << 17) + (Y << 3) + (int)Shape;

public bool Equals(Tile tile) => X == tile.X && Y == tile.Y && Shape == tile.Shape;

public override bool Equals(object obj) => obj is Tile tile && Equals(tile);

public int CompareTo(Tile other)
{
if (X != other.X) return X.CompareTo(other.X);
if (Y != other.Y) return Y.CompareTo(other.Y);
return Shape.CompareTo(other.Shape);
}

public override string ToString() => \$"({X},{Y},{Shape})";
}

enum TileShape : byte
{
Anticlockwise,
Clockwise,
NE,
SW,
NW,
SE
}

{
public TileSet(IEnumerable<Tile> tiles)
{
// Canonicalise
var ordered = _Canonicalise(new HashSet<Tile>(tiles));
int h = 1;
foreach (var tile in ordered) h = h * 37 + tile.GetHashCode();
_HashCode = h;

#if SUPERLIGHT

// Since we normalise to have minimum X and Y of 0, we can use unsigned coordinates.
// And since we're looking at connected graphs of on the order of 20 items, 6 bits per coordinate is plenty.
_Items = ordered.Select(tile => (short)((tile.X << 9) + (tile.Y << 3) + (int)tile.Shape)).ToArray();

#else

_Items = new HashSet<Tile>(ordered);

#endif
}

{
int ac = tiles.Count(tile => tile.Shape == TileShape.Anticlockwise);
int c = tiles.Count(tile => tile.Shape == TileShape.Clockwise);

if (ac < c) return _CanonicaliseRot(tiles);
if (ac > c) return _CanonicaliseRot(tiles.Select(tile => tile.Flip()));

return _Min(_CanonicaliseRot(tiles), _CanonicaliseRot(tiles.Select(tile => tile.Flip())));
}

{
for (int i = 0; i < tiles1.Count; i++)
{
int cmp = tiles1[i].CompareTo(tiles2[i]);
if (cmp < 0) return tiles1;
if (cmp > 0) return tiles2;
}

return tiles1;
}

{
//   Rotation anticlockwise by 90 maps NE -> NW -> SW -> SE -> NE. We rotate to get one of these necklaces (in rank order, not exact values):
//     Necklaces:
//     SE NE NW SW
//     0  0  0  0    ** Four positions to consider
//     1  0  0  0
//     1  0  1  0    ** Two positions to consider
//     1  1  0  0
//     1  1  1  0
//     2  0  0  1
//     2  0  1  0
//     2  0  1  1
//     2  1  0  0
//     2  1  0  1
//     2  1  1  0
//     2  1  2  0
//     2  2  0  1
//     2  2  1  0
//     3  0  1  2
//     3  0  2  1
//     3  1  0  2
//     3  1  2  0
//     3  2  0  1
//     3  2  1  0

int se = tiles.Count(tile => tile.Shape == TileShape.SE);
int ne = tiles.Count(tile => tile.Shape == TileShape.NE);
int nw = tiles.Count(tile => tile.Shape == TileShape.NW);
int sw = tiles.Count(tile => tile.Shape == TileShape.SW);
var sorted = new int[] { se, ne, nw, sw }.Distinct().OrderBy(x => x);
var index = 1000 * sorted.IndexOf(se) + 100 * sorted.IndexOf(ne) + 10 * sorted.IndexOf(nw) + sorted.IndexOf(sw);
switch (index)
{
case 0:
// All four positions need to be considered
var best = _Translate(tiles);
best = _Min(best, _Translate(tiles.Select(tile => tile.Rot())));
best = _Min(best, _Translate(tiles.Select(tile => tile.Rot().Rot())));
best = _Min(best, _Translate(tiles.Select(tile => tile.Rot().Rot().Rot())));
return best;

case 101:
// Two options need to be considered;
return _Min(_Translate(tiles.Select(tile => tile.Rot())), _Translate(tiles.Select(tile => tile.Rot().Rot().Rot())));

case 1010:
// Two options need to be considered;
return _Min(_Translate(tiles), _Translate(tiles.Select(tile => tile.Rot().Rot())));

case 1000:
case 1100:
case 1110:
case 2001:
case 2010:
case 2011:
case 2100:
case 2101:
case 2110:
case 2120:
case 2201:
case 2210:
case 3012:
case 3021:
case 3102:
case 3120:
case 3201:
case 3210:
// Already in the canonical rotation.
return _Translate(tiles);

case    1:
case 1001:
case 1101:
case   12:
case  102:
case  112:
case 1002:
case 1012:
case 1102:
case 1202:
case 2012:
case 2102:
case  123:
case  213:
case 1023:
case 1203:
case 2013:
case 2103:
// Needs one rotation.
return _Translate(tiles.Select(tile => tile.Rot()));

case   10:
case   11:
case 1011:
case  120:
case 1020:
case 1120:
case   21:
case  121:
case 1021:
case 2021:
case  122:
case 1022:
case 1230:
case 2130:
case  231:
case 2031:
case  132:
case 1032:
// Needs two rotations.
return _Translate(tiles.Select(tile => tile.Rot().Rot()));

case  100:
case  110:
case  111:
case 1200:
case  201:
case 1201:
case  210:
case 1210:
case  211:
case  212:
case 1220:
case  221:
case 2301:
case 1302:
case 2310:
case  312:
case 1320:
case  321:
// Needs three rotations.
return _Translate(tiles.Select(tile => tile.Rot().Rot().Rot()));

default:
throw new NotSupportedException("Case analysis failed");
}
}

{
int minX = tiles.Min(tile => tile.X);
int minY = tiles.Min(tile => tile.Y);
return tiles.
Select(tile => new Tile { Shape = tile.Shape, X = (sbyte)(tile.X - minX), Y = (sbyte)(tile.Y - minY) }).
OrderBy(tile => tile).
ToList();
}

#if SUPERLIGHT

public int Count => _Items.Length;

public IEnumerator<Tile> GetEnumerator()
{
foreach (var encoded in _Items)
{
yield return new Tile { X = (sbyte)((encoded >> 9) & 0x3f), Y = (sbyte)((encoded >> 3) & 0x3f), Shape = (TileShape)(encoded & 0x7) };
}
}

#else

public int Count => _Items.Count;

public IEnumerator<Tile> GetEnumerator() => _Items.GetEnumerator();

public bool Contains(Tile tile) => _Items.Contains(tile);

#endif

IEnumerator IEnumerable.GetEnumerator() => GetEnumerator();

public override int GetHashCode() => _HashCode;

public bool Equals(TileSet tileset) => tileset != null && tileset.Count == Count && tileset._HashCode == _HashCode && _Items.SequenceEqual(tileset._Items);

public override bool Equals(object obj) => obj is TileSet tileset && Equals(tileset);
}
}

static class Extensions
{
internal static int IndexOf<T>(this IEnumerable<T> elts, T elt)
where T : IEquatable<T>
{
int idx = 0;
foreach (var item in elts)
{
if (item.Equals(elt)) return idx;
idx++;
}
return -1;
}
}
}