Игра «Жизнь», придуманная Джоном Конвеем в 1970 году, до сих пор не теряет популярности. В 1984 году американский математик Билл Госпер опубликовал статью с алгоритмом, позже названным HashLife, ускоряющим симуляцию в триллионы и больше раз. Описание алгоритма и его реализации довольно сложны для понимания и отладки. Я хочу вам представить его версию, реализация которой на питоне влезает на один экран, а суть умещается в 40 строк кода. Если вам это интересно, а также то, как HashLife связан с динамическим программированием, персистентными структурами и системой контроля версий Git, добро пожаловать под кат.

Игра «Жизнь» – клеточный автомат на плоскости, в котором на каждом шаге времени клетка выживает при двух или трёх соседях, мёртвая клетка рождается при трёх соседях, иначе клетка остаётся или становится мёртвой.

Код на Питоне, эмулирующий игру
class SimpleLife:
  def __init__(self, points):
    self.field = set(points)

  def _evolve(self, x, y):
    own = 1 if (x,y) in self.field else 0
    alive = sum((x-1+i, y-1+j) in self.field for i in range(3) for j in range(3)) - own
    return alive == 3 or (alive == 2 and own == 1)
    
  def step(self):
    newfield = set()
    for x,y in self.field:
      for i in range(3):
        for j in range(3):
          if self._evolve(x-1+i, y-1+j):
            newfield.add((x-1+i, y-1+j))
    self.field = newfield

Код можно ускорить в сотни и тысячи раз, и даже в миллионы раз на GPU, но скорость каждого шага пропорциональна как минимум количеству живых клеток. Алгоритм Билла Госпера позволяет эмулировать миллиарды клеток на миллиарды шагов за секунды, переиспользуя результаты вычислений для повторяющихся частей. Про HashLife есть много информации в интернете, в том числе, и на Хабре, но я хочу рассказать немного отличающуюся версию, которую нигде не встречал. Основная идея – рекурсивное разбиение поля не на 4, а на 9 частей.

Разбиение поля на дерево

Представим, что наше поле – квадрат размера 3^n\times 3^n для некоторого большого n. Мы можем представить его в виде дерева: разобьём его на 9 квадратов размера 3^{n-1} \times 3^{n-1}, выложенных в решетку 3x3. Далее каждый квадрат рекурсивно разобьём на все меньшие квадраты, пока не дойдем до квадрата размера 1x1, то есть, одной клетки. Так как многие квадраты будут совпадать, мы можем их дедуплицировать: назначить одинаковым квадратам один id, и в узле дерева держать только массив 3x3 из девяти id квадратов нижнего уровня. Зарезервируем id=0,1 для мертвых и живых клеток поля, тогда узел дерева самого нижнего уровня будет содержать кусок поля 3x3. Например:

При этом можем дедуплицировать квадраты снизу вверх, от менших к большим: два узла дерева совпадают, если совпадают все девять id их дочерних узлов. Это первый из двух трюков в HashLife, позводяющий компактно хранить поля размером миллиард на миллиард клеток и больше.

