我用Brainfuck编写了一个光线追踪器

# 用Brainfuck编写光线追踪器 来源：https://epestr.com/blog/writing-a-ray-tracer-in-brainfuck/ 在准备C++系统编程比赛时，我重新学习CMake——毕竟Cargo已把我惯坏——注意到教程中一个有趣的说法： > 通常正确的做法是使用通用编程语言编写解决问题的工具，并在构建过程中通过CMake调用该工具。虽然有人用CMake语言编写过代码生成器、加密签名工具甚至光线追踪器，但这并非推荐做法。CMake语言基础（https://cmake.org/cmake/help/latest/guide/tutorial/CMake%20Language%20Fundamentals.html） 此前我写过光线追踪器，也为GPU重写过，这个说法引起了我的兴趣：究竟哪种语言更适合编写光线追踪器？ 上次重写时，代码很少基于第一性原理，主要依赖相对复杂的API集合。因此我选择了最熟悉的简单语言——Brainfuck，因为简单的语言自然能产生非常简洁的代码库。实际上，Brainfuck代码库通常只有几行。此外，Muller（https://en.wikipedia.org/wiki/Brainfuck#History）在README中的注释激发了我想展示一个反例的欲望。 代码托管在atmTvare6/rayfuck（https://github.com/mTvare6/rayfuck）。 ## 基础知识 Brainfuck是一种极其简单的语言，仅包含8个操作符和一种"数据结构"：一条单向无限磁带，每个单元可存储一个`u8`值。 遇到`>`字符时，指向磁带单元的数据指针右移，`<`则左移。 输入输出通过`,`和`.`管理。前者将输入字节存入数据指针指向的单元，后者将其输出。 除I/O外，唯一能修改值的操作是通过`+`和`-`在数据指针位置进行递增和递减。 其局限性显而易见：没有类似其他机器的多寄存器（n>1）、没有操作多个单元的指令、也没有加法或乘法指令。 使Brainfuck具备图灵完备性的最后要素是其循环结构，使用`[`和`]`实现。遇到`[`时，运行时检查数据指针指向的单元：若为零则跳转至匹配的`]`，否则进入循环。在`]`处，若单元非零则跳回匹配的`[`，否则退出循环。 一个快速练习是编写`cat`程序，尝试仅用5个字符实现，以熟悉环境。 ## 事前规划 在开始前阅读相关资料后，我决定避免查阅任何现有结果或实现细节，尽可能从第一性原理出发编写。为保持范围最小且程序是明显的光线追踪器，我决定让它渲染与[RIW（https://raytracing.github.io/books/RayTracingInOneWeekend.html#metal）]金属部分完全相同的图像。 无论如何，C代码过于复杂，编写C解析器显然超出范围。编写难以维护的C解析器更适合交给Anthropic（https://www.anthropic.com/engineering/building-c-compiler）处理。 我决定每个双精度数[及其他数据类型如布尔值]通过组合单元表示，半数位表示小数部分，另一半表示整数部分，实际上在它们之间设置了固定的二进制点。后来我得知这被称为Q格式。采用更经济的带符号Q8.8格式可提供`1/256`的分辨率，范围约为`[-128, 128)`。但这显然不够，因为场景中地面球体需要`r=1000`才能显得平坦，所以我选择了更昂贵的带符号Q16.16格式。其分辨率为`1/2^16`，范围是`[-2^15, 2^15)`，已足够使用。 我决定将代码转换为类SSA格式[并决定这将是LLM的唯一任务]，将递归代码转为迭代形式，并用匈牙利式记法为函数中定义的变量添加前缀，避免名称查找时的地址冲突。 类似地，解析与代码生成分离是必要的，将复杂度分为两个代码区域，使用中间"DSL"作为IR。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`。最初我计划采用双状态方案： ``` A = (A-B) % 256 B = (B+1) % 256 或 B = (B+A+p) % 256 # p为某质数 ``` 但大多数变体周期性较差。我决定采用更简单的方案： ``` A = (A + B) % 256 B = (B + 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`。 最后，整个变量映射和对应的Brainfuck地址将用字典维护。 ## 实现 理论部分完成后，只剩下清晰简单的实现细节。两个重要原语是`move`和`copy`。 ``` [a, 0] ``` 移动通过持续减小值直到初始数据指针处单元为零，同时等量递增另一单元来实现。 ``` [ # 开始循环 - # 递减 >+ # 右移并递增 < # 返回，此单元用于控制循环 ] ``` 写成单行： ``` [->+<] ``` 复制从数组开始： ``` [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语言的真值倾向。构造的表示中，真值最低位为1，假值所有位为零。`and`、`or`和`not`等布尔运算使用该位表示。 取反使用`-x = ~x + 1`技巧。每个单元`x`取补[通过`255-x`]，然后向最低单元加一并进位。`abs`仅检查最后一个单元的最高位，若设置则执行此取反。 鉴于之前决定按函数限定变量名并使用SSA式代码，函数极其简单。代码生成将函数体记录在其名称下，看到`call func`时内联发出。 有了循环，`if`很简单。 对于`else`，另一个标志从1开始，被第一个主体清除。 `while`同理。 ``` 条件 [ 主体 移至条件并计算 ] ``` ## 最终产物 C版本输出 程序最终`23MB`，大于图像本身[约`0.9MB`]，使其成为相当差的压缩技术选择。 通过粗略计算，我发现其每分钟执行100次光线计算，即每分钟一个像素。鉴于图像为`400x225`，初始估计在我的笔记本上无进一步优化应约62.5天，但我意识到只看到天空，球体周围的弹射会大幅推迟预计时间。因此上图是用C代码生成的渲染近似值。在撰写时生成的`1229`[共90k]像素中，仅`10`个不同，大多差值为一。 可进行一些优化。若地面可更粗糙近似，减少触及单元数是首选。随机向量的归一化步骤也可跳过，尽管这会改变散射分布，因此将不再运行我旨在重现的精确代码。