13.1 求交的容差边界

|x_fx - x_exact| <= 1 + 2^16·δN/|det| + 2^16·|N_fx|·δdet/det² + |N_fx|·δinv/2^16

#!/usr/bin/env python3
# file: lantern/sim/cg_trace_13_1_tolerance_bound.py
# tier: S
# verify: /opt/homebrew/bin/python3 sim/cg_trace_13_1_tolerance_bound.py
"""13.1 求交的容差边界：Möller–Trumbore 的闭式误差界与「两宿主何时给出不同 t」。

第 2.2 节量了 det -> 0 时 (t, u, v) 的相对误差，第 2.3 节量了 EPSILON 被哪两条边界夹住。
本节把两件事接起来：给出一条**闭式的绝对误差界**，并回答**什么条件下两个宿主会给出不同
的 t**。MT 每一步都是一次 `fmul`（`(a*b) >> 16`，向下取整）或一次 `fdiv`（向零截断），
每次降精度丢掉不到 1 个 LSB；沿「叉积 -> 点积 -> 倒数 -> 最终乘」把损失加起来：

    |x_fx - x_exact| <= 1 + 2^16·δN/|det| + 2^16·|N_fx|·δdet/det² + |N_fx|·δinv/2^16

精确参照用**整数有理数**（不做任何截断的原始叉积/点积）算，无浮点、无 Fraction。在全部
43731 条命中主光线上核对「界 >= 实测」。输入全整数时运算序列完全确定，唯一的自由度是每
一步的**降精度规则**——这就是「两个宿主给出不同 t」的充要条件，用两个舍入变体量出来。
最后把 Δt 投影到**面法线**（`Δt·sinθ`，`sinθ = |det|/|e1×e2|`），接上第 2.3 节的 EPSILON。
"""

import json
import os
import sys

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
FMUL, FDIV = cglib.fmul, cglib.fdiv
BASELINE_SHA = "213b413069ca66c9f0cc091bce5328eac0f540d075afa3bab37788aea335ba71"


def raw_cross(a, b):                    # 无截断的叉积：整数 × 整数，保留全部位（Q32.32）
    return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0])


def raw_dot(a, b):
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]


def err_bound(det_fx, n_fx, d_n, d_det):
    """闭式误差界（Q16.16 的 LSB）：1（末次 fmul 的向下取整）+ 2^16·δN/|det|（分子截断）
    + 2^16·|N|·δdet/det²（det 截断被 1/det 放大）+ |N|·δinv/2^16（倒数本身的截断）。"""
    ad = abs(det_fx)
    return (1.0 + (ONE * d_n) / ad + (ONE * abs(n_fx) * d_det) / (ad * ad)
            + (abs(n_fx) * (1.0 + (1 << 32) * d_det / (ad * ad))) / ONE)


