Амёбы на 3D-сетке

Условие задачи

Рассмотрим трёхмерную кубическую решётку. Амёба в клетке $(x,y,z)$ может разделиться на три амёбы, которые займут клетки $(x+1,y,z)$, $(x,y+1,z)$ и $(x,y,z+1)$, если все они свободны. Изначально занята только клетка $(0,0,0)$. После $N$ делений на решётке находится $2N+1$ амёб.

Обозначим через $D(N)$ количество различных расположений, достижимых после $N$ делений. Разные последовательности делений, которые приводят к одному расположению, считаются одним вариантом. Из условия известны контрольные значения

$$D(2)=3,\qquad D(10)=44499,\qquad D(20)=9204559704,$$

а последние девять цифр $D(100)$ равны $780166455$. Требуется найти последние девять цифр $D(10000)$. Полное условие доступно на странице Project Euler 763.

Почему не нужно хранить порядок делений

Пусть $s_{x,y,z}\in\{0,1\}$ показывает, происходило ли деление в клетке $(x,y,z)$, а $a_{x,y,z}\in\{0,1\}$ — занята ли эта клетка в конечной конфигурации. Баланс амёб в каждой клетке имеет вид

$$ a_{x,y,z}= \mathbf 1_{(x,y,z)=(0,0,0)}-s_{x,y,z} +s_{x-1,y,z}+s_{x,y-1,z}+s_{x,y,z-1}, $$

где величины с отрицательной координатой считаются нулевыми. Всего произошло ровно $N$ делений, поэтому

$$\sum_{x,y,z\geq0}s_{x,y,z}=N.$$

Удобнее всего увидеть связь между набором делений и итоговым расположением в полиномиальной форме. Введём

$$ P(X,Y,Z)=\sum a_{x,y,z}X^xY^yZ^z, \qquad S(X,Y,Z)=\sum s_{x,y,z}X^xY^yZ^z. $$

Деление в клетке $(x,y,z)$ удаляет один соответствующий моном и добавляет три дочерних. Поэтому оно изменяет $P$ на

$$X^xY^yZ^z(X+Y+Z-1),$$

а для всей конфигурации получается равенство

$$\boxed{P=1+(X+Y+Z-1)S}.$$

Если два набора делений дали одно и то же $P$, то $(X+Y+Z-1)(S_1-S_2)=0$. Кольцо полиномов не имеет делителей нуля, следовательно, $S_1=S_2$. Значит, каждое достижимое расположение однозначно определяется множеством клеток, в которых произошло деление; различные допустимые порядки одного и того же набора отдельно считать не требуется.

Диагональные слои

Разобьём решётку на диагональные слои

$$L_k=\{(x,y,z)\in\mathbb Z_{\geq0}^3:x+y+z=k\}.$$

Все дочерние клетки деления из $L_k$ лежат в $L_{k+1}$. Поэтому процесс можно рассматривать послойно. До выполнения делений очередная клетка следующего слоя может получить амёб от трёх родителей; её промежуточная кратность равна $0$, $1$, $2$ или $3$.

Кратность $3$ недопустима: три амёбы пытаются занять одну клетку. Геометрический разбор переходов показывает, что среди клеток кратности $2$ нельзя допускать три типа локальных фигур:

  • треугольник;
  • Y-образное ветвление;
  • две параллельные части, соединённые поперечиной — «степлер».

Каждая такая фигура после вынужденного продолжения создаёт клетку кратности $3$. Если их исключить, клетки кратности $2$ образуют неветвящуюся ломаную — змейку. С точностью до поворотов и отражений её форма кодируется всего двумя неотрицательными числами $(a,b)$:

  • в состоянии $(a,b)$ находится $a+b+1$ активных клеток змейки;
  • $a$ и $b$ описывают длину и положение последнего сегмента после нормализации симметрий;
  • на границах $a=0$, $a=1$ или $b=0$ некоторые геометрические переходы совпадают, поэтому в формулах появляются коэффициенты $2$, $3$ и $4$.

Это ключевое сжатие задачи: вместо хранения произвольного рисунка целого слоя достаточно хранить пару $(a,b)$.

Двухпараметрическая рекуррентность

Пусть $F_{a,b}(m)$ — количество допустимых завершений из канонического состояния $(a,b)$ при бюджете $m$ дальнейших делений. Начальное условие:

$$F_{0,0}(0)=1.$$

Значения с отрицательным $m$ равны нулю. Чтобы построить состояние с $a+b+1$ активными клетками, необходимо выполнить как минимум

$$1+2+\dots+(a+b)=\frac{(a+b)(a+b+1)}2$$

предыдущих делений. Поэтому действует дополнительное граничное условие

$$ F_{a,b}(m)=0 \quad\text{при}\quad m<\frac{(a+b)(a+b+1)}2. $$

Чтобы не повторять один и тот же сдвиг аргумента, обозначим

$$G_m(p,q)=F_{p,q}(m-p-q-1).$$

Величина $p+q+1$ вычитается потому, что именно столько клеток делится в текущем слое, если следующим становится состояние $(p,q)$. После нормализации поворотов и отражений остаются шесть случаев.

