2.2 Möller–Trumbore 与数值稳定性
站内 /graphics-raytracing-hybrid/
那篇是实时混合光追，求交在硬件里做，开发者看不见 det 长什么样。本卷做离线光追的完整链路，用 Q16.16 定点数——一旦把浮点换成定点，Möller–Trumbore 的除法链就成了必须先算清的账。本节量三件事：除法链的位宽、det → 0 时相对误差的增长、以及原语为什么不做 epsilon 剔除。
为什么要算这条账
Q16.16 里一次乘法要右移 16 位才回到 Q16.16，移位之前的原始乘积决定了能不能用 32 位或 64 位承载。JavaScript 的 Number 做位运算时会被截成 32 位有符号整数，>> 越过 32 位会静默给出错值——不是抛异常，是安静地算错。所以卷三的 TS 侧全程 BigInt。这条账不是理论顾虑：本节与附录 A.2 用同一套抽样口径（320×240 每 8 像素一条主光线，1200 条 × 896 个三角形）把它量出来。
完整脚本
#!/usr/bin/env python3
# file: lantern/sim/cg_trace_02_2_stability.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_02_2_stability.py
"""2.2 Möller–Trumbore 的数值稳定性：det 趋近 0 的相对误差与除法链量级。
1. 量级链：与附录 A.2 同口径——在 320×240 网格上每 8 像素抽一条主光线（1200 条），
对全部 896 个三角形记录 `fmul` 移位前的原始乘积、`fdiv` 的分子、行列式 `det` 的
位宽。原始乘积摸到 2^53（离 Number 的精确整数上界只剩 1.57 倍余量），除法分子是
2^31 的 255 倍——两条独立理由都指向 TS 侧全程 BigInt。
2. det 扫描：命中点固定在三角形内部，方向逐渐共面使 `det -> 0`；用
`fractions.Fraction` 做精确有理数参照，量出定点 (t, u, v) 的相对误差。
3. 为什么不做 epsilon 剔除：统计真实命中光线的 |det| 分布。这些小 det 的命中在
几何上是真的，剔除它们等于凭空丢几何。
全整数/有理数，无浮点。产物：out/cg_trace_02_2_stability.json
"""
import json
import os
import sys
from fractions import Fraction as F
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cglib # noqa: E402
import cgtrace # noqa: E402
ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
OUT = os.path.join(ROOT, "out")
ONE = cglib.ONE
LIMIT_NUMBER, LIMIT_SHIFT = 1 << 53, 1 << 31
def sub(x, y):
return tuple(x[i] - y[i] for i in range(3))
def cross(x, y):
return (x[1] * y[2] - x[2] * y[1], x[2] * y[0] - x[0] * y[2], x[0] * y[1] - x[1] * y[0])
def dot(x, y):
return x[0] * y[0] + x[1] * y[1] + x[2] * y[2]
def world(p):
"""Q16.16 整数三元组 -> 精确有理数世界坐标。"""
return (F(p[0], ONE), F(p[1], ONE), F(p[2], ONE))
def mt_exact(orig, dirv, v0, v1, v2):
"""Möller–Trumbore 的精确有理数版：输入 Q16.16 整数，输出 Fraction 的 (t, u, v)。"""
o, d = world(orig), world(dirv)
a, b, c = world(v0), world(v1), world(v2)
e1, e2 = sub(b, a), sub(c, a)
p = cross(d, e2)
det = dot(e1, p)
if det == 0:
return None
tv = sub(o, a)
q = cross(tv, e1)
return dot(e2, q) / det, dot(tv, p) / det, dot(d, q) / det
def each_ray(spec):
"""逐像素产出 (px, py, dirv)，供两次全场扫描共用。"""
cam = spec["camera"]
eye = tuple(cam["eyeQ16"])
view = cglib.mat_look_at(eye, tuple(cam["targetQ16"]), tuple(cam["upQ16"]))
right, up, back = ((view[0], view[1], view[2]), (view[4], view[5], view[6]),
(view[8], view[9], view[10]))
w, h = spec["width"], spec["height"]
tan_half = cam["tanHalfFovYQ16"]
tan_x = cglib.fmul(tan_half, cglib.fdiv(cglib.f_from_int(w), cglib.f_from_int(h)))
for py in range(h):
sy = cglib.fdiv(cglib.f_from_int(2 * (h // 2 - py) - 1), cglib.f_from_int(h))
ay = cglib.fmul(sy, tan_half)
for px in range(w):
sx = cglib.fdiv(cglib.f_from_int(2 * (px - w // 2) + 1), cglib.f_from_int(w))
ax = cglib.fmul(sx, tan_x)
yield px, py, cglib.v_norm((
cglib.fmul(ax, right[0]) + cglib.fmul(ay, up[0]) - back[0],
cglib.fmul(ax, right[1]) + cglib.fmul(ay, up[1]) - back[1],
cglib.fmul(ax, right[2]) + cglib.fmul(ay, up[2]) - back[2]))
def magnitude_chain():
"""除法链的位宽上界：**复用附录 A.2 的扫描函数**，保证两边口径完全一致。
上界**取决于光线集合**：中间量随 `d × e2` 而变，只取一条探针光线会把上界算错。
故这里不自己造探针，直接调用附录 A 的 `scan_magnitudes()`——它按 `render_trace`
的公式在 320×240 上每 8 像素抽一条主光线（1200 条），对全部 896 个三角形扫描，
相机基取**行**。
"""
import cg_trace_appendix_a as appendix_a
nray, _ntri, _coord, max_mul, max_div, max_det = appendix_a.scan_magnitudes()
return nray, max_mul, max_div, max_det
def det_sweep():
"""命中点固定在三角形内部，方向逐渐共面：量 (t, u, v) 相对误差随 |det| 变化。"""
v0, v1, v2 = (0, 0, 0), (ONE, 0, 0), (0, ONE, 0)
center = (ONE // 4, ONE // 4, 0)
rows = []
for n, dn in [(1, 1), (1, 2), (1, 4), (1, 10), (1, 20), (1, 50), (1, 100), (1, 200),
(1, 500), (1, 1000)]:
dz = (n << 16) // dn
d = cglib.v_norm((int(0.3 * ONE), int(0.4 * ONE), dz))
orig = (center[0] - d[0], center[1] - d[1], center[2] - d[2])
e1, e2 = cglib.v_sub(v1, v0), cglib.v_sub(v2, v0)
det = cglib.v_dot(e1, cglib.v_cross(d, e2))
got = cgtrace.intersect_triangle(orig, d, v0, v1, v2)
exact = mt_exact(orig, d, v0, v1, v2)
if got is None or exact is None:
rows.append({"dzQ16": dz, "det": det, "hit": False})
continue
t_fx, u_fx, v_fx = got
t_ex, u_ex, v_ex = exact
errs = [abs(F(fx, ONE) - ex) / abs(ex)
for fx, ex in ((u_fx, u_ex), (v_fx, v_ex), (t_fx, t_ex)) if ex != 0]
rows.append({"dzQ16": dz, "det": det, "absDet": abs(det), "hit": True, "tQ16": t_fx,
"uQ16": u_fx, "vQ16": v_fx,
"relErrPermille": round(float(max(errs)) * 1000.0, 4)})
return rows
def real_det_distribution(spec, tris, nodes, order):
"""真实命中光线的 |det| 分布：最小 det 与几个阈值以下的命中数。"""
eye = tuple(spec["camera"]["eyeQ16"])
tmax = spec["camera"]["farQ16"]
min_det, hits = 1 << 62, 0
below = {1024: 0, 4096: 0, 16384: 0, 65536: 0}
for _px, _py, dirv in each_ray(spec):
hit = cgtrace.bvh_intersect(nodes, order, tris, eye, dirv, tmax)
if hit is None:
continue
hits += 1
t = tris[hit[1]]
e1 = cglib.v_sub((t[3], t[4], t[5]), (t[0], t[1], t[2]))
e2 = cglib.v_sub((t[6], t[7], t[8]), (t[0], t[1], t[2]))
det = abs(cglib.v_dot(e1, cglib.v_cross(dirv, e2)))
min_det = min(min_det, det)
for thr in below:
if det < thr:
below[thr] += 1
return {"hits": hits, "minAbsDet": min_det, "below": below}
def main():
spec = cglib.load_spec()
tris, _ = cgtrace.build_scene_tris(spec)
nodes, order = cgtrace.build_bvh(tris)
nray, max_mul, max_div, max_det = magnitude_chain()
print(f"== 2.2 除法链量级（{nray} 条抽样主光线 × {len(tris)} 个真实三角形）==")
print(" 口径与附录 A.2 相同：320×240 每 8 像素一条；相机基取行")
for name, value, note in (("fmul 原始乘积", max_mul, "上界来自 num × fdiv(ONE, det)"),
("fdiv 分子 |a| << 16", max_div, "重心坐标与 t 的定点除法"),
("行列式 det", max_det, "e1 · (d × e2)")):
v31 = "SHIFT-BROKEN" if value >= LIMIT_SHIFT else "OK"
v53 = "OVERFLOW" if value >= LIMIT_NUMBER else "OK"
print(f"{name:<20} {value:>20} 2^{value.bit_length():<3} 2^53={v53:<9} "
f"2^31={v31:<13} {note}")
print(f"2^31 = {LIMIT_SHIFT} 2^53 = {LIMIT_NUMBER}")
print(f"fmul 原始乘积 / 2^53 = {max_mul * 1000 // LIMIT_NUMBER}‰"
f"（余量只剩 {LIMIT_NUMBER * 100 // max_mul / 100:.2f} 倍）")
print(f"fdiv 分子 / 2^31 = {max_div // LIMIT_SHIFT} 倍")
print()
sweep = det_sweep()
print("== 2.2 det 扫描：方向逐渐共面，|det| 与相对误差 ==")
print(f"{'dz(Q16.16)':>10} {'|det|':>10} {'2^n':>6} {'u(Q16)':>10} {'v(Q16)':>10} "
f"{'相对误差‰':>12}")
for r in sweep:
if not r["hit"]:
print(f"{r['dzQ16']:>10} {r['det']:>10} {'-':>6} {'-':>10} {'-':>10} "
f"{'无交(det=0)':>12}")
continue
print(f"{r['dzQ16']:>10} {r['absDet']:>10} {('2^' + str(r['absDet'].bit_length())):>6} "
f"{r['uQ16']:>10} {r['vQ16']:>10} {r['relErrPermille']:>12}")
print()
dist = real_det_distribution(spec, tris, nodes, order)
print("== 2.2 真实命中光线的 |det| 分布（为什么不做 epsilon 剔除）==")
print(f"命中光线 : {dist['hits']}")
print(f"最小 |det| : {dist['minAbsDet']}（2^{dist['minAbsDet'].bit_length()}）")
for thr in sorted(dist["below"]):
print(f"|det| < {thr:<8} 的命中 : {dist['below'][thr]}")
print()
os.makedirs(OUT, exist_ok=True)
payload = {"magnitude": {"fmulRaw": max_mul, "fdivNum": max_div, "det": max_det},
"detSweep": sweep, "realDet": dist}
with open(os.path.join(OUT, "cg_trace_02_2_stability.json"), "w", encoding="utf-8") as f:
json.dump(payload, f, ensure_ascii=False, indent=2)
f.write("\n")
print("落盘 out/cg_trace_02_2_stability.json")
return 0
if __name__ == "__main__":
sys.exit(main())
跑一次（本机约 6 秒）：
# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_02_2_stability.py
== 2.2 除法链量级（1200 条抽样主光线 × 896 个真实三角形）==
口径与附录 A.2 相同：320×240 每 8 像素一条；相机基取行
fmul 原始乘积 5714018048147456 2^53 2^53=OK 2^31=SHIFT-BROKEN 上界来自 num × fdiv(ONE, det)
fdiv 分子 |a| << 16 549345492992 2^39 2^53=OK 2^31=SHIFT-BROKEN 重心坐标与 t 的定点除法
行列式 det 671983 2^20 2^53=OK 2^31=OK e1 · (d × e2)
2^31 = 2147483648 2^53 = 9007199254740992
fmul 原始乘积 / 2^53 = 634‰（余量只剩 1.57 倍）
fdiv 分子 / 2^31 = 255 倍
== 2.2 det 扫描：方向逐渐共面，|det| 与相对误差 ==
dz(Q16.16) |det| 2^n u(Q16) v(Q16) 相对误差‰
65536 58618 2^16 16384 16385 0.061
32768 46342 2^16 16384 16385 0.061
16384 29310 2^15 16385 16385 0.061
6553 12852 2^14 16389 16389 0.3052
3276 6520 2^13 16394 16394 0.6104
1310 2618 2^12 16396 16396 0.7324
655 1309 2^11 16421 16471 5.3101
327 654 2^10 16534 16534 9.1553
131 262 2^9 16759 16509 22.8882
65 130 2^8 17140 17140 46.1426
== 2.2 真实命中光线的 |det| 分布（为什么不做 epsilon 剔除）==
命中光线 : 43731
最小 |det| : 1102（2^11）
|det| < 1024 的命中 : 0
|det| < 4096 的命中 : 12
|det| < 16384 的命中 : 378
|det| < 65536 的命中 : 2775
落盘 out/cg_trace_02_2_stability.json
量级链：原始乘积已摸到 2^53
fmul 的原始乘积（叉积与点积的逐分量积，移位之前）最大到 5714018048147456，即 2^53——离 Number 能精确表示整数的上界 9007199254740992 只剩 1.57 倍余量（占比 634‰）。fdiv 的分子（|a| << 16）最大到 549345492992，即 2^39，是 2^31 的 255 倍；行列式 det 本身只到 671983（2^20）——问题全在乘法的中间量上，不在 det 本身。两条独立理由指向同一个结论：Number 的 >> 先 ToInt32 再移位，越过 2^31 会静默出错；原始乘积又已逼近 2^53，任何用浮点承载中间量的写法都站在悬崖边上。
上界取决于光线集合：中间量随 d × e2 而变，只取一条探针光线会把上界算错。所以本节不自己造探针，而是直接复用附录 A.2 的扫描函数——320×240 每 8 像素一条主光线（1200 条），对全部 896 个三角形扫描，相机基取行。
det → 0：相对误差从 0.061‰ 涨到 46.1426‰
扫描把命中点固定在三角形内部（重心 (0.25, 0.25)），让方向逐渐与三角形平面共面，det 从 2^16 掉到 2^8。用 fractions.Fraction 做精确有理数参照，量定点 (t, u, v) 的相对误差：2^16 时 0.061‰，2^12 时 0.7324‰，2^9 时 22.8882‰，2^8 时 46.1426‰——误差随 |det| 单调增长：定点除法的分子被左移 16 位后是固定的，分母越小，商的绝对误差越大。det = 0 时原语直接返回 None，不会给出被 0 除的噪声值。
为什么原语不做 epsilon 剔除
真实分布否掉了一个自然的「防御」——在 det 太小时直接判无交：43731 条命中光线里最小 |det| 是 1102（2^11），|det| < 1024 的命中是 0 条，< 4096 有 12 条，< 16384 有 378 条，< 65536 有 2775 条。这些小 det 的命中在几何上是真的——光线确实以掠射角打到了那个三角形；按 det < ε 剔除，这 12 条（乃至 378 条、2775 条）会凭空消失，画面出现漏光的洞。原语的正确做法是如实返回，把「要不要信这条命中」交给调用方：调用方知道自己能容忍多大误差，原语不知道。这也是 2.3 的伏笔——容差属于调用方，不属于求交本身。
跨语言核对
P 层在 TS 侧用 BigInt 重跑同一套 1200 条抽样主光线与同一次 det 扫描，逐位比对：
// snippet: 节选自 lantern/core/test/cg_trace_02_2_stability.test.ts
// 与附录 A.2 同口径：320×240 每 8 像素一条主光线（1200 条），相机基取**行**
for (let py = 0; py < height; py += 8) {
const ay = fmul(fdiv(fFromInt(2 * (halfH - py) - 1), fFromInt(height)), tanHalf);
for (let px = 0; px < width; px += 8) {
const ax = fmul(fdiv(fFromInt(2 * (px - halfW) + 1), fFromInt(width)), tanX);
const dirv = vNorm([
fmul(ax, right[0]) + fmul(ay, up[0]) - back[0],
fmul(ax, right[1]) + fmul(ay, up[1]) - back[1],
fmul(ax, right[2]) + fmul(ay, up[2]) - back[2],
]);
nray += 1;
for (const t of tris) {
const e1 = vSub([t[3], t[4], t[5]], [t[0], t[1], t[2]]);
const e2 = vSub([t[6], t[7], t[8]], [t[0], t[1], t[2]]);
const p = vCross(dirv, e2);
for (let i = 0; i < 3; i += 1) {
const m = e1[i] * p[i]; // 移位前的原始乘积（还有叉积的六个）
const am = m < 0n ? -m : m;
if (am > maxMul) maxMul = am;
}
const det = vDot(e1, p); // vDot 走 fmul，已 >>16
const ad = det < 0n ? -det : det;
if (ad > maxDet) maxDet = ad;
if (det === 0n) continue;
const num = vDot(vSub(eye, v0), p);
const div = (num < 0n ? -num : num) << 16n; // fdiv 的分子，是 2^31 的 255 倍
if (div > maxDiv) maxDiv = div;
const mm = num * fdiv(ONE, det); // 上界正出在这一项：2^53
const amm = mm < 0n ? -mm : mm;
if (amm > maxMul) maxMul = amm;
}
}
}
assert.equal(nray, 1200, "抽样主光线数应为 1200");
assert.equal(Number(maxMul), s.magnitude.fmulRaw, "fmul 原始乘积上界不符");
assert.ok(maxMul < (1n << 53n) && maxMul > (1n << 52n), "原始乘积应逼近 2^53");
三组断言都过：1200 条抽样光线的 fmulRaw / fdivNum / det 与 S 侧相同，且前两者都 > 2^31、原始乘积落在 (2^52, 2^53) 内（证明「必须 BigInt」不是空话）；det 扫描每一行的 (det, u, v, t) 与 S 侧逐位相同，且相对误差随 |det| 单调不减；真实命中的最小 |det| 与四个阈值计数逐条相同，并断言「真实命中里应存在 |det| < 4096 的光线」——这条断言直接钉住「不能剔除」。
阅读导航：上一节：2.1 球、平面与三角形 · 下一节：2.3 容差与自交 。
继续阅读
探索更多技术文章
浏览归档，发现更多关于系统设计、工具链和工程实践的内容。