[Перевод] Пишем трассировщик лучей на Brainfuck

В процессе подготовки к соревнованиям по системному программированию на C++ я начал заново изучать CMake, потому что Cargo сильно меня избаловал. При этом я заметил в туториале интересное утверждение:

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

CMake Language Fundamentals

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

Для последней версии я писал код, почти не связанный с базовыми алгоритмами и в основном зависевший от довольно сложного набора API. Поэтому я выбрал самый простой из известных мне языков — Brainfuck, ведь простота языка очевидным образом приводит к созданию простой кодовой базы. На самом деле, кодовые базы на BF обычно состоят всего из нескольких строк. Кроме того, комментарий Урбана Мюллера в README заставил меня написать контрпример.

Код выложен на Github.

Основы

BF — крайне простой язык, в нём есть всего восемь операций и одна «структура данных»: односторонняя лента ячеек, в каждой из которых можно хранить u8.

Встретив символ >, указатель данных, указывающий на ячейку ленты, перемещается вправо, и наоборот в случае символа <.

Ввод‑вывод реализован через , и .. Первый сохраняет байт ввода в ячейку, на которую установлен указатель данных, а второй выводит его.

Кроме ввода‑вывода, единственные примитивы, позволяющие менять значение — это инкремент (+) и декремент (-) по указателю данных.

Ограничения такой системы очевидны: здесь нет регистров n > 1, как в других машинах, нет команд, работающих с несколькими ячейками, и нет команд сложения и умножения.