1. Пустая каноническая граница

$$\boxed{F_{0,0}(m)=3G_m(0,0)+3G_m(1,0)}.$$

Коэффициент $3$ соответствует трём равноправным направлениям координатных осей.

2. Граница $a=0$, $b\geq1$

$$ \boxed{ F_{0,b}(m)=G_m(0,b)+G_m(b-1,1)+G_m(0,b+1)+G_m(b,1) }. $$

Этот граничный случай важно выписать отдельно: формально он получается из общей формулы после зануления состояний с отрицательным первым индексом, но без явной записи его легко пропустить при реализации.

3. Граница $(1,0)$

$$ \boxed{ F_{1,0}(m)=G_m(0,0)+4G_m(1,0)+G_m(2,0)+2G_m(1,1) }. $$

4. Случай $a=1$, $b\geq1$

$$ \boxed{ \begin{aligned} F_{1,b}(m)={}&G_m(b,0)+G_m(b+1,0)+G_m(b,1)\\ &+G_m(1,b)+G_m(b+1,1)+G_m(1,b+1). \end{aligned} } $$

5. Случай $a\geq2$, $b=0$

$$ \boxed{ \begin{aligned} F_{a,0}(m)={}&G_m(a-1,0)+2G_m(a,0)+2G_m(a-1,1)\\ &+G_m(a+1,0)+2G_m(a,1). \end{aligned} } $$

6. Общий случай $a\geq2$, $b\geq1$

$$ \boxed{ \begin{aligned} F_{a,b}(m)={}&G_m(a-1,b)+G_m(a,b)+G_m(a-1,b+1)\\ &+G_m(a+b-1,1)+G_m(a,b+1)+G_m(a+b,1). \end{aligned} } $$

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

Первое деление исходной амёбы обязательно. После него остаётся бюджет $N-1$, поэтому искомое количество конфигураций равно

$$\boxed{D(N)=F_{0,0}(N-1)}.$$

Как вычислить динамику

Значения вычисляются по возрастанию $m$. Для каждого $m$ достаточно обрабатывать только пары $(a,b)$, которые проходят треугольное ограничение. При $m\leq9999$ выполняется

$$a+b\leq140,$$

поскольку

$$\frac{140\cdot141}{2}=9870, \qquad \frac{141\cdot142}{2}=10011.$$

Наивно можно сохранить все значения $F_{a,b}(m)$, но это расходует слишком много памяти. Рекуррентность читает состояние $(a,b)$ только со сдвигом $a+b+1$, поэтому для каждой пары достаточно кольцевого буфера длины $a+b+2$. Новое значение перезаписывает элемент, который уже не понадобится.

Для одного $m$ существует $O(m)$ допустимых пар $(a,b)$, а полный расчёт требует $O(N^2)$ операций. Благодаря кольцевым буферам память уменьшается до $O(N^{3/2})$. Все значения для основной задачи можно сразу вычислять по модулю $10^9$: сложение и умножение на небольшие коэффициенты совместимы с взятием остатка.

Эскиз алгоритма без готового кода

Публично я не привожу готовую реализацию, но порядок вычислений можно описать так:

  1. Создать состояние для каждой пары $(a,b)$ с $a+b\leq140$ и выделить кольцевой буфер длины $a+b+2$.
  2. Записать базу $F_{0,0}(0)=1$; остальные ячейки изначально равны нулю.
  3. Перебирать $m$ от $1$ до $N-1$.
  4. Для текущего $m$ найти максимальное $a+b$, удовлетворяющее треугольному ограничению.
  5. Перебрать допустимые пары $(a,b)$ и применить соответствующий из шести случаев рекуррентности.
  6. После каждого сложения сохранять результат по модулю $10^9$.
  7. Вернуть $F_{0,0}(N-1)$.

Проверки

Для небольших $N$ можно независимо перебрать все достижимые множества занятых клеток. Такой прямой перебор для $N=0,1,\ldots,8$ даёт последовательность

$$1,\ 1,\ 3,\ 9,\ 30,\ 99,\ 336,\ 1134,\ 3855.$$

Динамика воспроизводит эту последовательность и все четыре контрольных значения из условия:

  • $D(2)=3$;
  • $D(10)=44499$;
  • $D(20)=9204559704$;
  • $D(100)\equiv780166455\pmod{10^9}$.

Совпадение прямого перебора с динамикой проверяет переходы на маленьких состояниях, а значения $D(20)$ и $D(100)$ дополнительно ловят ошибки в границах и модульной арифметике.

Итог

Главная идея решения — отказаться от перебора последовательностей делений и перейти к диагональным слоям. После исключения невозможных локальных фигур весь сложный рисунок слоя сжимается до состояния змейки $(a,b)$. Треугольная нижняя граница оставляет лишь небольшое число состояний, а кольцевые буферы делают расчёт практичным по памяти.

В результате последние девять цифр числа конфигураций после $10000$ делений равны

$$\boxed{D(10000)\equiv798443574\pmod{10^9}}.$$

Вверх Вниз