Free polyominoes enumeration

Free polyominoes enumeration is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

A Polyomino is a plane geometric figure formed by joining one or more equal squares edge to edge. Free polyominoes are distinct when none is translation, rotation, reflection or glide reflection of another: wp:Polyomino.

Task: generate all the free polyominoes with n cells.

You can visualize them just as a sequence of the coordinate pairs of their cells (rank 5):

[(0, 0), (0, 1), (0, 2), (0, 3), (0, 4)]
[(0, 0), (0, 1), (0, 2), (0, 3), (1, 0)]
[(0, 0), (0, 1), (0, 2), (0, 3), (1, 1)]
[(0, 0), (0, 1), (0, 2), (1, 0), (1, 1)]
[(0, 0), (0, 1), (0, 2), (1, 0), (1, 2)]
[(0, 0), (0, 1), (0, 2), (1, 0), (2, 0)]
[(0, 0), (0, 1), (0, 2), (1, 1), (2, 1)]
[(0, 0), (0, 1), (0, 2), (1, 2), (1, 3)]
[(0, 0), (0, 1), (1, 1), (1, 2), (2, 1)]
[(0, 0), (0, 1), (1, 1), (1, 2), (2, 2)]
[(0, 0), (0, 1), (1, 1), (2, 1), (2, 2)]
[(0, 1), (1, 0), (1, 1), (1, 2), (2, 1)]

But a better basic visualization is using ASCII art (rank 5):

#    ##   #    ##  ##  ###  #     #    #    #    #      #
#    #    ##   ##  #   #    ###   #    ###  ##   ###   ###
#    #    #    #   ##  #    #     ##    #    ##    #    #
#    #    #                        #
#

For a slow but clear solution see this Haskell Wiki page: http://www.haskell.org/haskellwiki/The_Monad.Reader/Issue5/Generating_Polyominoes

Bonus Task: you can create an alternative program (or specialize your first program) to generate very quickly just the number of distinct free polyominoes, and to show a sequence like:

1, 1, 1, 2, 5, 12, 35, 108, 369, 1285, 4655, 17073, 63600, 238591, 901971, 3426576, ...

Number of free polyominoes (or square animals) with n cells: http://oeis.org/A000105

D

Translation of: Haskell

<lang d>import std.stdio, std.range, std.algorithm, std.typecons, std.conv;

alias Coord = byte; alias Point = Tuple!(Coord,"x", Coord,"y"); alias Polyomino = Point[];

/// Finds the min x and y coordiate of a Polyomino. enum minima = (in Polyomino poly) pure @safe =>

   Point(poly.map!q{ a.x }.reduce!min, poly.map!q{ a.y }.reduce!min);

Polyomino translateToOrigin(in Polyomino poly) {

   const minP = poly.minima;
   return poly.map!(p => Point(cast(Coord)(p.x - minP.x), cast(Coord)(p.y - minP.y))).array;

}

enum Point function(in Point p) pure nothrow @safe @nogc

   rotate90  = p => Point( p.y, -p.x),
   rotate180 = p => Point(-p.x, -p.y),
   rotate270 = p => Point(-p.y,  p.x),
   reflect   = p => Point(-p.x,  p.y);

/// All the plane symmetries of a rectangular region. auto rotationsAndReflections(in Polyomino poly) pure nothrow {

   return only(poly,
               poly.map!rotate90.array,
               poly.map!rotate180.array,
               poly.map!rotate270.array,
               poly.map!reflect.array,
               poly.map!(pt => pt.rotate90.reflect).array,
               poly.map!(pt => pt.rotate180.reflect).array,
               poly.map!(pt => pt.rotate270.reflect).array);

}

enum canonical = (in Polyomino poly) =>

   poly.rotationsAndReflections.map!(pl => pl.translateToOrigin.sort().release).reduce!min;

auto unique(T)(T[] seq) pure nothrow {

   return seq.sort().uniq;

}

/// All four points in Von Neumann neighborhood. enum contiguous = (in Point pt) pure nothrow @safe @nogc =>

   only(Point(cast(Coord)(pt.x - 1), pt.y), Point(cast(Coord)(pt.x + 1), pt.y),
        Point(pt.x, cast(Coord)(pt.y - 1)), Point(pt.x, cast(Coord)(pt.y + 1)));

/// Finds all distinct points that can be added to a Polyomino. enum newPoints = (in Polyomino poly) nothrow =>

   poly.map!contiguous.joiner.filter!(pt => !poly.canFind(pt)).array.unique;

enum newPolys = (in Polyomino poly) =>

   poly.newPoints.map!(pt => canonical(poly ~ pt)).array.unique;

/// Generates polyominoes of rank n recursively. Polyomino[] rank(in uint n) {

   static immutable Polyomino monomino = [Point(0, 0)];
   static Polyomino[] monominoes = [monomino]; // Mutable.
   if (n == 0) return [];
   if (n == 1) return monominoes;
   return rank(n - 1).map!newPolys.join.unique.array;

}

