Free polyominoes enumeration
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
<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: # # # # # ## # # # # ## # # ## ## # ## # ## ### # # # ### # # # ## # # ### # # ## ## # ### # # ### #
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
<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)))
- 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: ##### #### # #### # ### ## ### # # ### # # ### # # ### ## ## ## # ## ## # ## # ## # ### #