Я хочу сыграть с вами в одну игру...

— Дж. Конвей, про игру "Жизнь"

Эту статью я хотел написать ещё полтора года назад, но почему-то замешкался и стёр черновик. Но вот наконец дошли руки.

Суть игры

Вы задаёте список неотрицательных рациональных дробей (программу) и некоторое натуральное число (состояние).
Ищете в программе первую дробь, такую, что произведение её на состояние остаётся натуральным. Это значение становится новым состоянием.
Повторяете процедуру.
Если ни одна из дробей не подошла, программа останавливается.

Этот список дробей является программой на эзотерическом языке Фрактран, о котором на хабре была статья.

Название fractran происходит от fraction tran...что-нибудь (translation, transformation, и т.п.)
Поэтому на русский его можно перевести, как "дробедробилка".

Ну и раз в этой статье упомянута гипотеза Коллаца, то понятно, какой именно хелловорлд мы будем писать на фрактране.

Хелловорлд: последовательность Коллаца

Напомню: последовательность Коллаца задаётся рекуррентным выражением:

Примем за постулат, что последовательность Коллаца всегда сходится к циклу 1-4-2-1. (Это не доказано, но проверено для очень большого диапазона).

Поэтому нас будут интересовать не сами элементы последовательности, а их количество до 1. Обозначим эту функцию как L(x), а заодно как K(x)

L(1) = 1 (сразу встретили)
L(2) = 2 (x1=2, x2=1, всё)
L(3) = || [3, 10, 5, 16, 8, 4, 2, 1] || = 8
L(7) = || [7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1] || = 17
L(27) = || [27, 82, ..., 9232, ..., 1] || = 112

K(3) = 16, K(7) = 52, K(27) = 9232.

Вернёмся к фрактрану

Очевидно, что, поскольку мы работаем с умножением дробей, то удобно представлять числа как разложения на простые множители (и делители):

, то есть, , где все показатели - целые (положительные - в числителе, отрицательные - в знаменателе). И дальше иметь дело только с этими векторами, причём основания степеней для нас несущественны.

Пример программы, которая делит число на 3 с остатком: из [x, 0, 0] получает [0, x div 3, x mod 3]:

1. [-3, +1, 0] — если можем вычесть из делимого 3, то вычитаем, и добавляем к частному единичку

2. [-1, 0, +1] — очевидно, что делимое уже меньше 3, то есть, это остаток; перекладываем в соответствующую ячейку по единичке

А вот пример программы, которая бесконечно прыгает между тремя состояниями [1,0,0], [0,1,0], [0,0,1]:

1. [-1, +1, 0]

2. [0, -1, +1]

3. [+1, 0, -1]

Возвращаемся к хелловорлду

Начальное состояние нашей программы, естественно, должно содержать входное число x, как показатель степени простого числа, а окончательное - вычисленную функцию, также как показатель. Понятно, что с другим основанием, иначе как мы (и как интерпретатор) отличит вход от выхода.

Теперь мы оговорили все исходные условия задачи, и я предлагаю самостоятельно подумать над тем, как написать программу на фрактране, которая из первого состояния приводит к последнему.

Ну а дальше, конечно же, сплошные спойлеры и мысли вслух.

Как устроена фрактран-программа

Во-первых, удобно представить программу, как граф некоторого автомата, а состояние машины как четвёрку:
<вход, выход, переменные, вершина>.
В принципе, вход и выход - это тоже переменные, просто подчеркнём, что они лежат на самом видном месте :)

Во-вторых, каждой вершине соответствует набор исходящих из неё рёбер:
<вершина откуда, вершина куда, условия и изменения переменных>.
Находясь на данной вершине, интерпретатор за один такт выбирает первое подходящее ребро и переходит по нему на другую вершину.

Идентифицировать вершины можно, как минимум, двумя способами. Первый и самый напрашивающийся - пронумеровать. Другой - сопоставить каждой вершине флажок, а всему семейству вершин - вектор, в котором только один флажок взведён. Давайте посмотрим, к чему приведут оба способа.

Нумерованные вершины

Фрактран-код для перехода из вершины X в вершину Y будет выглядеть так. Для начала, надо убедиться, что ребро исходит из X, и мы в состоянии X. Сделать это можно с помощью вычитания. [(тут переменные), -X, ...]. Это значит, что в программе все рёбра расположены строго по убыванию номеров вершин, ведь иначе у нас были бы ложные срабатывания для состояний X+dX. Это также значит, что рёбра каждой вершины сгруппированы вместе.