Последний ингредиент, необходимый для превращения BF в Тьюринг‑полный язык — это цикл, записываемый при помощи [ и ]. Когда в коде встречается токен [, среда проверяет ячейку по указателю данных: если она равна нулю, исполнение перескакивает через соответствующий ], а в противном случае входит в цикл. На токене ], если ячейка равна нулю, среда исполнения возвращается к соответствующему [, а в противном случае выходит из цикла.

В качестве небольшого упражнения можно написать программу cat, попробовав реализовать её всего в пяти символах, чтобы освоиться с окружением.

Замысел

Перед тем, как приступать к проекту, я решил не искать никаких результатов и подробностей реализации и писать всё по возможности с нуля. Чтобы масштаб проекта оставался минимальным, а программа была очевидным трассировщиком лучей, я выбрал для рендеринга изображение, которое рендерится в разделе Metal RIW.

Код на C всё равно слишком сложный, а создание парсера C явно сильно выходит за рамки масштабов проекта. С написанием неподдерживаемого парсера C лучше справится Anthropic.

Я решил, что все double [и другие типы данных наподобие bool] будут представлены комбинациями ячеек: половина бит будет обозначать дробную часть, а другая половина — целую; таким образом, мы, по сути, ставим между ними фиксированную запятую. Позже я узнал, что это называется Q format. Выбор Q8.8 со знаком дал бы мне разрешение 1/256 и диапазон, приблизительно равный [-128, 128). Но очевидно, что этого было бы недостаточно, потому что сфера, используемая в качестве основания сцены, должна иметь r=1000, чтобы казаться плоской, поэтому я выбрал более затратный формат Q16.16. Он имеет разрешение 1/2^16 и диапазон [-2^15, 2^15), чего вполне достаточно.

Код будет преобразовываться в формат наподобие SSA [я решил, что это единственная задача, которую я отдам LLM]: рекурсивный код будет превращаться в итеративный, а определяемые в функциях переменные будут иметь префиксы в стиле венгерской нотации, чтобы избежать коллизий имён при поиске адреса для имени.

Также мне показалось необходимым отделить парсинг от генерации кода: сложность кода будет разделена на две области, где в качестве промежуточного представления используется DSL. Он содержит простые операции наподобие abs, add, and, call, copy, div, else, end, eq, func, ge, gt, if, int, le, lt, mul, neg, not, or, print2, print3, set, sqrt, sub, text, var и while.

Следующим сложным моментом была замена вызовов функций. Я использовал sqrt, rand и abs. Изначально я планировал использовать для rand решение из двух состояний вот такого вида:

A = (A-B) % 256
B = (B+1) % 256
или
B = (B+A+p) % 256 # для некоторого простого p

но у большинства таких вариантов был плохой период. Я выбрал более простое решение:

A = (5*A + 1) % 256

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

Для sqrt был один очевидный кандидат — формула Герона [о которой мне напомнил тот же туториал по CMake]. Но учитывая деление, было очевидно, что она плохо подойдёт из‑за используемого деления. Многократное вычитание, позволяющее к тому же генерировать меньший размер кода [что предпочтительнее, потому что интерпретатору придётся работать меньше], всё равно было относительно затратным в реализации. Оставшимися кандидатами были ряды Тейлора и «школьный алгоритм» с делением столбиком.

sqrt(1 + x) = 1 + x/2 - x^2/8
y = 2^16*x
sqrt(2^16 + y) / 2^8  = (1 + y/2^17 - y^2/2^35 )
sqrt(y) / 2^8  = (1 + (y - 2^16)/2^17 - (y - 2^16)^2/2^35 )

Создав график в Desmos, я увидел, что выравнивание было плохим ниже закодированного значения 20k, то есть ниже примерно 0,305, а это довольно важная область. Остался только алгоритм с делением в столбик, а он достаточно прост. Если реальное значение равно x, то представленное в системе значение имеет вид:

N = x * 2^16

Чтобы выразить sqrt(x), нужно сделать следующее:

N' = sqrt(x) * 2^16
isqrt(N) = sqrt(x) * 2^8
isqrt(2^16 * N) = sqrt(x) * 2^16 = N'

Использование isqrt здесь оправдано, потому что различие в закодированном результате изменяет декодированный квадратный корень менее чем на 1/2^16, или приблизительно на 0,00001526.

И наконец, вся map переменных и соответствующие ей адреса BF поддерживаются с помощью словаря.

Реализация

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

[a, 0]

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

[ # начало цикла
    - # декремент
    >+ # сдвиг вправо и инкремент
    < # возвращаемся назад, эта ячейка используется для управления циклом
]

В виде одной строки это выглядит так:

[->+<]

copy берёт массив:

[a, 0, 0]

и при помощи кода

[->+>+<<]

превращает его в следующий:

[0, a, a]

Затем при необходимости последнюю a можно переместить внутрь.

Учитывая, что для сложения, а потом и деления потребуется многократное применение временных значений, которые должны находиться рядом со значением, чтобы избежать чрезмерно далёкого перемещения указателя данных, у каждого значения в map есть поблизости свои слоты временных значений. Также у него есть ячейки переноса и другое полезное временное пространство, позволяющее ограничить область работы с копированием.

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

[a, b, a->0, b->0, a + ... + a]

Здесь a исполняется каждый раз, а для b выполняется декремент при нулевом a, что создаёт b копий a.

В нашем случае значения хранятся в четырёх ячейках, поэтому каждая ячейка одного по очереди обрабатывается с каждой ячейкой другого. Умножение ячеек в i м j прибавляется в виде i+j к результату, состоящему из восьми ячеек.

[a0, a1, a2, a3] * [b0, b1, b2, b3]

result[i+j] += a_i*b_j

Так как оба входных значения уже имеют 2^16 в своём представлении, две младшие ячейки при копировании результата отбрасываются.

N_1 = x_1 * 2^16
N_2 = x_2 * 2^16

( N_1 * N_2 ) / 2^16 = N = (x_1 * x_2) * 2^16

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

R = R * 10 + A[next] # школьный алгоритм
R = R * 256 + A[next] # наш алгоритм

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

while R >= D:
    R -= D
    result += 1

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

Аналогично тому, как представление сдвигается вправо, увеличивая себя во время умножения, так и деление теряет информацию из‑за сдвига влево, и перед делением необходимо выполнение сдвига представления.

N_1 = x_1 * 2^16
N_2 = x_2 * 2^16

(N_1 / N_2) * 2^16 = (x_1 / x_2) * 2^16 # биты уже потеряны
(N_1 * 2^16) / N_2 = (x_1 / x_2) * 2^16

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

В сравнениях используется одна и та же маленькая операция. Выполняется одновременный декремент двух ячеек, пока одна из них не станет равна нулю, и этот процесс продолжается, пока существует разность или пока временные копии полностью не обнулятся. Кроме того, происходит небольшая обработка, связанная с прибавлением 2^7 к старшей ячейке, поскольку в противном случае отрицательные числа в байтовом виде будут иметь большее значение.

00 ... 7f  -> положительная половина
80 ... ff  -> отрицательная половина

после прибавления 128 и циклического переноса

80 ... ff  -> положительная половина
00 ... 7f  -> отрицательная половина

Проверка логических значений заключается в чтении ячеек и установке единицы в выходном значении, если хотя бы одна из них была ненулевой, чтобы учесть свойственное C восприятие ненулевых значений как истинных. В полученном представлении для значения true устанавливается младший бит, а для false все биты равны нулю. Логические операции, такие как and, or и not, выполняются с использованием этого битового представления.

В отрицании применяется трюк -x = ~x + 1. Каждая ячейка x инвертируется [путём вычисления 255-x], после чего к младшей ячейке прибавляется единица с переносом. abs проверяет только старший бит последней ячейки и выполняет эту инверсию, если он установлен.

Учитывая принятое ранее решение ограничивать область видимости имён переменных рамками функции и использовать код в стиле SSA, функции реализуются чрезвычайно просто. Генератор кода сохраняет тело функции под её именем и вставляет его непосредственно в код при встрече конструкции call func.

Благодаря циклам if тоже было реализовать легко.

[- body ]

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

Аналогично с while.

condition
[ 
    body
    move to condition and calculate
]

Артефакт

Итоговая программа имеет размер 23 МБ, что больше самого изображения [примерно 0,9 МБ], поэтому в качестве техники сжатия она не подходит.

Я приблизительно подсчитал, что в минуту выполняется 100 вычислений лучей, по одному пикселю за минуту. Учитывая размер изображения 400x225, изначально я предполагал, что без оптимизаций на моём ноутбуке процесс займёт примерно 62,5, но понял, что увидел бы только небо: отражения лучей между сферами ещё сильнее бы увеличили сроки. Поэтому показанное выше изображение — написанная на языке C аппроксимация того, что могло отрендериться. Из 1229 пикселей [всего их 90 тысяч], сгенерированных на момент написания статьи, от аппроксимации отличались только 10, чаще всего на единицу.

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

Дополнение

На Reddit мне задали вопрос о том, как можно улучшить код при помощи примитивов fork/join. Эта задача показалась мне интересной, поэтому я занялся совершенствованием используемого JIT‑интерпретатора, что привело к огромному росту производительности. Итоговый рендер из‑за погрешностей точности немного походит на картину ван Гога.

Если эта публикация вас вдохновила и вы хотите поддержать автора — не стесняйтесь нажать на кнопку