/// Generates a textual representation of a Polyomino. char[][] textRepresentation(in Polyomino poly) pure @safe {

   immutable minPt = poly.minima;
   immutable maxPt = Point(poly.map!q{ a.x }.reduce!max, poly.map!q{ a.y }.reduce!max);
   auto table = new char[][](maxPt.y - minPt.y + 1, maxPt.x - minPt.x + 1);
   foreach (row; table)
       row[] = ' ';
   foreach (immutable pt; poly)
       table[pt.y - minPt.y][pt.x - minPt.x] = '#';
   return table;

}

void main(in string[] args) {

   iota(1, 11).map!(n => n.rank.length).writeln;

   immutable n = (args.length == 2) ? args[1].to!uint : 5;
   writefln("\nAll free polyominoes of rank %d:", n);

   foreach (const poly; n.rank)
       writefln("%-(%s\n%)\n", poly.textRepresentation);

}</lang>

Output:

[1, 1, 2, 5, 12, 35, 108, 369, 1285, 4655]

All free polyominoes of rank 5:
#
#
#
#
#

##
#
#
#

#
##
#
#

##
##
#

##
#
##

###
#
#

#
###
#

#
#
##
 #

#
###
 #

#
##
 ##

#
###
  #

 #
###
 #

D: Count Only

Translated and modified from C code: http://www.geocities.jp/tok12345/countomino.txt

<lang d>import core.stdc.stdio: printf; import core.stdc.stdlib: atoi;

__gshared ulong[] g_pnCountNH; __gshared uint[] g_pnFieldCheck, g_pnFieldCheckR; __gshared uint g_nFieldSize, g_nFieldWidth; __gshared uint[4] g_anLinkData; __gshared uint[8] g_anRotationOffset, g_anRotationX, g_anRotationY;

void countMain(in uint n) nothrow {

   g_nFieldWidth = n * 2 - 2;
   g_nFieldSize = (n - 1) * (n - 1) * 2 + 1;
   g_pnCountNH = new ulong[n + 1];

   auto pnField = new uint[g_nFieldSize];
   auto pnPutList = new uint[g_nFieldSize];
   g_pnFieldCheck = new uint[n ^^ 2];
   g_pnFieldCheckR = new uint[n ^^ 2];
   g_anLinkData[0] = 1;
   g_anLinkData[1] = g_nFieldWidth;
   g_anLinkData[2] = -1;
   g_anLinkData[3] = -g_nFieldWidth;

   initOffset(n);

   countSub(n, 0, pnField, pnPutList, 0, 1);

}

void countSub(in uint n, in uint lv, uint[] field, uint[] putlist,

             in uint putno, in uint putlast) nothrow /*@nogc*/ {
   check(field, n, lv);
   if (n == lv) {
       return;
   }

   foreach (immutable uint i; putno .. putlast) {
       immutable pos = putlist[i];
       field[pos] |= 1;
       uint k = 0;
       foreach (immutable uint j; 0 .. 4) {
           immutable pos2 = pos + g_anLinkData[j];
           if (0 <= pos2 && pos2 < g_nFieldSize && !field[pos2]) {
               field[pos2] = 2;
               putlist[putlast + k] = pos2;
               k++;
           }
       }
       countSub(n, lv + 1, field, putlist, i + 1, putlast + k);
       foreach (immutable uint j; 0 .. k)
           field[putlist[putlast + j]] = 0;
       field[pos] = 2;
   }

   foreach (immutable uint i; putno .. putlast) {
       immutable pos = putlist[i];
       field[pos] &= -2;
   }

}

void initOffset(in uint n) nothrow /*@nogc*/ {

   g_anRotationOffset[0] = 0;
   g_anRotationX[0] = 1;
   g_anRotationY[0] = n;
   // 90
   g_anRotationOffset[1] = n - 1;
   g_anRotationX[1] = n;
   g_anRotationY[1] = -1;
   // 180
   g_anRotationOffset[2] = n ^^ 2 - 1;
   g_anRotationX[2] = -1;
   g_anRotationY[2] = -n;
   // 270
   g_anRotationOffset[3] = n ^^ 2 - n;
   g_anRotationX[3] = -n;
   g_anRotationY[3] = 1;

   g_anRotationOffset[4] = n - 1;
   g_anRotationX[4] = -1;
   g_anRotationY[4] = n;
   // 90
   g_anRotationOffset[5] = 0;
   g_anRotationX[5] = n;
   g_anRotationY[5] = 1;
   // 180
   g_anRotationOffset[6] = n ^^ 2 - n;
   g_anRotationX[6] = 1;
   g_anRotationY[6] = -n;
   // 270
   g_anRotationOffset[7] = n ^^ 2 - 1;
   g_anRotationX[7] = -n;
   g_anRotationY[7] = -1;

}