При этом сама по себе нумерация вершин в графе может быть произвольной. Перенумеруем по-другому - просто переставим местами фрагменты программы (и разумеется, поправим номера в рёбрах).

Единственная финальная вершина должна иметь минимальный номер (0). А из всех остальных вершин обязательно должны существовать исходящие рёбра на все случаи жизни. В противном случае, если интерпретатор не найдёт подходящее ребро из вершины X, то он перескочит к вершине X-dx, где случайным образом найдёт ребро, да ещё и не обнулит, а только уменьшит текущий номер.

Итак, переход "откуда" - понятно, как сделать, вычитанием. Но вот как сделать переход "куда"? Прибавлением? Но мы же не можем написать [..., -X+Y, ...], эти два числа сократятся. И вот тут мы прибегнем к технике пинг-понга.

Пусть вершина идентифицируется не одним, а двумя числами: (X,0) и (0,X). То есть, рядом с рабочей вершиной автомата лежит её зеркало.

Разрежем каждое ребро пополам: X->Y превратим в (X,0) -> (0,Y) -> (Y,0). Причём первое ребро проверяет какие-то условия переменных (об этом позже), а второе делает безусловный переход. В программе это будут [..., -X, +Y] и [..., +Y, -Y]. Единственная проблема - в том, что теперь каждый полный переход по ребру исходного автомата занимает не один, а два такта.

Можем поступить иначе: отзеркалим наш граф полностью. И из "левых" вершин с номерами X (с идентификатором (X,0)) протянем рёбра в "правые" (0,Y), а из правых - точно такие же в левые. И тогда у нас будут [..., -X, +Y] и [..., +Y, -X]. Теперь мы тратим на переход один такт.

Вершины с флажками

Номеру X тут соответствует вектор вида

Ребро X->Y - это сбросить флажок номер X и взвести номер Y.
То есть,

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

Теперь из каждой вершины есть неявное ребро в финальную вершину с нулевым вектором флажков. И, как и с нумерованными, - код этой вершины должен быть в самом конце программы. А вот остальные фрагменты кода можно даже не группировать, а перетасовывать вообще как угодно, - сохраняя лишь относительный порядок рёбер, исходящих из общей вершины.

Единственная незадача - это петли, то есть, рёбра X->X. Как и с номерами, мы не можем одновременно проверить (сбросить) и взвести флажок, -1+1 = 0.

Если у какой-то вершины есть петля, то снова сделаем зеркальную. То есть, вместо единственной X у нас будут вершины X и X+1. При этом даже не обязательно зеркалить весь граф!

Для двухтактного перехода - петля будет кодироваться как [..., -1, +1, ...] и [..., +1, -1, ...].

Для однотактного - сделаем полноценное зеркало вершины. Продублируем все исходящие рёбра, для каждого X->Y добавим (X+1)->Y. Входящие рёбра можно произвольным образом повесить хоть на X, хоть на X+1, - поскольку с логической точки зрения эти вершины неразличимы.

Ветвления

С переходами между вершинами мы разобрались. Но что же содержательного могут делать рёбра?

Всё, что умеет интерпретатор фрактрана, - это проверить, что некоторое подмножество переменных содержит значения, не менее, чем заданные константы, - и немедленно выпить вычесть! И тут же к другому подмножеству переменных что-нибудь прибавить.

Как, например, ранее в программе деления на три с остатком. У нас там единственная вершина и петли в неё, поэтому не потребовалось городить нумерацию.

Три переменных: делимое, частное, остаток.

Петля номер один: "если делимое >= 3, то +1 к частному, и какая удача, -3 к делимому".
Петля номер два: "если делимое >= 1, то +1 к остатку, и опять же, -1 к делимому".

Если бы нам захотелось сделать неразрушающее сравнение, или просто вычесть и тут же прибавить, то, как и с номерами вершин, мы просто разрежем ребро и введём промежуточную вершину: первое ребро честно вычитает, второе прибавляет константы.

Заметим, что рёбра одной вершины образуют "лесенку решений": проверяем первое условие (и идём по первому ребру), иначе проверяем второе, иначе третье...

Арифметика на фрактране

Пока мы напрочь забыли про Коллаца, просто научимся разным трюкам.

