11.1 SDF 的构造与布尔运算

球：    f(p) = |p - c| - r
盒：    f(p) = |max(q, 0)| + min(max(q.x, q.y, q.z), 0)，  q = |p - c| - b
平面：  f(p) = p.y - h
高度场：f(p) = p.y - h(x, z)

# snippet: sim/cg_trace_11_sdf.py（节选：定点长度与四个基本场）
# verify: /opt/homebrew/bin/python3 sim/cg_trace_11_1_sdf_construction.py
def len_q16(v):
    """Q16.16 向量的**精确**长度：`isqrt(v0^2 + v1^2 + v2^2)`（结果仍是 Q16.16）。"""
    return math.isqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])


def sdf_sphere(p, c, r):
    return len_q16(v_sub(p, c)) - r


def sdf_box(p, c, b):
    d = v_sub(p, c)
    q = (abs(d[0]) - b[0], abs(d[1]) - b[1], abs(d[2]) - b[2])
    return len_q16((max(q[0], 0), max(q[1], 0), max(q[2], 0))) + min(max(q[0], q[1], q[2]), 0)


def sdf_plane(p, height):
    """过 `y = height` 的平面（半空间）：`p.y - height`。精确距离场。"""
    return p[1] - height


def sdf_terrain(p, seed, hs):
    """高度场：`p.y - h(x,z)`。**不是**距离场——陡坡处 `|grad f| > 1`。"""
    return p[1] - cglib.terrain_height(p[0], p[2], seed, hs)

并集：   min(a, b)
交集：   max(a, b)
差集：   max(a, -b)          ← a 减去 b，不是 b 减去 a
平滑并： smin(a, b, k) = mix(b, a, h) - k·h·(1-h)，h = clamp(0.5 + 0.5·(b-a)/k, 0, 1)

// file: lantern/core/test/cg_trace_11_1_sdf_construction.test.ts
// tier: P
// verify: node --test 'core/test/cg_trace_11_1_sdf_construction.test.ts'
// 卷三第 11.1 节《SDF 的构造与布尔运算》的 P 层（上）：Lipschitz 界与下界性。
//
// 两个断言对象：
//   1. 八个场的 Lipschitz 估计（同一条 SplitMix64 采样序列，必须与 S 侧逐位相同）；
//   2. 四种可解析场的下界性 `max(f - d_true)`（球/盒/平面/并集都必须是 0）。
// 平滑并 `smin` 与高度场高估的近邻断言在 cg_trace_11_1_sdf_bounds.test.ts。
// 金标准由 S 侧产出：
//   /opt/homebrew/bin/python3 sim/cg_trace_11_1_sdf_construction.py

import assert from "node:assert/strict";
import { readFileSync } from "node:fs";
import { dirname, join } from "node:path";
import { test } from "node:test";
import { fileURLToPath } from "node:url";

import {
  type SdfConfig,
  lenQ16,
  sdfBox,
  sdfPlane,
  sdfSphere,
  sdfTerrain,
} from "../src/lantern/cg_trace_11_sdf.ts";
import { rngNext } from "../src/lantern/cgtrace.ts";
import { ONE, fdiv, fmul } from "../src/lantern/fixed.ts";
import { type Vec3, vAdd, vNorm, vSub } from "../src/lantern/vec.ts";

const root = join(dirname(fileURLToPath(import.meta.url)), "..", "..");
const cfg = JSON.parse(readFileSync(join(root, "config", "cg_trace_11_sdf.json"), "utf8")) as SdfConfig;
const golden = JSON.parse(readFileSync(join(root, "out", "cg_trace_11_1_sdf.json"), "utf8"));

const pr = cfg.primitives;
const SC: Vec3 = [BigInt(pr.sphere.centerQ16[0]), BigInt(pr.sphere.centerQ16[1]), BigInt(pr.sphere.centerQ16[2])];
const SR = BigInt(pr.sphere.radiusQ16);
const BC: Vec3 = [BigInt(pr.box.centerQ16[0]), BigInt(pr.box.centerQ16[1]), BigInt(pr.box.centerQ16[2])];
const BB: Vec3 = [BigInt(pr.box.halfQ16[0]), BigInt(pr.box.halfQ16[1]), BigInt(pr.box.halfQ16[2])];
const PH = BigInt(pr.groundPlane.heightQ16);
const SEED = BigInt(pr.terrain.noiseSeed);
const HS = BigInt(pr.terrain.heightScaleQ16);

type Field = (p: Vec3) => bigint;

const sphere: Field = (p) => sdfSphere(p, SC, SR);
const box: Field = (p) => sdfBox(p, BC, BB);
const plane: Field = (p) => sdfPlane(p, PH);

function opUnion(a: bigint, b: bigint): bigint {
  return a < b ? a : b;
}
function opIntersect(a: bigint, b: bigint): bigint {
  return a > b ? a : b;
}
function opDifference(a: bigint, b: bigint): bigint {
  return opIntersect(a, -b);
}

const fields: Record<string, Field> = {
  sphere,
  box,
  plane,
  terrain: (p) => sdfTerrain(p, SEED, HS),
  union: (p) => opUnion(sphere(p), box(p)),
  intersect: (p) => opIntersect(sphere(p), box(p)),
  difference: (p) => opDifference(sphere(p), box(p)),
};

