用 Mathematica 搜索生命游戏中的静物（二）

之前写的那篇《用 Mathematica 搜索生命游戏中的静物》，由于自己写的部分多了一点（虽然关键的一步还是用 Mathematica 自带的 FindPath 函数），特别慢，还特别耗内存。果然对我这种完全不懂算法的人，就应该把所有的事情都交给 Mathematica 才对。

不过这里写出来的还是特别慢，和原来相比快不了多少，只是内存不再会爆炸。

我在 Mathematica StackExchange 上提了一个问题：Searching for still lifes in Conway's Game of Life。大家如果有什么改进代码的建议，欢迎去那里回答。

---

生命游戏里的一个（大小有限的）图样，可以用一个由布尔值组成的二维数组来表示，比如说，Array[b, {w, h}]，这里每个 b[i,j] 表示一个细胞。静物指的是稳定的图样，于是它需要满足下面的条件：

• 每个活细胞周围的活细胞的个数必须是2或者3，
• 每个死细胞周围的活细胞的个数不能是3。

这些条件可以写成这样的 Mathematica 代码：

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}]);

此外，这个图样是有限的，边界之外的细胞总是死的。于是，整个图样要满足的条件可以写成：

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

于是，我们把找静物的问题转化成了一个布尔可满足性问题。Mathematica 里有个叫 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}]]][[1]], {w, h}];

然后我们就可以用这个 SearchStillLife 函数来找静物。不过，它出奇地慢，找个10乘10大小的静物都要花个6秒钟：

ArrayPlot[Boole@SearchStillLife[10, 10], Mesh -> All] // AbsoluteTiming

找出来的静物长这样

不过在这里要提高速度并不难：前面数一个细胞周围的活细胞的个数的时候，用到了 BooleanCountingFunction 函数。它可以指定返回的结果的形式。我测试了各种可能的形式，最快的是 "BFF"，或者说，返回一个 BooleanFunction 对象。于是，我们可以把前面的 NeighborCount 函数改写成：

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

现在再找一遍10乘10大小的静物，只需要0.5秒。找出来的结果没有变化。

不过在找比较大的静物的时候，它还是很慢。比如说，找个30乘30的静物，竟然花了 1246 秒：

ArrayPlot[Boole@SearchStillLife[30, 30], Mesh -> All] // AbsoluteTiming

找出来的静物长这样

然后我就不知道怎样进一步加快它了。

另一个问题是这个代码里没有涉及到任何随机性，因此对同样的输入，找出来的静物也是同样的。一个简单的解决方案是把它异或上一个随机的数组：

SearchRandomStillLife[w_, h_] := 
 Block[{r = RandomChoice[{True, False}, {w, h}]},
  MapThread[Xor, {r,
    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}]]][[1]], {w, h}]}, 2]]

现在它花的时间取决于异或上的那个随机数组。有时很慢，有时会快一些。慢的情况好像比较多。找下面这个10乘10的静物花了1.5秒：

SeedRandom[233];
ArrayPlot[Boole@SearchRandomStillLife[10, 10], Mesh -> All] // AbsoluteTiming

找出来的静物长这样

用这种方法找更大的静物的时候，就会变得特别慢。我试着找了以下30乘30的静物，花了……算了两个多小时也没算出来，我懒得继续算下去了。

有什么能够提速的办法欢迎留言。

---

顺便，一个相关的问题是：在指定的长宽之内，可能的静物有多少个？

如果把旋转、平移或反射之后相同的哪些看成不同的静物的话，代码并不难写，把前面的 SatisfiabilityInstances 换成 SatisfiabilityCount 就行：

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

在长和宽相等的时候，算出来前面几项是：1, 2, 12, 83, 417, 3928, 127425……这就是 OEIS 上的数列 A136278。