Обнуление: (a) -> (0).
Подпрограмма элементарная: одна вершина и петля [-1]. Просто будет вычитать единички, пока не дойдёт до нуля.
Можем ускорить исполнение, сделаем несколько петель: [-1000], [-100], [-10], [-1]. Естественно, расположим их по убыванию, иначе петля [-1] подойдёт наравне с другими, но подойдёт раньше других, и мы ничего не выиграем.

Перенос: (a, 0) -> (0, a).
Тоже из одной вершины и петли: [-1, +1]. В одном месте убыло, в другом прибыло.
Более общо, это инкремент с обнулением: (a, b) -> (0, a+b).
И тоже можем ускорить: [-1000, +1000] и т.д.

Перенос-дублирование: (a, 0, 0, 0) -> (0, a, a, a).
Или такой же инкремент: (a, b, c, d) -> (0, a+b, a+c, a+d).
Оказывается, фрактран умеет в параллелизм: [-1, +1, +1, +1].
Ускорители - [-1000, +1000, +1000, +1000] и т.д.

Перенос с масштабированием: (a, 0, 0, 0) -> (0, ka, ma, pa).
Или (a, b, c, d) -> (0, ka+b, ma+c, pa+d).
Тоже очевидно: [-1, +k, +m, +p]. Ускорители - [-1000, +1000k, +1000m, +1000p] и т.д.

Деление на константу: (a, 0) -> (a mod k, a div k).
Это уже на примере k=3 было показано выше. [-k, +1].
Ускорители - [-1000k, +1000] и т.д.

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

Копирование (или инкремент): (a, b) -> (a, a+b).
Тут уже понадобится два состояния и промежуточная переменная.
START(a, b, 0) -> TEMP(0, a+b, a) -> FINISH(a, a+b, 0)
START(a, b, 0) -> TEMP(0, b, a) -> FINISH(a, a+b, 0)
Заметим, что мы можем объединить второе и третье состояния (то есть, там будет петля), но не должны возвращаться в первое. Иначе, как только там появится ненулевое значение в переменной a, мы снова бросимся его складывать.

Вычитание: (a, b) -> (a-b, 0).
Если у нас есть гарантии, что a>=b, то достаточно петли [-1,-1].
Фактически, это (a, b) -> (min(0, a-b), min(0, b-a)).

Отрицательные числа, увы, придётся кодировать как-нибудь изобретательно. Ведь состояния - это натуральные числа, то есть, показатели там неотрицательные. Так что - либо в виде пары (min(0,x), max(0,x)), либо в виде вынесенного знака (abs(x), bool(x<0)).
Ну, не буду сейчас голову забивать.

Умножение двух чисел: (a, b, 0) -> (0, 0, a*b).
Тут уже будет полноценная программа и одна промежуточная переменная.
Нам понадобится несколько состояний для организации двух вложенных циклов.

1. LOOP_A (a, b, m, 0) ==> в перспективе FINISH(0, b, m+a*b, 0)
   [-1, 0, 0, 0] -> LOOP_B
   [0, 0, 0, 0] -> FINISH

2. LOOP_B (a, b, m, 0) ==> LOOP_C(a, 0, m+b, b)
   [0, -1, +1, +1] -> LOOP_B (помним про трюки с зеркалами)
   [0, 0, 0, 0] -> LOOP_C

3. LOOP_C (a, 0, m, b) ==> LOOP_A(a, b, m, 0)
   [0, +1, 0, -1] -> LOOP_C
   [0, 0, 0, 0] -> LOOP_A

4. FINISH (0, b, m, 0) ==> (0, 0, m, 0)
   [0, -1, 0, 0] -> FINISH

И снова к хелловорлду. Решение "на пальцах".

Итак, у нас есть неотрицательное число x. Будет число n = L(x). И потребуется ещё как минимум одна переменная y, потому что нам нужно проверять x на чётность, делить пополам и всё такое. На самом деле, ровно одна переменная.
Так что состояние описывается четвёркой VERTEX(x, n, y)

Стартовое состояние, по условиям задачи, одновременно является и финишным. У него идентификатор вершины (номер или флажки) нулевой. Это значит, что из финишной вершины есть рёбра куда-то в тело программы при ненулевом x, в противном случае мы остановимся.
Кстати, замечу, что это ограничивает класс программ на фрактране - мы не можем стартовать с полностью нулевыми входными данными.

Итак, распишем наш алгоритм более подробно.