/** 点对采样估计 Lipschitz 常数：`max |f(p)-f(q)| / |p-q|`（Q16.16）。 */
function lipschitz(f: Field): bigint {
  const ext = BigInt(cfg.gradient.extentQ16);
  const delta = BigInt(cfg.gradient.deltaQ16);
  const samples = cfg.gradient.lipschitzSamples;
  let state = BigInt(cfg.gradient.lipschitzSeed);
  const span = 2n * ext + 1n;
  const twoOne = 2n * ONE + 1n;
  let mx = -1n;
  for (let i = 0; i < samples; i += 1) {
    const pt: Vec3 = [0n, 0n, 0n];
    for (let k = 0; k < 3; k += 1) {
      const [s, z] = rngNext(state);
      state = s;
      pt[k] = (z % span) - ext;
    }
    const comp: bigint[] = [];
    for (let k = 0; k < 3; k += 1) {
      const [s, z] = rngNext(state);
      state = s;
      comp.push((z % twoOne) - ONE);
    }
    const u = vNorm([comp[0], comp[1], comp[2]]);
    if (u[0] === 0n && u[1] === 0n && u[2] === 0n) continue;
    const q: Vec3 = vAdd(pt, [fmul(delta, u[0]), fmul(delta, u[1]), fmul(delta, u[2])]);
    const dist = lenQ16(vSub(q, pt));
    if (dist === 0n) continue;
    const df = f(pt) - f(q);
    const r = fdiv(df < 0n ? -df : df, dist);
    if (r > mx) mx = r;
  }
  return mx;
}

test("11.1 Lipschitz：八个场的点对估计与 S 层逐位相同", () => {
  for (const name of Object.keys(fields)) {
    assert.equal(lipschitz(fields[name]), BigInt(golden.lipschitz[name].lipschitzQ16), name);
  }
  // 可解析场（球/盒/平面/并集）的 Lipschitz 常数必须落在 1.0 附近
  for (const name of ["sphere", "box", "plane", "union"]) {
    const L = lipschitz(fields[name]);
    assert.ok(L >= ONE && L <= ONE + 64n, `${name} L=${L}`);
  }
  // 高度场必须明显越界，否则 11.2 的「穿透」一列无从解释
  assert.ok(lipschitz(fields.terrain) > ONE + 4096n, "高度场必须明显越界");
});

test("11.1 下界性：球/盒/平面/并集的 max(f - d_true) 必须为 0", () => {
  const ext = BigInt(cfg.gradient.extentQ16);
  const n = 12;
  const step = fdiv(2n * ext, BigInt(n) << 16n);
  const dTrue = (name: string, p: Vec3): bigint => {
    if (name === "sphere") return sphere(p);
    if (name === "box") return box(p);
    if (name === "plane") return plane(p);
    const a = sphere(p);
    const b = box(p);
    return a < b ? a : b;
  };
  for (const name of ["sphere", "box", "plane", "union"]) {
    let over = 0;
    let worst = -(1n << 60n);
    for (let iz = 0; iz <= n; iz += 1) {
      for (let iy = 0; iy <= n; iy += 1) {
        for (let ix = 0; ix <= n; ix += 1) {
          const p: Vec3 = [
            -ext + fmul(step, BigInt(ix) << 16n),
            -ext + fmul(step, BigInt(iy) << 16n),
            -ext + fmul(step, BigInt(iz) << 16n),
          ];
          const excess = fields[name](p) - dTrue(name, p);
          if (excess > 0n) over += 1;
          if (excess > worst) worst = excess;
        }
      }
    }
    assert.equal(worst, BigInt(golden.lowerBound[name].maxExcessQ16), `${name} 最大超出`);
    assert.equal(over, golden.lowerBound[name].overCount, `${name} 高估点数`);
    assert.equal(over, 0, `${name} 不允许高估`);
  }
});

# snippet: 本机实跑（工作区 lantern/）
cd lantern
/opt/homebrew/bin/python3 sim/cg_trace_11_1_sdf_construction.py
node --test 'core/test/cg_trace_11_1_*.test.ts'

== 11.1 SDF 的构造与布尔运算 ==
Lipschitz 点对采样 40000 对，步长 16384 Q16.16
场            Lipschitz L（Q16.16 / 千分比）   max(f - d_true)
  sphere         65536 /   1000.0‰      0
  box            65540 / 1000.061‰      0
  plane          65536 /   1000.0‰      0
  terrain        71048 / 1084.1064‰      —
  union          65540 / 1000.061‰      0
  intersect      65536 /   1000.0‰      —
  difference     65536 /   1000.0‰      —
  smin           65540 / 1000.061‰      —
高度场高估：426/648 个采样点 f > d_true，最大高估 3909 Q16.16（比值上界 1056.6406‰）
smin <= min 违反 0 次 / 4913 点；最大间隙 min-smin = 8187 Q16.16
csv/json          : out/cg_trace_11_1_{lipschitz,sdf_slice}.csv / _sdf.json