def analyze(orig, dirv, v0, v1, v2):
    """一次 MT：定点结果、精确有理参照、误差界，以及两个舍入变体的 t。"""
    e1 = (v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2])
    e2 = (v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2])
    tvec = cglib.v_sub(orig, v0)
    p = cglib.v_cross(dirv, e2)                 # 定点链，与 cgtrace.intersect_triangle 同式
    det = cglib.v_dot(e1, p)
    if det == 0:
        return None
    inv = FDIV(ONE, det)
    q = cglib.v_cross(tvec, e1)
    nu_fx, nv_fx, nt_fx = (cglib.v_dot(tvec, p), cglib.v_dot(dirv, q), cglib.v_dot(e2, q))
    P, Q = raw_cross(dirv, e2), raw_cross(tvec, e1)   # 精确参照：同一套整数输入，全无截断
    det_r, nu_r, nv_r, nt_r = (raw_dot(e1, P), raw_dot(tvec, P),
                               raw_dot(dirv, Q), raw_dot(e2, Q))
    s = lambda v: 3.0 + 2.0 * (abs(v[0]) + abs(v[1]) + abs(v[2])) / ONE
    d_det = s(e1)

    def abs_err(x_fx, n_r):
        return abs(x_fx * det_r - (n_r << 16)) / abs(det_r)   # 整数分子，最后一步才除

    area2 = cglib.v_len(cglib.v_cross(e1, e2))     # |e1×e2|（Q16.16，2 倍面积）
    t_err = abs_err(FMUL(nt_fx, inv), nt_r)
    return {
        "det": det, "tErr": t_err, "uErr": abs_err(FMUL(nu_fx, inv), nu_r),
        "vErr": abs_err(FMUL(nv_fx, inv), nv_r),
        "tBound": err_bound(det, nt_fx, s(e2), d_det),
        "uBound": err_bound(det, nu_fx, s(tvec), d_det),
        "vBound": err_bound(det, nv_fx, s(dirv), d_det),
        "perpErr": t_err * abs(det) / area2 if area2 else 0.0,
        "tFx": FMUL(nt_fx, inv),
        "tFloor": FMUL(nt_fx, (ONE * ONE) // det),          # 倒数改 floor（Python `//`）
        "tRound": (nt_fx * inv + (ONE >> 1)) >> 16,         # 末步四舍五入
    }


def primary_hits(spec, tris, nodes, order):
    """跑一遍主光线，对每条命中重做一次带误差界的 MT 分析。"""
    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"]
    th, tmax = cam["tanHalfFovYQ16"], cam["farQ16"]
    tx = FMUL(th, FDIV(cglib.f_from_int(w), cglib.f_from_int(h)))
    rows = []
    for py in range(h):
        ay = FMUL(FDIV(cglib.f_from_int(2 * (h // 2 - py) - 1), cglib.f_from_int(h)), th)
        for px in range(w):
            ax = FMUL(FDIV(cglib.f_from_int(2 * (px - w // 2) + 1), cglib.f_from_int(w)), tx)
            dirv = cglib.v_norm((
                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]))
            hit = cgtrace.bvh_intersect(nodes, order, tris, eye, dirv, tmax)
            if hit is None:
                continue
            rec = tris[hit[1]]
            r = analyze(eye, dirv, (rec[0], rec[1], rec[2]), (rec[3], rec[4], rec[5]),
                        (rec[6], rec[7], rec[8]))
            if r is not None:
                rows.append(r)
    return rows


def main():
    spec = cglib.load_spec()
    tris, _ = cgtrace.build_scene_tris(spec)
    nodes, order = cgtrace.build_bvh(tris)
    rows = primary_hits(spec, tris, nodes, order)

    n = len(rows)
    max_t = max(r["tErr"] for r in rows)
    hold = all(r["tErr"] <= r["tBound"] and r["uErr"] <= r["uBound"] and r["vErr"] <= r["vBound"]
               for r in rows)
    worst = max(max(r["tErr"] / r["tBound"], r["uErr"] / r["uBound"], r["vErr"] / r["vBound"])
                for r in rows)
    min_det = min(abs(r["det"]) for r in rows)
    n_floor = sum(1 for r in rows if r["tFloor"] != r["tFx"])     # 舍入变体：t 会变的命中数
    n_round = sum(1 for r in rows if r["tRound"] != r["tFx"])
    n_negdet = sum(1 for r in rows if r["det"] < 0)
    max_perp = max(r["perpErr"] for r in rows)                   # Δt vs 垂直误差 Δt·sinθ
    n_t_eps = sum(1 for r in rows if r["tErr"] > cgtrace.EPSILON)
    n_perp_eps = sum(1 for r in rows if r["perpErr"] > cgtrace.EPSILON)

    buckets = {}                                                 # |det| 分箱
    for r in rows:
        cur = buckets.setdefault(abs(r["det"]).bit_length(),
                                 {"count": 0, "maxErr": 0.0, "maxBound": 0.0, "minDet": 1 << 62})
        cur["count"] += 1
        cur["maxErr"] = max(cur["maxErr"], r["tErr"])
        cur["maxBound"] = max(cur["maxBound"], r["tBound"])
        cur["minDet"] = min(cur["minDet"], abs(r["det"]))
    bucket_rows = [(b, buckets[b]) for b in sorted(buckets)]

    print(f"== 13.1 求交的容差边界（{n} 条命中主光线）==")
    print(f"命中主光线        : {n}")
    print(f"最小 |det|        : {min_det}（2^{min_det.bit_length()}）")
    print(f"|Δt| 最大         : {max_t:.6f} LSB")
    print(f"闭式界是否恒成立   : {hold}（实测误差 / 界 的最大比值 {worst:.6f}）")
    print(f"\n== 舍入变体：改一处规则，t 会变的命中数 ==\n命中里 det < 0 的  : {n_negdet} / {n}")
    print(f"倒数改 floor 除法  : {n_floor} 条 t 改变\n末步改四舍五入     : {n_round} 条 t 改变")
    print("\n== 容差的真正尺度：沿光线的 Δt vs 垂直的 Δt·sinθ ==")
    print(f"Δt > EPSILON(64)   : {n_t_eps} / {n}\n垂直误差 最大      : {max_perp:.6f} LSB")
    print(f"垂直误差 > 64      : {n_perp_eps} / {n}")
    print(f"\n|det| 分箱表（{len(bucket_rows)} 个箱）落盘 out/cg_trace_13_1_det_buckets.csv")

    os.makedirs(OUT, exist_ok=True)
    cglib.write_csv(
        os.path.join(OUT, "cg_trace_13_1_det_buckets.csv"),
        ["det_bits", "hits", "min_abs_det", "max_abs_dt_lsb", "bound_lsb"],
        [[b, cur["count"], cur["minDet"], round(cur["maxErr"], 6), round(cur["maxBound"], 6)]
         for b, cur in bucket_rows])
    payload = {
        "hits": n, "minAbsDet": min_det, "maxAbsDtLsb": round(max_t, 6),
        "boundHolds": hold, "maxErrOverBound": round(worst, 6),
        "detNegative": n_negdet, "tChangedByFloorDiv": n_floor, "tChangedByRoundHalfUp": n_round,
        "epsilon": cgtrace.EPSILON, "tErrGtEpsilon": n_t_eps,
        "tErrGt256": sum(1 for r in rows if r["tErr"] > 256),
        "maxPerpErrLsb": round(max_perp, 6), "perpErrGtEpsilon": n_perp_eps,
        "baselineSha": BASELINE_SHA,
        "buckets": [{"detBits": b, "hits": c["count"], "minAbsDet": c["minDet"],
                     "maxAbsDtLsb": round(c["maxErr"], 6), "boundLsb": round(c["maxBound"], 6)}
                    for b, c in bucket_rows],
    }
    with open(os.path.join(OUT, "cg_trace_13_1_tolerance.json"), "w", encoding="utf-8") as fh:
        json.dump(payload, fh, ensure_ascii=False, indent=2)
        fh.write("\n")
    print("\n落盘 out/cg_trace_13_1_tolerance.json、out/cg_trace_13_1_det_buckets.csv")
    return 0


if __name__ == "__main__":
    sys.exit(main())

# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_13_1_tolerance_bound.py

== 13.1 求交的容差边界（43731 条命中主光线）==
命中主光线        : 43731
最小 |det|        : 1102（2^11）
|Δt| 最大         : 1442.726772 LSB
闭式界是否恒成立   : True（实测误差 / 界 的最大比值 0.746916）

== 舍入变体：改一处规则，t 会变的命中数 ==
命中里 det < 0 的  : 41254 / 43731
倒数改 floor 除法  : 41254 条 t 改变
末步改四舍五入     : 22043 条 t 改变

== 容差的真正尺度：沿光线的 Δt vs 垂直的 Δt·sinθ ==
Δt > EPSILON(64)   : 4582 / 43731
垂直误差 最大      : 241.471057 LSB
垂直误差 > 64      : 909 / 43731

|det| 分箱表（10 个箱）落盘 out/cg_trace_13_1_det_buckets.csv

落盘 out/cg_trace_13_1_tolerance.json、out/cg_trace_13_1_det_buckets.csv

// snippet: 节选自 lantern/core/test/cg_trace_13_1_tolerance_bound.test.ts（BigInt 重实现）
// verify: node --test 'core/test/cg_trace_13_1_tolerance_bound.test.ts'
  assert.equal(rows.length, golden.hits);
  assert.equal(rows.filter((r) => r.det < 0n).length, golden.detNegative);
  assert.equal(rows.filter((r) => r.tFloor !== r.tFx).length, golden.tChangedByFloorDiv);
  assert.equal(rows.filter((r) => r.tRound !== r.tFx).length, golden.tChangedByRoundHalfUp);