Ториангуляции
Условие задачи
Название Toriangulations склеено из слов torus и triangulations. Берётся прямоугольник $a\times b$, его противоположные стороны склеиваются, и получившийся тор нужно замостить единичными равносторонними треугольниками. Два замощения считаются эквивалентными, если одно можно непрерывно превратить в другое без разрывов и наложений.
Пусть $F(n)$ — число классов замощений всех возможных торов ровно $n$ треугольниками, а
$$G(N)=\sum_{n=1}^{N}F(n).$$
Известно, что $G(6)=14$, $G(100)=8090$ и $G(10^5)\equiv645124048\pmod{1\,000\,000\,007}$. Требуется найти $G(10^9)$ по тому же модулю. Полное условие и исходные иллюстрации находятся на странице Project Euler 780.
Почему я на этой задаче застрял
Первое желание — рисовать треугольники и сравнивать, какие картинки можно непрерывно деформировать друг в друга. Но такой перебор ничего не упрощает: размеры тора вещественные, прямоугольник склеен по двум направлениям, а один и тот же рисунок можно разрезать в разных местах. Искать инвариант непосредственно по картинке оказалось бесполезно — частных случаев становилось только больше.
В итоге задача распалась для меня на три отдельных препятствия:
- понять, какую форму вообще может иметь любое замощение;
- закодировать эту форму целыми числами, несмотря на вещественные размеры тора;
- найти и исправить замощения, которые один и тот же подсчёт видит несколько раз.
Первая действительно полезная зацепка нашлась в обсуждении решивших задачу: ecnerwala предложил перестать классифицировать картинки и считать направления полос. После этого геометрия наконец начала сжиматься до арифметики.
Первая рабочая идея: развернуть тор
Возьмём прямоугольник тора и приклеим к нему бесконечно много копий по горизонтали и вертикали. Замощение продолжится до замощения всей плоскости единичными равносторонними треугольниками. Высота такого треугольника равна
$$h=\frac{\sqrt3}{2}.$$
В плоскости треугольники неизбежно собираются в бесконечные параллельные полосы ширины $h$. Внутри каждой полосы ориентации чередуются: вверх, вниз, вверх, вниз. Причина локальная: каждая сторона должна полностью совпасть со стороной соседа, а в каждой вершине сходятся шесть углов по $60^\circ$. Продолжив один ряд через общие стороны, мы получаем полосу постоянной ширины; её граница вынуждает соседний ряд быть параллельным.
Это был главный перелом. Вместо огромного зоопарка триангуляций остались замкнутые полосы. Из чередования ориентаций сразу следует, что число треугольников $n$ обязательно чётно. Но здесь же спрятана будущая проблема:
- у обычного сдвинутого замощения направление полос единственно;
- у регулярной треугольной решётки полосы можно провести в трёх направлениях.
Второй пункт я сначала недооценил. Именно он позже испортил первую почти готовую формулу.
Как одна полоса обходит тор
Чтобы не путать размеры прямоугольника с параметрами обхода, обозначим стороны тора через $u$ и $v$. За один полный оборот полоса проходит $a$ горизонтальных и $b$ вертикальных периодов. На развёртке её направляющий вектор равен
$$(au,bv),$$
где можно считать $a>0$, $b\geq0$ и $\gcd(a,b)=1$. Взаимная простота здесь не украшение: без неё мы несколько раз описывали бы один более короткий обход.
Пусть после склейки получилось $k$ параллельных полос. В каждой полосе находится
$$\frac{n}{2k}$$
пар треугольников, поэтому $k\mid n/2$, а длина замкнутого обхода равна $n/(2k)$. Теперь геометрия даёт два уравнения. Первое — по площади:
$$uv=\frac{\sqrt3}{4}n.$$
Второе — по длине направляющего вектора:
$$(au)^2+(bv)^2=\left(\frac{n}{2k}\right)^2.$$
Сначала я воспринимал эти равенства только как необходимые проверки. Оказалось, что они ещё и достаточны: суммарная длина всех полос равна $n/2$, а отношение площади тора к этой длине равно $\sqrt3/2$ — ровно ширине одной полосы. Значит, подходящими полосами действительно можно заполнить весь тор.
Простой случай: полосы параллельны стороне
Если $b=0$, взаимная простота заставляет взять $a=1$. Два уравнения сразу дают
$$u=\frac{n}{2k},\qquad v=\frac{\sqrt3}{2}k.$$
Для каждого делителя $k\mid n/2$ существует одно замощение. Полосы могут быть горизонтальными или вертикальными, поэтому вклад равен
$$2\tau\!\left(\frac n2\right),$$
где $\tau(r)$ — число положительных делителей $r$. Для суммы по всем $n\leq N$ положим
$$M=\left\lfloor\frac N2\right\rfloor, \qquad D(X)=\sum_{d=1}^{X}\left\lfloor\frac Xd\right\rfloor.$$
Тогда весь вклад полос, параллельных сторонам, равен
$$S_0(N)=2D(M).$$
Наклонные полосы и квадратный корень
Случай $a,b>0$ поначалу выглядел самым неприятным: размеры $u$ и $v$ вещественные, а считать нужно целочисленные классы. Помогла замена
$$A=au,\qquad B=bv,\qquad m=\frac{n}{2k}.$$
После сложения и вычитания уравнений площади и длины получаются два квадрата:
$$ (A+B)^2=m^2+\sqrt3\,abkm, \qquad (A-B)^2=m^2-\sqrt3\,abkm. $$
Положительное решение существует тогда и только тогда, когда
$$m>\sqrt3\,abk.$$
Равенство невозможно из-за иррациональности $\sqrt3$. Знак $A-B$ даёт два различных отношения сторон, а наклон полосы может быть положительным или отрицательным. Поэтому каждый допустимый набор $(a,b,k,m)$ даёт четыре замощения.
В этот момент стало понятно, зачем была нужна взаимная простота. Зафиксируем произведение $x=ab$. Если
$$x=\prod_{i=1}^{r}p_i^{e_i},$$
то каждая простая степень должна целиком попасть либо в $a$, либо в $b$. Значит, число упорядоченных взаимно простых пар с произведением $x$ равно
$$R(x)=2^{\omega(x)},$$
где $\omega(x)$ — число различных простых делителей $x$. Для фиксированных $k$ и $x$ допустимые $m$ лежат в промежутке
$$\sqrt3\,kx<m\leq\left\lfloor\frac Mk\right\rfloor.$$
Первая почти правильная формула
Обозначим $[y]_+=\max(y,0)$. Подсчёт по направлениям полос приводит к формуле
$$ \boxed{ \begin{aligned} S(N)={}&2D(M)\\ &+4\sum_{k\geq1}\sum_{x\geq1}2^{\omega(x)} \left[ \left\lfloor\frac Mk\right\rfloor- \left\lfloor\sqrt3\,kx\right\rfloor \right]_+. \end{aligned} } $$
Сумма конечна: для данного $k$ достаточно брать
$$ x\leq \left\lfloor \frac{\lfloor M/k\rfloor}{\sqrt3\,k} \right\rfloor. $$
Именно здесь у меня был самый неприятный затуп. Формула выглядела законченной: осевые полосы учтены, наклоны и отражения учтены, условие взаимной простоты учтено. Но контрольные значения получались больше нужных. Я перепроверял коэффициент $4$, границы по $m$ и делители, хотя ошибка была вообще не в арифметике.
Откуда взялся тройной подсчёт
Обычное замощение действительно имеет одно направление полос. Но в регулярной треугольной решётке один и тот же набор треугольников можно разрезать на полосы вдоль каждого из трёх семейств параллельных рёбер. Поэтому $S(N)$ считает такое замощение три раза вместо одного.
Это тот случай, когда ошибка выглядит как неверный коэффициент, а на самом деле вызвана дополнительной симметрией. Для каждого регулярного решёточного замощения нужно вычесть две лишние копии. Если $E(n)$ — число таких замощений с ровно $n$ треугольниками, то
$$G(N)=S(N)-2\sum_{n\leq N}E(n).$$
Целые Эйзенштейна: недостающий кусок
Самостоятельно догадаться до целых Эйзенштейна на этом месте мне не удалось. Подсказка снова нашлась в том же обсуждении: вершины регулярной треугольной решётки образуют именно решётку целых Эйзенштейна. После этого странная геометрическая поправка стала обычной задачей о делителях.
Положим
$$\rho=e^{i\pi/3}=\frac12+\frac{\sqrt3}{2}i, \qquad \mathbb Z[\rho]=\{r+s\rho:r,s\in\mathbb Z\}.$$
После переноса один угол прямоугольника можно поместить в $0$. Из-за периодичности остальные углы тоже должны быть точками этой решётки, поэтому стороны задаются двумя перпендикулярными элементами $U$ и $V$. Уникальность разложения в кольце целых Эйзенштейна даёт параметризацию
$$U=az,\qquad V=b\sqrt{-3}\,z,$$
либо ту же пару с переставленными сторонами. Здесь $a,b$ — положительные целые, а $z$ рассматривается с точностью до умножения на единицу кольца. Площадь прямоугольника равна $ab\sqrt3\,|z|^2$, откуда число треугольников
$$n=4ab|z|^2.$$
В частности, регулярное решёточное замощение возможно только при $4\mid n$. Поведение простых в кольце зависит от остатка по модулю $3$:
- если $p\equiv1\pmod3$, простое распадается на два сопряжённых множителя и даёт вклад $(e+1)^2$;
- если $p\equiv2\pmod3$, простое остаётся простым и даёт вклад $e+1$;
- простое $3$ разветвляется, но в этом подсчёте также даёт вклад $e+1$.
Перестановка двух перпендикулярных сторон добавляет множитель $2$. Поэтому
$$ E(n)= \begin{cases} \displaystyle 2\prod_{p^e\parallel n/4} \begin{cases} (e+1)^2,&p\equiv1\pmod3,\\ e+1,&p\not\equiv1\pmod3, \end{cases} &4\mid n,\\[1.2em] 0,&4\nmid n. \end{cases} $$
Удобно вынести произведение в мультипликативную функцию
$$ H(m)= \prod_{p^e\parallel m} \begin{cases} (e+1)^2,&p\equiv1\pmod3,\\ e+1,&p\not\equiv1\pmod3. \end{cases} $$
Тогда $E(4m)=2H(m)$.
Итоговая формула
Пусть
$$M=\left\lfloor\frac N2\right\rfloor, \qquad L=\left\lfloor\frac N4\right\rfloor.$$
После удаления двух лишних копий каждой регулярной решётки получаем
$$ \boxed{ \begin{aligned} G(N)={}&2D(M)\\ &+4\sum_{k\geq1}\sum_{x\geq1}2^{\omega(x)} \left[ \left\lfloor\frac Mk\right\rfloor- \left\lfloor\sqrt3\,kx\right\rfloor \right]_+\\ &-4\sum_{m=1}^{L}H(m). \end{aligned} } $$
Теперь у трёх строк есть ясный геометрический смысл:
- полосы параллельны стороне прямоугольника;
- полосы наклонены, а $2^{\omega(x)}$ считает примитивные типы обхода;
- вычитаются две лишние копии регулярных треугольных решёток.
Как считать это без новых сюрпризов
Делительная сумма
В $D(X)$ одинаковые частные можно обрабатывать блоками. Если $q=\lfloor X/d\rfloor$, то это значение сохраняется до $r=\lfloor X/q\rfloor$. После учёта всего отрезка $d,\ldots,r$ переходим к $r+1$. Таких блоков порядка $\sqrt X$.
Решета для двух мультипликативных функций
Значения $2^{\omega(x)}$ удобно получать решетом: каждый новый простой делитель удваивает значение у всех своих кратных. Функция $H(m)$ тоже считается решетом по минимальному простому делителю и показателю его степени. Число реально встречающихся пар $(k,x)$ линейно по $N$, поскольку
$$ \sum_{k\geq1}\frac{M}{\sqrt3\,k^2} <\frac{M}{\sqrt3}\cdot\frac{\pi^2}{6}. $$
Ловушка с $\sqrt3$
Здесь я едва не сделал обычный вариант через double.
В обсуждении несколько участников отдельно предупреждают, что такой код проходит маленькие проверки,
но ломается именно на $N=10^9$ из-за округления числа, почти совпавшего с целым.
Один из них потратил на поиск этой причины около восьми часов — это чужой опыт, но очень полезное предупреждение.
Вещественные вычисления здесь вообще не нужны. Для положительного целого $t$
$$ \left\lfloor t\sqrt3\right\rfloor =\left\lfloor\sqrt{3t^2}\right\rfloor. $$
Правая часть вычисляется целочисленным квадратным корнем. Число $3t^2$ не бывает полным квадратом при $t>0$, поэтому пограничной неоднозначности нет. Остаток по модулю $1\,000\,000\,007$ можно брать после каждого сложения и умножения, но границы сумм и целые части нужно считать до модульного сокращения.
Проверки
Я проверял не только готовое $G(N)$, но и две части, из которых оно собирается: сырой полосовой подсчёт $S(N)$ и накопленное число регулярных решёток $E_{\leq N}=\sum_{n\leq N}E(n)$. Получилось:
- для $N=6$: $S=18$, $E_{\leq N}=2$, поэтому $G=14$;
- для $N=100$: $S=8494$, $E_{\leq N}=202$, поэтому $G=8090$;
- для $N=10^5$: $S=22647276698$, $E_{\leq N}=1076248$, поэтому $G=22645124202\equiv645124048\pmod{1\,000\,000\,007}$.
Последняя проверка особенно полезна: она достаточно большая, чтобы поймать ошибки в решете, поправке за регулярные решётки и вычислении $\lfloor\sqrt3\,kx\rfloor$.
Что в итоге оказалось самым красивым
Больше всего мне нравится, что непрерывные деформации здесь не уничтожают дискретность, а проявляют её. Полосы можно двигать, но числа обходов $(a,b)$, число компонент $k$ и выбор ветви $A-B$ непрерывно не меняются. Поэтому классы эквивалентности кодируются целыми числами.
А единственное место, где описание полосами неоднозначно, — идеальная треугольная решётка. Её лишняя симметрия сначала портит ответ тройным подсчётом, а затем неожиданно приводит к целым Эйзенштейна и разбиению простых по модулю $3$. В одной задаче встретились топология тора, евклидова геометрия, делители, решета и алгебраическая теория чисел.
Для требуемого $N=10^9$ остаётся подставить $M=500000000$ и $L=250000000$ в итоговую формулу, выполнить сравнения с $\sqrt3$ через целочисленный корень и взять остаток по заданному модулю. Готовый код и итоговый численный остаток здесь не привожу.