код на питоне с примером создания пустого узла
class Hashlife:
    def __init__(self):
        self.node_to_id = {}
        self.nodes = [None, None] # reserve for leaves

    def get_id(self, node) -> int:
        if node not in self.node_to_id:
            self.node_to_id[node] = len(self.nodes)
            self.nodes.append(node)
        return self.node_to_id[node]

    def make_empty(self, size):
        e = 0 if size == 3 else self.make_empty(size//3)
        return self.get_id(((e,e,e),(e,e,e),(e,e,e)))

Связь с Git

Мы от двумерного поля перешли к базе данных узлов, каждый из которых либо лист, либо ссылается на дочерние узлы; в этой реализации пока нету никакой специфики игры. Такая структура часто встречается, самый известный пример – система контроля версий Git. Он хранит файлы и директории как узлы делева, которые ссылаются на id дочерних. Только вместо числовых id, они используют хеш данных:

id     name       type     data
--------------------------------------------
ca2f5  main.py    file   print("Hello World!")
dd35a  readme.md  file   Sample python project
3af4e  src        dir    [ ca2f5 ]
69fed  .          dir    [ dd35a, 3af4e ]

Такая структура позволяет гиту эффективно создавать новые коммиты, отличающиеся от родителей только несколькими файлами, переиспользуя уже существующие объекты. Тогда, если вы знаете id объекта верхнего уровня (в гите они называются ref-ами), вы можете быстро восстановить всё дерево. Это классическая персистентная структура данных, где вместо модификации текущего состояния создаётся новое, переиспользующее неизменившиеся части старого, и все состояния хранятся вместе. Мы эту идею и будем использовать для HashLife.

Эволюция

Возьмем узел дерева первого уровня: квадрат 3x3. Мы можем однозначно понять, как эволюционирует его центральная клетка за 1 шаг:

    def evolve3x3(self, node) -> int:
        own = node[1][1]
        alive = sum(node[i][j] for i in range(3) for j in range(3)) - own
        return 1 if alive == 3 or (alive == 2 and own == 1) else 0

Давайте теперь возьмем узел второго уровня, представляющий квадрат 9x9 клеток. За один шаг мы можем эволюционировать его внутренний подквадрат 7x7: каждый из 49 квадратов 3x3 мы превратим в узел первого уровня, и рекурсивно вызовем evolve3x3. Если мы повторим эту процедуру еще дважды, для квадрата 7x7 и потом 5x5, то в конце мы получим квадрат 3x3, из которого мы можем создать узел 1-го уровня, являющимся эволюцией центра нашего узла второго уровня через три шага.

Точно также мы можем для любого узла 3^n \times 3^n получить эволюцию его центра 3^{n-1} \times 3^{n-1} через 3^{n-1} шагов: возьмем девять его дочерних подузлов, и в каждом подузле возьмем девять его под-подузлов, расположив их в решётке 9x9. Каждый из 49 квадратов размера 3x3 в решетке превратим в узел размера 3^{n-1} \times 3^{n-1}, и вызовем для него рекурсивно evolve, получив 49 узлов 3^{n-2} \times 3^{n-2} после 3^{n-2} шагов в форме решетки 7x7. Повторим еще дважды, получив решетки 5x5 и 3x3, и из последней создадим узел 3^{n-1} \times 3^{n-1} после 3*3^{n-2} = 3^{n-1} шагов.

Если мы один раз вычислили эволюцию, мы можем запомнить результат в словаре, и в следующий раз взять оттуда, не вычисляя заново. Это второй трюк в HashLife, позволяющий делать эволюцию на миллиарды поколений за секунды! Это классический пример динамического программирования: мы разбиваем задачу на меньшие подзадачи, часто рекурсивно, и запоминаем результат для решенных подзадач, чтобы взять его из кеша. когда мы встречаем уже виденную раньше подзадачу.

Еще один технический трюк, упрощающий реализацию: при эволюции решетки 9x9 -> 7x7 -> 5x5, левые верхние ячейки квадратов 3x3 используются только один раз. Поэтому, вместо того, чтобы создавать колию решетки 7x7, мы можем записывать результат в самy решетку, в первые ее 49 клеток.

код эволюции
class Hashlife:
  def __init__(self):
      self.node_to_id = {}
      self.nodes = [None, None]
      self.evolved = {}

  def get_id(self, node) -> int:
      if node not in self.node_to_id:
          self.node_to_id[node] = len(self.nodes)
          self.nodes.append(node)
      return self.node_to_id[node]
    
  def evolve(self, node) -> int:
      if node[1][1] in [0,1]:
        own = node[1][1]
        alive = sum(node[i][j] for i in range(3) for j in range(3)) - own
        return 1 if alive == 3 or (alive == 2 and own == 1) else 0

      if node not in self.evolved:
        field = [[0]*9 for _ in range(9)]
        def fget(x, y): 
          return (tuple(field[x][y:y+3]), tuple(field[x+1][y:y+3]), tuple(field[x+2][y:y+3]))
        def fset(x, y, node):
          field[x][y:y+3],field[x+1][y:y+3],field[x+2][y:y+3] = node

        for i in range(3):
          for j in range(3):
            fset(i*3, j*3, self.nodes[node[i][j]])
        for _ in range(3):
          for x in range(7):
            for y in range(7):
              field[x][y] = self.evolve(fget(x, y))
        self.evolved[node] = self.get_id(fget(0, 0))

      return self.evolved[node]

Это и есть основа алгоритма HashLife, занимающая всего 35 строчек кода на Питоне (сравните с классической реализацией алгоритма, имеющей десятки строчек типа ad = life(m.a.a, m.a.b, m.b.a, m.a.c, m.a.d, m.b.c, m.c.a, m.c.b, m.d.a) ).

Эволюция корня

Эволюция корня дерева, то есть, узла, представляющего всё поле 3^n\times 3^n, сложнее: нам не достаточно получить его центр из метода evolve; наоборот, нам надо в худшем случае увеличить поле в три раза (куда докатится волна изменений за 3^n шагов). Самый простой способ – создадим поле 5x5 ячеек, и поместим наш узел в его центр, заполнив оставшиеся ячейки пустыми узлами нужного размера, и запустим только один шаг evolve 5x5->3x3, и центральный квадрат и будет полем размера 3^{n+1} \times 3^{n+1}, являющимся результатом эволюции всего поля за 3^n шагов. То есть, с каждой эволюцией поля его размер и количество шагов за одну эволюцию будет экспоненциально увеличиваться в три раза: 3, 9, 27, 81, итд.

Вот полный код с проверкой результатов:

полный код на Питоне
class Hashlife:
    def __init__(self):
        self.node_to_id = {}
        self.nodes = [None, None]
        self.evolved = {}
        self.population = {}

    def get_id(self, node) -> int:
        if node not in self.node_to_id:
            self.node_to_id[node] = len(self.nodes)
            self.nodes.append(node)
        return self.node_to_id[node]
    
    def evolve(self, node) -> int:
        if node[1][1] in [0,1]:
            own = node[1][1]
            alive = sum(node[i][j] for i in range(3) for j in range(3)) - own
            return 1 if alive == 3 or (alive == 2 and own == 1) else 0

        if node not in self.evolved:
            field = [[0]*9 for _ in range(9)]
            fget = lambda x, y: (tuple(field[x][y:y+3]), tuple(field[x+1][y:y+3]), tuple(field[x+2][y:y+3]))
            def fset(x, y, node):
                field[x][y:y+3],field[x+1][y:y+3],field[x+2][y:y+3] = node

            for i in range(3):
                for j in range(3):
                    fset(i*3, j*3, self.nodes[node[i][j]])
            for _ in range(3):
                for x in range(7):
                    for y in range(7):
                        field[x][y] = self.evolve(fget(x, y))
            self.evolved[node] = self.get_id(fget(0, 0))

        return self.evolved[node]
    
    def evolve_root(self, node):
        e = self.make_empty(self.get_size(node))
        field = [[e]*5 for _ in range(5)]
        field[2][2] = self.get_id(node)
        fget = lambda x, y: (tuple(field[x][y:y+3]), tuple(field[x+1][y:y+3]), tuple(field[x+2][y:y+3]))

        for x in range(3):
            for y in range(3):
                field[x][y] = self.evolve(fget(x, y))
        return fget(0, 0)

    def make_empty(self, size):
        e = 0 if size == 3 else self.make_empty(size//3)
        return self.get_id(((e,e,e),(e,e,e),(e,e,e)))

    def _makefld(self, field: set, x: int, y: int, size: int) -> int:
        if size == 1: return 1 if (x,y) in field else 0
        field = {(p,q) for p,q in field if x <= p < x+size and y <= q < y+size}
        if not field: return self.make_empty(size)
        size //= 3
        node = tuple(tuple(self._makefld(field, x + size*i, y+size*j, size) for i in range(3)) for j in range(3))
        return self.get_id(node)
    
    def init(self, field: set, size: int):
        field = {(x + size//2, y + size//2) for x,y in field}
        return self.nodes[self._makefld(field, 0, 0, size)]
    
    def get_size(self, node):
        if node[1][1] in [0,1]: return 3
        return 3 * self.get_size(self.nodes[node[1][1]])
    
    def count_cells(self, node):
        if node[1][1] in [0,1]: return sum(sum(row) for row in node)
        if node not in self.population:
            count = sum(self.count_cells(self.nodes[node[i][j]]) for i in range(3) for j in range(3))
            self.population[node] = count
        return self.population[node]

class SimpleLife:
    def __init__(self, points):
        self.field = set(points)

    def _evolve(self, x, y):
        own = 1 if (x,y) in self.field else 0
        alive = sum((x-1+i, y-1+j) in self.field for i in range(3) for j in range(3)) - own
        return alive == 3 or (alive == 2 and own == 1)
    
    def step(self):
        newfield = set()
        for x,y in self.field:
            for i in range(3):
                for j in range(3):
                    if self._evolve(x-1+i, y-1+j):
                        newfield.add((x-1+i, y-1+j))
        self.field = newfield


life = SimpleLife([(1,0), (2,0), (0,1), (1,1), (1,2)])
for i in range(30000):
    if i in {0, 9, 36, 117, 360, 1089, 3276, 9837, 29520}:
        print(i, len(life.field))
    life.step()


life = Hashlife()
node = life.init({(1,0), (2,0), (0,1), (1,1), (1,2)}, 3**2)
steps = 0
for i in range(9):
    print(steps, life.count_cells(node))
    steps += life.get_size(node)
    node = life.evolve_root(node)

Если его запустить на тестовой конфигурации, развивающейся хаотично первые 1000 шагов, а потом превращающейся в почти статичную картинку с несколькими осцилляторами и движущимися глайдерами, то можно заметить, что стандартная имплементация симулирует каждый шаг с примерно одной скоростью, зависящей от количества живых клеток, а HashLife сначала работает медленно, пока конфигурация развивается хаотично, а потом экспоненциально ускоряется, на каждом этапе утраивая и размер поля и количество шагов.

Этот пример является полнофункциональным, тем не менее, для реальной реализации можно добавить несколько вещей: сборку мусора (периодически удалять узлы, которе больше не используются, кстати, гит это тоже делает), и в эволюции запоминать не только центральный узел через n^k шагов, но также через n^{k-1}, 2n^{k-1} шагов. Тогда, комбинируя частичные эволюции на каждом уровне дерева, можно очень быстро переходить на любое произвольное число шагов в эволюции. Также нужна визуализация, где можно динамически менять масштаб и смотреть на произвольные части поля.

Заключение

Алгоритм HashLife является довольно нишевым, применимым в основном, к клеточным автоматам. Тем не менее, он комбинирует сразу несколько концепций из computer science (деревья, рекурсия, динамическое программирование, персистентные структуры данных). Попробуйте на этой структуре данных реализовать визуализацию или, например, подсчет живых клеток – это будет хорошее упражнение на динамическое программирование в самой классической его форме! Но по отдельности эти концепции очень важны в теории и на практике: динамическое программирование используется, например, в git diff, основанном на алгоритме Майерса; разбиение плоскости на дерево (QuadTree в оригинальном алгоритме, 9-ary tree в нашем) используется в обработке изображений и индексации геоданных; персистентные структуры используются, как мы видели, в Git, а также, например, в текстовых редакторах для Undo/Redo.

Ну и, конечно, с исторической точки зрения, это пример идеи, неожиданно оказавшейся настолько изящной и эффективной, что ее полезно знать каждому, кто интересуется алгоритмами и их практическим применением.