void check(in uint[] field, in uint n, in uint lv) nothrow /*@nogc*/ {

   g_pnFieldCheck[0 .. n ^^ 2] = 0;

   uint x, y;
   outer:
   for (x = n; x < n * 2 - 2; x++)
       for (y = 0; y + x < g_nFieldSize; y += g_nFieldWidth)
           if (field[x + y] & 1)
               break outer;

   immutable uint x2 = n - x;
   foreach (immutable uint i; 0 .. g_nFieldSize) {
       x = (i + n - 2) % g_nFieldWidth;
       y = (i + n - 2) / g_nFieldWidth * n;
       if (field[i] & 1)
           g_pnFieldCheck[x + x2 + y] = 1;
   }

   uint of1;
   for (of1 = 0; of1 < g_pnFieldCheck.length && !g_pnFieldCheck[of1]; of1++) {}

   bool c = true;
   for (uint r = 1; r < 8 && c; r++) {
       for (x = 0; x < n; x++) {
           for (y = 0; y < n; y++) {
               immutable pos = g_anRotationOffset[r] +
                               g_anRotationX[r] * x + g_anRotationY[r] * y;
               g_pnFieldCheckR[pos] = g_pnFieldCheck[x + y * n];
           }
       }

       uint of2;
       for (of2 = 0; of2 < g_pnFieldCheckR.length && !g_pnFieldCheckR[of2]; of2++) {}
       of2 -= of1;
       immutable ed = (of2 > 0) ? (n ^^ 2 - of2) : (n ^^ 2);

       foreach (immutable uint i; of1 .. ed) {
           if (g_pnFieldCheck[i] > g_pnFieldCheckR[i + of2])
               break;
           if (g_pnFieldCheck[i] < g_pnFieldCheckR[i + of2]) {
               c = false;
               break;
           }
       }
   }

   if (c) {
       uint parity;
       if (!(lv & 1)) {
           parity = (lv & 2) >> 1;
           for (x = 0; x < n; x++)
               for (y = 0; y < n; y++)
                   parity ^= (x + y) & g_pnFieldCheck[x + y * n];
           parity &= 1;
       } else
           parity = 0;

       g_pnCountNH[lv]++;
   }

}

int main(in string[] args) {

   immutable n = (args.length == 2) ? (args[1] ~ '\0').ptr.atoi : 11;
   if (n < 1)
       return 1;

   if (n == 1)
       countMain(2);
   else
       countMain(n);

   foreach (immutable i; 1 .. n + 1)
       printf("%llu\n", g_pnCountNH[i]);

   return 0;

}</lang>

Output:

1
1
2
5
12
35
108
369
1285
4655
17073

Output with n=14 (run-time about 36 seconds):

1
1
2
5
12
35
108
369
1285
4655
17073
63600
238591
901971

Haskell

This Haskell solution is relatively slow, it's meant to be readable and as manifestly correct as possible.

Code updated and slightly improved from: http://www.haskell.org/haskellwiki/The_Monad.Reader/Issue5/Generating_Polyominoes <lang haskell>import Data.List (sort) import Data.Set (toList, fromList) import System.Environment (getArgs)

type Coord = Int type Point = (Coord, Coord) type Polyomino = [Point]

-- Finds the min x and y coordiate of a Polyomino. minima :: Polyomino -> Point minima (p:ps) = foldr (\(x, y) (mx, my) -> (min x mx, min y my)) p ps

translateToOrigin :: Polyomino -> Polyomino translateToOrigin p =

   let (minx, miny) = minima p in
       map (\(x, y) -> (x - minx, y - miny)) p

rotate90, rotate180, rotate270, reflect :: Point -> Point rotate90 (x, y) = ( y, -x) rotate180 (x, y) = (-x, -y) rotate270 (x, y) = (-y, x) reflect (x, y) = (-x, y)

-- All the plane symmetries of a rectangular region. rotationsAndReflections :: Polyomino -> [Polyomino] rotationsAndReflections p =

   [p,
    map rotate90 p,
    map rotate180 p,
    map rotate270 p,
    map reflect p,
    map (rotate90 . reflect) p,
    map (rotate180 . reflect) p,
    map (rotate270 . reflect) p]

canonical :: Polyomino -> Polyomino canonical = minimum . map (sort . translateToOrigin) . rotationsAndReflections

unique :: (Ord a) => [a] -> [a] unique = toList . fromList

-- All four points in Von Neumann neighborhood. contiguous :: Point -> [Point] contiguous (x, y) = [(x - 1, y), (x + 1, y), (x, y - 1), (x, y + 1)]

