用 Brainfuck 写光线追踪器

为了准备一场 C++ 系统编程比赛，我开始重新学习 CMake——这期间 Cargo 把我惯得太厉害了——然后在教程里看到一句很有意思的说法。

很多时候，正确的做法是用一门通用编程语言写一个工具来解决问题，再教 CMake 如何在构建过程中调用这个工具。代码生成、密码学签名工具，甚至光线追踪器，都有人用 CMake Language 写过，但这并不是推荐做法。CMake Language Fundamentals

我之前写过一个光线追踪器，后来又为 GPU 重写了一遍，所以这句话格外扎眼，让我开始琢磨：还有什么语言比它更适合写光线追踪器？上一次重写时，我写的代码几乎没有从第一性原理出发，基本依赖一套相对复杂得多的 API。于是我挑了自己知道的最简单的语言——BF，因为简单的语言显然会带来非常简单的代码库。事实上，BF 代码库往往只有几行。而且 Muller 在 README 里的那句评论让我想举一个反例。代码在 mTvare6/rayfuck。

入门

BF 是一门极其简单的语言，只有 8 个操作和一个“数据结构”：一条单侧无限的单元带，每个单元可以存一个 u8。看到字符 > 时，指向带上某个单元的数据指针向右移动，看到 I/O is managed through 时则向左移动，and。. The first stores the input byte where the data pointer points and the latter prints it out.The only primitives other than I/O which allow changing a value are increment and decrement at the data pointer, through +and-.The limitations should be obvious: there are no n > 1 registers as other machines tend to have, no instruction operating on more than one cell, and no instructions for addition or multiplication.The last ingredient required to make BF Turing-complete is its loop, written using [and]. When the token [is met, the runtime checks the cell at the data pointer: if it is zero, execution jumps past the matching], otherwise it enters the loop. At ], it returns to the matching [if the cell is non-zero and exits the loop otherwise.A quick exercise would be writing acat` 程序，试着只用 5 个字符写一个，感受一下这个环境。

预谋

在动手之前读过相关资料后，我决定不去查任何结果或实现细节，尽量从第一性原理出发把东西写下来。为了让范围尽可能小，同时让程序明显是个光线追踪器，我决定让它渲染 RIW 的 Metal 一节所渲染的那张一模一样的图。但那段 C 代码还是太复杂了，写一个 C 解析器显然超出了范围。写一个没人维护的 C 解析器，这种事还是交给 Anthropic 更合适。我决定每个 double（以及其他类似 bool 的数据类型）都用多个单元组合来表示，一半比特表示小数部分，另一半表示整数部分，相当于在两者之间放一个固定的二进制小数点。后来我才知道这叫 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
or
B = (B+A+p) % 256 # for some prime p

但这些变体的周期大多很差。我最终选了更简单的

A = (5*A + 1) % 256

，因为它能保证在完整的 256 个值之后才会重复，对这个用途（也就是超采样抗锯齿）来说已经够用了。

sqrt有一个显而易见的候选方案：Heron 公式（我正好在同一篇 CMake 教程里复习了一遍）。但很明显它效果会很差，因为涉及除法。反复减法虽然生成的代码更小（这更好，解释器要移动的指令更少），但执行起来仍然相当昂贵。 其余候选是 Taylor 级数和“学校方法”，后者要用长除法。

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，解码后的平方根变化不到1/2^16，约等于0.00001526。

最后，整个变量映射以及对应的 BF 地址用一个字典来维护。

实现

理论部分搭好之后，剩下的实现细节都很简单。两个重要的原语是move和copy。

[a, 0]

Move 的做法是不断减小一个值，直到初始数据指针所在的单元变为零，同时每次把另一个单元等量增加。

[ # start loop
 - # decrement
 >+ # move right and increment
 ++>+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，拷回结果时最低的 2 个单元会被丢弃。

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] # school
R = R * 256 + A[next] # here

然后从这个余数中反复减去除数，每减一次就给结果单元加一。当余数变为负数时，撤销这一步，移到下一个单元。

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 # bits already lost
(N_1 * 2^16) / N_2 = (x_1 / x_2) * 2^16

在这些运算中，首先去掉符号，如果只有一个输入为负，则将结果设为负。

比较运算共用同一个更小的操作。两个单元一起递减，直到至少一个变为零，然后继续直到出现差异或临时副本完全归零。还有一些小处理涉及向最高单元加上 2^7，否则负数在直接按字节看待时技术上具有更高的值。

00 ... 7f -> positive half
80 ... ff -> negative half

加上 128 并回绕之后

80 ... ff -> positive half
00 ... 7f -> negative half

布尔检查涉及读取各单元，如果其中任何一个非零则将输出设为一，以照顾 C 的真值倾向。构造出的表示中，最低位为一表示真，所有位为零表示假。诸如 and、or 和 not 这样的布尔运算就是基于该位表示工作的。

取负使用 -x = ~x + 1 技巧。每个单元 x 被取反[通过 255-x]，然后向最低单元加一，并进位。abs 只检查最后一个单元的最高位，如果该位被设置就执行这个取负。

鉴于之前决定按函数限定变量名并使用 SSA 风格的代码，函数的处理极其直接。代码生成器在函数名之下记录函数体，并在看到 call func 时将其内联发出。

有了循环之后，if 就很直接了。

[- body ]

对于 else，另一个标志从一始，由第一个循环体清除。

而 while 同样如此。

condition
[ 
 body
 move to condition and calculate
]

制品

这个程序最终有 23MB，比图像本身（约 0.9MB）还大，作为压缩手段实在不怎么样。

粗略估算下来，它每分钟做 100 次光线计算，也就是每分钟一个像素。图像是 400x225，照此推算，不做任何优化，在我的笔记本上大概要 62.5 天，但我意识到自己只看到了天空部分，光线在球体之间的来回反弹会把预计时间拖长很多。所以上面那张图是用 C 代码渲染出的近似结果。写作时生成的 1229 个像素（总共 9 万）中，只有 10 个不同，且大多只差 1。

还有优化的余地。如果地面可以近似得更粗糙，那么牺牲精度、减少访问的格子数是首选。随机向量的归一化步骤也可以跳过，不过这会改变散射分布，就不再是我这里想要复现的那段代码了。

更新

Reddit 帖子里有人问，用 fork/join 原语能怎么改进。这个问题挺有意思，我就一头扎进去优化了所用的 JIT 解释器（借助了一些早先的工作），性能因此大幅提升。实际渲染出来有点像梵高的画，大概是精度误差导致的。

← 所有文章