• START(x, n, y=0)

  • если x = 0, то всё, конец. n = L(x)

  • если x = 1, то переходим к START(0, n+1, 0)

  • если x >= 2, то переходим к проверке делимости DIV(x, n, 2)

• DIV(x, n, 0)

  • находим частное и остаток от деления пополам: (x mod 2, n, x div 2)

  • если остаток x = 0, то идём по чётной ветке EVEN(0, n, y)

  • если остаток x = 1, то идём по нечётной ODD(0, n, y)

• EVEN(0, n, y)

  • выполняем x' = x/2 = y

  • переходим к START(x', n+1, 0)

• ODD(0, n, y)

  • выполняем x' = 3x+1 = 6y+4

  • переходим к START(x', n+1, 0)

Оказывается, всё очень просто! Нужно только учесть, что наши проверки разрушающие, а порядок ветвлений имеет значение. Перепишем тщательнее. И отследим значения. Далее x и n - это константы в пределах одной итерации.

• START(x, n, 0) где n - количество уже сделанных итераций

  • если x>=2, перейдём к DIV(x-2, n, 0)
    вектор изменений [-2, 0, 0]

  • если x=1, перейдём к START(0, n, 0)

  • иначе стоп

• DIV(x-2k-2, n, k) где k бежит от 0
  обратите внимание, уходя из START, мы вычли 2, поэтому в начале забега там x-2, а не x

  • если (x-2k-2) >= 2,
    петля DIV(x-2k-4, n, k+1) - сохраняем наш инвариант
    вектор изменений [-2, 0, +1]

  • если (x-2k-2) = 1,
    то есть, 2k = (x-3), k = (x-3)/2,
    привычное нам y = x div 2 = (x-1)/2 = k+1
    перейдём к ODD(0, n, y)
    вектор изменений [-1, 0, +1]

  • иначе (x-2k-2) = 0,
    то есть, 2k = (x-2),
    k = (x-2)/2,
    y = x div 2 = x/2 = k+1,
    перейдём к EVEN(0, n, y)
    вектор изменений [0, 0, +1]

• ODD(6k, n, y-k) где k бежит от 0
  обратите внимание, в инварианте у нас загадочное 6 - потому что нам надо будет получить x'=6y+1

  • если (y-k) >= 1,
    петля ODD(6k+6, n, y-k-1)
    вектор изменений [+6, 0, -1]

  • иначе (y-k) = 0,
    k = y = (x-1)/2,
    x' = 3x+1 = 6k+4
    перейдём к START(x', n+1, 0)
    вектор изменений [+4, +1, 0]

• EVEN(k, n, y-k) где k бежит от 0
  тут всё просто

  • если (y-k) >= 1,
    петля EVEN(k+1, n, y-k-1)
    вектор изменений [+1, 0, -1]

  • иначе (y-k) = 0,
    k = y
    x' = x/2 = y = k
    перейдём к START(x', n+1, 0)
    вектор изменений [0, +1, 0]

Итоговое решение

Как мы увидели, все рабочие вершины содержат петли. Поэтому зеркалить нужно их все, поэтому можем спокойно выбрать идентификацию вершин через нумерацию. (Через пары чисел).

START = 0, это мы уже обсудили. Раздадим остальные номера. Например, так: ODD = 1, EVEN = 2, DIV = 3. Это значит, что в программе сперва будут инструкции для DIV, затем EVEN, затем ODD, и наконец, START.

Ну, поехали. Вектор состояния у нас теперь из 5 элементов: [x, n, y, v1, v2].

Можно ли написать программу короче, чем эта? Подозреваю, что нет, но вдруг у вас получится.

Можно ли написать программу быстрее? О, да! Ведь здесь вся арифметика делается на +1/-1. А, как мы выяснили, K(7) = 27, K(27) = 9232. И в диапазоне до 500 ещё встретятся такие числа, как

Поэтому прибегнем к ускорителям (для краткости я не буду писать зеркальные рёбра):

Можно пойти дальше: захардкодить часто встречающиеся значения, чтобы сразу находить из

Выглядеть это будет примерно так.

Например, для x=9232. Нам понадобятся две вершины, на этот раз не зеркальные, но смежные. Назовём их A9232 (с номером (0,9232)) и B9232 (с номером (9232,0)).

И, чтобы показать принцип, добавим сюда, например, 17 и 13.
Для чётных и нечётных чисел вершины будут устроены по-разному. Связано это с тем, что уходить из них мы будем не только обратно на START, но и на DIV, а он ожидает чётный декремент.

• START(x, n, 0):

  • если x >= 9232
    перейдём к A9232(x-9232, n, 0)
    (мы помним, что проверка у нас - на неравенство и разрушающая)

  • если x >= 17
    перейдём к A17(x-17, n, 0)

  • если x >= 13
    перейдём к A13(x-13, n, 0)

  • и т.п. в порядке убывания особых чисел

  • далее - как обычно

• A9232(x-9232, n, 0):

  • если (x-9232) >= 2, то есть, x >= 9234
    это какой-то общий случай
    перейдём к DIV(x-9232-2, n, k),
    где k находится из (x-9232-2) = x-2k-2, то есть, k=9232/2 = 4616

  • если (x-9232) = 1, то есть, x = 9233
    то почти сразу найдём x', - но для этого надо сделать пинг-понг:
    перейдём к B9232(0, n, 1)

  • иначе x-9232 = 0,
    перейдём к B9232(0, n, 0)

• B9232(0, n, b), где b - флажок: мы имеем дело с 9232 или 9233

  • если b >= 1, то START(9233*3+1 = 27700, n+1, 0)

  • иначе START(9232/2 = 4616, n+1, 0)

• A17(x-17, n, 0):

  • если (x-17) >= 1, то есть, x >= 18
    это какой-то общий случай
    перейдём к DIV(x-18, n, k)
    где (x-18) = x-2k-2, то есть, k=8

  • иначе B17(0, n, 0)

• B17(0, n, 0):

  • безусловный переход на START(17*3+1 = 52, n+1, 0)

Пара вершин нам понадобилась по двум причинам. Во-первых, мы обратно изменяем x, а за один такт невозможно и декрементировать (проверить), и инкрементировать его. Во-вторых, для чётных чисел приходится различать ветви x, x+1, x+2+... - и следовательно, у нас там логически вообще три вершины (или два ребра из второй вершины, и аргумент-флажок).

Что интересно, тем самым мы ускорили не только обработку указанных чисел, но и частично - деление больших чисел. Если на вход подать 10000, то мы сразу вычтем 9234 и перейдём к DIV(766, n, 4616).

Интерпретатор фрактрана

Если мы сразу договоримся работать с разложениями, то, например, код на питоне займёт всего пару десятков строк

# вещественное число представлено вектором с показателями
# простых множителей-делителей
# 0 - не участвует
# +n - множитель
# -n - делитель

# целое число - строго неотрицательные показатели

def is_int(qq):
    '''
    проверка на то, что qq - разложение целого числа
    '''
    return all(q >= 0 for q in qq)

def product(aa, bb):
    '''
    произведение двух разложений - это просто поэлементная сумма
    '''
    assert len(aa) == len(bb)
    return [a+b for (a, b) in zip(aa, bb)]

def step_v(state, program):
    '''
    один такт программы:
    - находим первую подходящую дробь
    - возвращаем новое состояние
    - если ничего не нашли, возвращаем None
    '''
    for case in program:
        new_state = product(state, case)
        if is_int(new_state):
            return new_state
    return None

def run_v(state, program):
    '''
    генератор состояний программы
    '''
    while True:
        state = step_v(state, program)
        if not state:
            break
        yield state

def run_v_until_stop(state, program):
    '''
    забег до финального состояния
    '''
    while True:
        result = step_v(state, program)
        if not result:
            return state
        state = result

Если с дробями - чуть-чуть меньше.

from fractions import Fraction
from typing import List, Optional, Generator

def step_f(state: int, program: List[Fraction]) -> Optional[int]:
    for p in program:
        s = state * p
        if s.denominator == 1:
            return s.numerator
    return None

def run_f(state: int, program: List[Fraction]) -> Generator[int]:
    while True:
        state = step_f(state, program)
        if state is None:
            return
        yield state

def run_f_until_stop(state: int, program: List[Fraction]) -> int:
    while True:
        result = step_f(state, program)
        if result is None:
            return state
        state = result

Разумеется, лучше работать с разложениями, потому что, например, 11^9232 - это ну очень большое число.

Что дальше?

Вычислять длину последовательности L(x) мы научились.

Попробуйте теперь написать программу, вычисляющую максимальное значение K(x).

Ну, не буду лишать вас удовольствия сделать это самостоятельно.

Исходники

https://github.com/nickolaym/fractran/blob/main/fractran-kollatz.ipynb

Юпитер-ноутбук с работающим кодом на питоне.