# Searching for still lifes in Conway's Game of Life

I'm trying to search for still lifes in Conway's Game of Life.

An object in Conway's Game of Life can be represented by an array of booleans, say, Array[b, {w, h}]. A still life is a stable object, so it needs to satisfy the following conditions:

• For each living cell, the number of its living neighbors must be 2 or 3,
• For each dead cell, the number of its living neighbors must NOT be 3.

Written in Mathematica code, these conditions are:

NeighborCount[k_, {i_, j_}] :=
BooleanCountingFunction[{k},
Delete[Catenate[Array[b, {3, 3}, {i, j} - 1]], 5]];

StillLifeCondition[i_, j_] :=
(b[i, j] && NeighborCount[{2, 3}, {i, j}]) ||
(! b[i, j] && ! NeighborCount[{3}, {i, j}]);


Moreover, all cells outside the boundary are dead. So the condition for the whole object can be written as:

Array[StillLifeCondition[##] /.
b[i_, j_] /; i < 1 || i > w || j < 1 || j > h :> False &,
{w, h} + 2, 0, And]


Now this is a Boolean satisfiability problem, so we can apply SatisfiabilityInstances:

SearchStillLife[w_, h_] :=
ArrayReshape[
SatisfiabilityInstances[
Array[StillLifeCondition[##] /.
b[i_, j_] /; i < 1 || i > w || j < 1 || j > h :> False &,
{w, h} + 2, 0, And],
Catenate[Array[b, {w, h}]]][], {w, h}];


But it is slow. Searching for a still life of size 10x10 takes 6 seconds:

ArrayPlot[Boole@SearchStillLife[10, 10], Mesh -> All] // AbsoluteTiming We can specify the return form of BooleanCountingFunction. I tested all the possible forms. Among them, "BFF" is the fastest. So we can rewrite the NeighborCount function:

NeighborCount[k_, {i_, j_}] :=
BooleanCountingFunction[{k},
Delete[Catenate[Array[b, {3, 3}, {i, j} - 1]], 5], "BFF"];


Now searching for a still life of size 10x10 takes ~0.5 seconds.

But it is still very slow for large sizes. Searching for a still life of size 30x30 takes 1246 seconds:

ArrayPlot[Boole@SearchStillLife[30, 30], Mesh -> All] // AbsoluteTiming How can I speed up the code?

Another problem is that there is no randomness in the code, so the result is completely determined by the input. The simplest way to solve it is to Xor it with a random array:

SearchRandomStillLife[w_, h_] :=
Block[{r = RandomChoice[{True, False}, {w, h}]},
ArrayReshape[
SatisfiabilityInstances[
Array[StillLifeCondition[##] /.
b[i_, j_] /; i < 1 || i > w || j < 1 || j > h :> False /.
b[i_, j_] :> Xor[b[i, j], r[[i, j]]] &,
{w, h} + 2, 0, And],
Catenate[Array[b, {w, h}]]][], {w, h}]}, 2]]


Now the timing depends on this random array. Usually it is slower. The following takes 1.5 seconds.

SeedRandom;
ArrayPlot[Boole@SearchRandomStillLife[10, 10], Mesh -> All] // AbsoluteTiming Is there a better way to generate random result?

### Another algorithm:

Here is the code I wrote in the first version of this question. It is extremely slow, but works for oscillators, and generates random results:

(* An association that, given the first two rows of an 3x3 grid, and \
the next generation of the central cell, return a list of all \
possible third rows. *)
a = Merge[
k[#[[;; 2]], CellularAutomaton["GameOfLife", #][[2, 2]]] -> #[[
3]] & /@ Tuples[Tuples[{0, 1}, 3], 3], Identity];