-- Finds all distinct points that can be added to a Polyomino. newPoints :: Polyomino -> [Point] newPoints p =

   let notInP = filter (not . flip elem p) in
       unique . notInP . concatMap contiguous $ p

newPolys :: Polyomino -> [Polyomino] newPolys p = unique . map (canonical . flip (:) p) $ newPoints p

monomino = [(0, 0)] monominoes = [monomino]

-- Generates polyominoes of rank n recursively. rank :: Int -> [Polyomino] rank 0 = [] rank 1 = monominoes rank n = unique . concatMap newPolys $ rank (n - 1)

-- Generates a textual representation of a Polyomino. textRepresentaton :: Polyomino -> String textRepresentaton p =

   unlines  x <- [0 .. maxx - minx
            | y <- [0 .. maxy - miny]]
   where
       maxima :: Polyomino -> Point
       maxima (p:ps) = foldr (\(x, y) (mx, my) -> (max x mx, max y my)) p ps
       (minx, miny) = minima p
       (maxx, maxy) = maxima p

main = do

   print $ map (length . rank) [1 .. 10]
   args <- getArgs
   let n = if null args then 5 else read $ head args :: Int
   putStrLn ("\nAll free polyominoes of rank " ++ show n ++ ":")
   mapM_ putStrLn $ map textRepresentaton $ rank n</lang>

Output:

[1,1,2,5,12,35,108,369,1285,4655]

All free polyominoes of rank 5:
#
#
#
#
#

##
# 
# 
# 

# 
##
# 
# 

##
##
# 

##
# 
##

###
#  
#  

#  
###
#  

# 
# 
##
 #

#  
###
 # 

#  
## 
 ##

#  
###
  #

 # 
###
 # 

Python

Translation of: Haskell

<lang python>from itertools import imap, imap, groupby, chain, imap from operator import itemgetter from sys import argv from array import array

def concat_map(func, it):

   return list(chain.from_iterable(imap(func, it)))

def minima(poly):

   """Finds the min x and y coordiate of a Polyomino."""
   return (min(pt[0] for pt in poly), min(pt[1] for pt in poly))

def translate_to_origin(poly):

   (minx, miny) = minima(poly)
   return [(x - minx, y - miny) for (x, y) in poly]

rotate90 = lambda (x, y): ( y, -x) rotate180 = lambda (x, y): (-x, -y) rotate270 = lambda (x, y): (-y, x) reflect = lambda (x, y): (-x, y)

def rotations_and_reflections(poly):

   """All the plane symmetries of a rectangular region."""
   return (poly,
           map(rotate90, poly),
           map(rotate180, poly),
           map(rotate270, poly),
           map(reflect, poly),
           [reflect(rotate90(pt)) for pt in poly],
           [reflect(rotate180(pt)) for pt in poly],
           [reflect(rotate270(pt)) for pt in poly])

def canonical(poly):

   return min(sorted(translate_to_origin(pl)) for pl in rotations_and_reflections(poly))

def unique(lst):

   lst.sort()
   return map(next, imap(itemgetter(1), groupby(lst)))

1. All four points in Von Neumann neighborhood.

contiguous = lambda (x, y): [(x - 1, y), (x + 1, y), (x, y - 1), (x, y + 1)]

def new_points(poly):

   """Finds all distinct points that can be added to a Polyomino."""
   return unique([pt for pt in concat_map(contiguous, poly) if pt not in poly])

def new_polys(poly):

   return unique([canonical(poly + [pt]) for pt in new_points(poly)])

monomino = [(0, 0)] monominoes = [monomino]

def rank(n):

   """Generates polyominoes of rank n recursively."""
   assert n >= 0
   if n == 0: return []
   if n == 1: return monominoes
   return unique(concat_map(new_polys, rank(n - 1)))

def text_representation(poly):

   """Generates a textual representation of a Polyomino."""
   min_pt = minima(poly)
   max_pt = (max(p[0] for p in poly), max(p[1] for p in poly))
   table = [array('c', ' ') * (max_pt[1] - min_pt[1] + 1)
            for _ in xrange(max_pt[0] - min_pt[0] + 1)]
   for pt in poly:
       table[pt[0] - min_pt[0]][pt[1] - min_pt[1]] = '#'
   return "\n".join(row.tostring() for row in table)

def main():

   print [len(rank(n)) for n in xrange(1, 11)]

   n = int(argv[1]) if (len(argv) == 2) else 5
   print "\nAll free polyominoes of rank %d:" % n

   for poly in rank(n):
       print text_representation(poly), "\n"

main()</lang>

Output:

[1, 1, 2, 5, 12, 35, 108, 369, 1285, 4655]

All free polyominoes of rank 5:
##### 

####
#    

####
 #   

###
##  

###
# # 

###
#  
#   

###
 # 
 #  

### 
  ## 

## 
 ##
 #  

## 
 ##
  # 

## 
 # 
 ## 

 # 
###
 #