(* Given two rows of cells and the next generation of the second row, \
return a list of all possible third rows. *)
NextRow[x_, y_, z_] :=
FindPath[Flatten@
MapIndexed[Outer[v[#2 - 1, Most@#] -> v[#2, Rest@#1] &, ##, 1] &,
Lookup[a,
k, {First[
Partition[{x, y}, {2, 3}, {2, 1}, {{1, -2}, {-1, 2}}, 0]],
z}], {}]], v[0, {0, 0}], v[Length@y, {0, 0}], Length@y + 2,
All][[;; , 2 ;;, 2, 1]];

(* Search for an oscillator with width w, height h, and period p. *)
SearchOscillators[w_, h_, p_] := Module[{o = Table[0, w + 2], i, f},
i = NextRow[o, o, o];
f[l : {___, x_, y_}] :=
If[Length@l == h + 2,
If[SameQ[##, o] & @@ y && SubsetQ[i, x], l, Missing["NotFound"]],
Catch[Do[
If[MissingQ[#], , Throw[#]] &[f[Append[l, z]]], {z,
Tuples[RandomSample /@
Missing["NotFound"]]];
ArrayPlot[#, Mesh -> All] & /@
Transpose[
NestWhile[f[{Table[o, p], RandomChoice[i, p]}] &,
Missing["NotFound"], MissingQ]]]

• Maybe (big maybe) you could use some flavor of genetic algorithm or (even bigger maybe) reinforcement learning (think of AlphaGo). I'd actually be really interest to learn if this problem is similar to playing Go in a complexity sense. – Sascha Dec 5 '17 at 17:59
• Please consider making a post at community.wolfram.com We have a relevant group Wolfram Science – Vitaliy Kaurov Sep 4 '18 at 13:44

I don't know what algorithm Mathematica uses in SatisfiabilityInstances, but many SAT solvers are based on the DPLL algorithm, which takes conjunctive normal forms (CNF) as input.

So I tried converting the condition to a CNF before applying SatisfiabilityInstances, i.e., rewriting the StillLifeCondition in the question as:

StillLifeCondition[i_, j_] :=
BooleanConvert[
(b[i, j] && NeighborCount[{2, 3}, {i, j}]) ||
(! b[i, j] && ! NeighborCount[{3}, {i, j}]), "CNF"];


Now the complete code becomes:

NeighborCount[k_, {i_, j_}] :=
BooleanCountingFunction[{k},
Delete[Catenate[Array[b, {3, 3}, {i, j} - 1]], 5]];

StillLifeCondition[i_, j_] :=
BooleanConvert[
(b[i, j] && NeighborCount[{2, 3}, {i, j}]) ||
(! b[i, j] && ! NeighborCount[{3}, {i, j}]), "CNF"];

SearchStillLife[w_, h_] :=
Block[{r = RandomChoice[{True, False}, {w, h}]},
ArrayReshape[
SatisfiabilityInstances[
Array[StillLifeCondition[##] /.
b[i_, j_] /; i < 1 || i > w || j < 1 || j > h :> False /.
b[i_, j_] :> Xor[b[i, j], r[[i, j]]] &,
{w, h} + 2, 0, And],
Catenate[Array[b, {w, h}]]][], {w, h}]}, 2]]


(If I don't Xor it with a random array, it will only return an empty still life.)

It is much faster. Finding a still life of size 30x30 only takes 9 seconds:

ArrayPlot[Boole@SearchStillLife[30, 30], Mesh -> All] // AbsoluteTiming And size 50x50 takes 29 seconds:

ArrayPlot[Boole@SearchStillLife[50, 50], Mesh -> All] // AbsoluteTiming As for oscillators, we can add the time as an additional dimension, i.e., consider the pattern as a 3-dimensional array Array[b, {p, w, h}], where p is the period. So we get:

NeighborCount[k_, {t_, i_, j_}] :=
BooleanCountingFunction[{k},
Delete[Flatten[Array[b, {1, 3, 3}, {t, i - 1, j - 1}]], 5]];

OscillatorCondition[t_, i_, j_] :=
BooleanConvert[
(b[t, i, j] && (b[t + 1, i, j] ⧦
NeighborCount[{2, 3}, {t, i, j}])) ||
(! b[t, i, j] && (b[t + 1, i, j] ⧦
NeighborCount[{3}, {t, i, j}])), "CNF"];

SearchOscillator[p_, w_, h_] :=
Block[{r = RandomChoice[{True, False}, {p, w, h}]},
Xor, {r,
ArrayReshape[
SatisfiabilityInstances[
Array[OscillatorCondition[##] /.
{b[t_, i_, j_] /; i < 1 || i > w || j < 1 || j > h :> False,
b[0, i_, j_] :> b[p, i, j]} /.
b[t_, i_, j_] :> Xor[b[t, i, j], r[[t, i, j]]] &,
{p, w + 2, h + 2}, 0, And],
Flatten[Array[b, {p, w, h}]]][], {p, w, h}]}, 3]]


It is still slow, but no longer unacceptably slow. Finding this oscillator of period 3 and size 16x16 takes 133 seconds:

ArrayPlot[#, Mesh -> All] & /@
Boole@SearchOscillator[3, 16, 16] // AbsoluteTiming I put the codes in a package: https://github.com/AlephAlpha/LifeFind

• Thank you for self-answering to the benefit of the community. – Mr.Wizard Jul 18 '18 at 7:56