Symbolic Rendering

Appendix B — The Machine-Checked Archive

An explainer comparing numerical rendering with symbolic rendering,
in the spirit of Utah graphics meets Macsyma (today: Python + sympy).
Dr. L. Van Warren — Original Idea & Mathematical Development
Kimi K3 — Assistance with Mathematical Development
Single HTML5 file · the archive itself is the notebook (⬇ below) · Requires an internet connection for mathematical typesetting (MathJax).
❦
How to use this file. This appendix is notebook-native by design: the artifact is the ⬇ full archive notebook below — twelve cells, run in order, top to bottom, in a fresh Colab session. §1 states the rules of evidence (three tolerance classes, named and fenced). §2 is the cell map. The laboratory hosts the twelve cells with their three run paths each. §3 records what the archive cannot certify. Then the book closes.
☰  contents
  1. The rules of evidence
  2. The cell map
  3. The laboratory: the archive, cell by cell
  4. What the archive cannot certify
  5. The book closes

§1 · The rules of evidence

Every check in the archive declares which court it is tried in. There are three, and conflating them is how verification suites learn to lie:

ClassToleranceUsed forExamples
Exact==, or contradiction of a nonzero constantidentities in exact arithmetic; determinism; set membershipstrata reproducibility; gamma non-commutativity; bitmap sync
Machine≤ 10⁻¹² … 10⁻⁹equalities that hold in ℝ, evaluated in float64closed forms (ν, r·V̂); ∂N′/∂β; mirror flip; bokeh law
Statisticalscaled to the discretization, printed with the numberquadrature comparisons; scan verifications; mass-config censuseschecker integral vs 64² quadrature; switching times vs 400-pt scans; the Shadow Lemma's 200,000 configurations

Statistical-class tolerances are never silently small: each check prints its measured number and its skip policy (the Shadow Lemma's marginal-band exclusion is printed, not implied). All seeds are fixed; a rerun reproduces the archive bit-for-bit on the same runtime. A red line means the environment or the edit broke something — the theorems themselves live in the chapters' §3 sections, and no amount of green here substitutes for them.

§2 · The cell map

CellChapter batteryChecks
0Framework & the shared scene — check registry, compact renderer, both shadow formulations—
1One sphere, two computations — the IEEE trap; ν and r·V̂ closed forms; the silhouette radius4
2Half-vectors & two lights — disable-by-substitution; color affinity; Blinn's constant-h special case3
3Shadows — the Shadow Lemma, 200,000 configurations; chamber partition; the boundary is a curve3
4Reflection — mirror identities; the factoring self-hit; mirror flip ≡ bounce path; the ε-footprint exists4
5Texture — chart roundtrip; longitude theorem; Nz(r²); sRGB roundtrip; bitmap sync; gamma order; the Nyquist knee7
6Antialiasing — the checker box integral; the thin band; filter/encode non-commutativity; supersampling converges4
7Area lights — the on-axis anchor; the ρ→0 regression; the conic classified; variance has a seed, bias doesn't6
8Gloss — exact Lambert weights; the closed-form normalizer; partition of unity; mirror-grade concentration4
9Bump — the U-curve's two slopes; ∂N′/∂β vs FD; the β→0 regression; the bump is not a no-op4
10Motion & depth — switching times vs scans; the bokeh law, exact; the v→0 regression; quartic shadows and the quartic outline, certified by sympy5
11The tally — the archive counts itself and prints its verdict44 total

§B·3 · The laboratory: the archive, cell by cell

Run top to bottom in a fresh session (cells are cumulative, per series convention). Three run paths per cell — 📋 Copy code, ⬇ notebook (File → Upload notebook in Colab), or paste into colab.new.

🐍 Cell 0 — the framework and the shared scene

· · paste into colab.new → Run
import numpy as np
import math

# --- the check registry: every battery writes here -------------------------------
RESULTS = []
def check(ch, name, passed, detail=""):
    RESULTS.append((ch, name, bool(passed), detail))
    print(f"  [{'PASS' if passed else 'FAIL'}] ch{ch}: {name}"
          + (f"  ({detail})" if detail else ""))

def banner(ch, title):
    print(f"\n=== Chapter {ch} — {title} ===")

# --- the series' default scene, compact ------------------------------------------
SC = dict(R=1.0, e=4.0, Cb=np.array([0.5, 0.5, 1.3]), R2=0.32,
          lights=[(np.array([5., 5., 10.]), np.array([1., 1., 1.])),
                  (np.array([-6., 2., 8.]), np.array([0.9, 0.65, 0.25]))],
          ambient=np.array([0.15, 0.15, 0.16]),
          matA=np.array([0.63, 0.19, 0.16]), matB=np.array([0.33, 0.44, 0.62]),
          ka=0.1, kd=0.7, ks=0.5, shin=32.0, alpha=0.25, ku=8.0, kv=4.0)

def srgb(c):
    c = np.maximum(np.asarray(c, float), 0.0)
    return np.where(c <= 0.0031308, 12.92*c, 1.055*c**(1/2.4) - 0.055)

def srgb_inv(c):
    c = np.asarray(c, float)
    return np.where(c <= 0.04045, c/12.92, ((c + 0.055)/1.055)**2.4)

def chartA(P):
    N = P/SC['R']
    return ((math.atan2(N[1], N[0]) + math.pi)/(2*math.pi),
            math.acos(max(-1.0, min(1.0, N[2])))/math.pi)

def texA(u, v):
    s = math.sin(math.pi*u*SC['ku'])*math.sin(math.pi*v*SC['kv'])
    return SC['matA']*(1.0 if s >= 0 else 0.2)

def floor_tex(x, y):
    base = 0.85 if math.sin(math.pi*x)*math.sin(math.pi*y) >= 0 else 0.15
    return np.array([base*(1 - SC['alpha'])]*2 + [base])

# --- BOTH shadow formulations (the archive cross-examines them against each other) --
def shadowed_roots(P, Lp, C, Rr):       # numeric formulation: explicit roots
    seg = Lp - P
    a2 = float(seg @ seg); b2 = float((P - C) @ seg)
    c2 = float((P - C) @ (P - C)) - Rr*Rr
    if a2 < 1e-18: return False
    dsc = b2*b2 - a2*c2
    if dsc >= 0:
        sq = math.sqrt(dsc)
        s1, s2 = (-b2 - sq)/a2, (-b2 + sq)/a2
        if (0 < s1 < 1) or (0 < s2 < 1): return True
    return False

def shadowed_signs(P, Lp, C, Rr):       # symbolic formulation: the Shadow Lemma
    seg = Lp - P
    a2 = float(seg @ seg); b2 = float((P - C) @ seg)
    c2 = float((P - C) @ (P - C)) - Rr*Rr
    if a2 < 1e-18: return False
    dsc = b2*b2 - a2*c2
    return dsc >= 0 and b2 < 0 and a2 + b2 > 0

def shade(P, N, v_hat, alb, blockers, lights=None):
    col = SC['ka']*SC['ambient']
    for pos, color in (SC['lights'] if lights is None else lights):
        if color is None: continue
        if any(shadowed_roots(P, pos, C, Rr) for C, Rr in blockers): continue
        w = pos - P; w_hat = w/np.linalg.norm(w)
        mu = float(N @ w_hat)
        h = w_hat + v_hat; h = h/np.linalg.norm(h)
        spec = max(0.0, float(N @ h))**SC['shin']
        col += color*(SC['kd']*alb*max(0.0, mu) + SC['ks']*spec)
    return col

def render(W, lights=None, beta=0.0, grad=None):
    """Compact point-light two-sphere textured renderer (the Ch.5 scene)."""
    img = np.zeros((W, W, 3)); code = np.zeros((W, W), np.uint8)
    e, R, Cb, R2 = SC['e'], SC['R'], SC['Cb'], SC['R2']
    half = 1.25*e*R/math.sqrt(e*e - R*R)
    E = np.array([0., 0., e])
    for j in range(W):
        y = half*(1 - 2*(j + .5)/W) - 0.15*half
        for i in range(W):
            x = half*(2*(i + .5)/W - 1)
            Aq = x*x + y*y + e*e
            dA = e*e*R*R - (e*e - R*R)*(x*x + y*y)
            tA = (e*e - math.sqrt(dA))/Aq if dA >= 0 else math.inf
            bB = -Cb[0]*x - Cb[1]*y - (e - Cb[2])*e
            cB = float(Cb @ Cb) - R2*R2
            dB = bB*bB - Aq*cB
            tB = (-bB - math.sqrt(dB))/Aq if dB >= 0 else math.inf
            tm = min(tA, tB, 1.0)
            if tm == math.inf: continue
            P = np.array([tm*x, tm*y, e*(1 - tm)])
            if tm == tA:
                N = P/R; u, v = chartA(P); alb = texA(u, v)
                if beta and grad is not None:
                    hu, hv = grad(u, v)
                    sp_ = math.sqrt(max(0.0, 1 - N[2]**2))
                    if sp_ > 1e-4:
                        eT = np.array([-N[1]/sp_, N[0]/sp_, 0.])
                        eP = np.array([N[2]*N[0]/sp_, N[2]*N[1]/sp_, -sp_])
                    else:
                        eT = np.array([1., 0., 0.]); eP = np.array([0., 1., 0.])
                    M = N + beta*(hu*eT + hv*eP)
                    N = M/np.linalg.norm(M)
                blockers = [(Cb, R2)]; kind = 1
            elif tm == tB:
                N = (P - Cb)/R2; alb = SC['matB']
                blockers = [(np.zeros(3), R)]; kind = 2
            else:
                N = np.array([0., 0., 1.]); alb = floor_tex(P[0], P[1])
                blockers = [(np.zeros(3), R), (Cb, R2)]; kind = 3
            v_hat = (E - P)/np.linalg.norm(E - P)
            img[j, i] = shade(P, N, v_hat, alb, blockers, lights)
            code[j, i] = kind
    return np.clip(img, 0, 1), code

print("framework loaded. scene: two spheres, two point lights, checker floor.")
print(f"scene check: silhouette rho = "
      f"{SC['e']*SC['R']/math.sqrt(SC['e']**2 - SC['R']**2):.6f} (expect 1.032796)")
Expect: the framework banner and ρ = 1.032796 — the Chapter-1 constant, printed as a handshake.

🐍 Cell 1 — Chapter 1 battery: the pitfall, the closed forms, the silhouette

· · paste into colab.new → Run
banner(1, "one sphere, two computations")
rng = np.random.default_rng(1)

with np.errstate(invalid='ignore'):
    trap = 0.0*np.sqrt(-1.0)
check(1, "IEEE: 0*sqrt(-1) is NaN (the raw-Heaviside trap is real)",
      bool(np.isnan(trap)), f"value = {trap}")

e, R = SC['e'], SC['R']
x = rng.uniform(-0.9, 0.9, 20000); y = rng.uniform(-0.9, 0.9, 20000)
Aq = x*x + y*y + e*e
disc = e*e*R*R - (e*e - R*R)*(x*x + y*y)
m = disc >= 0
x, y, disc, Aq = x[m], y[m], disc[m], Aq[m]
t = (e*e - np.sqrt(disc))/Aq
P = np.stack([t*x, t*y, e*(1 - t)], -1)
N = P/R
L = SC['lights'][0][0]
w = L[None, :] - P; w /= np.linalg.norm(w, axis=1, keepdims=True)
V = np.stack([-x, -y, np.full_like(x, e)], -1)
V /= np.linalg.norm(V, axis=1, keepdims=True)
mu = np.einsum('ij,ij->i', N, w)
nu_direct = np.einsum('ij,ij->i', N, V)
nu_closed = np.sqrt(disc)/(R*np.sqrt(Aq))
d_nu = np.abs(nu_direct - nu_closed).max()
check(1, "closed form: nu = sqrt(Delta)/(R*sqrt(A))", d_nu < 1e-12,
      f"worst {d_nu:.2e} over {m.sum()} hits")
rdv_direct = np.einsum('ij,ij->i', 2*mu[:, None]*N - w, V)
rdv_closed = 2*mu*nu_closed - np.einsum('ij,ij->i', w, V)
d_rv = np.abs(rdv_direct - rdv_closed).max()
check(1, "closed form: r.V = 2*mu*nu - w.V (no new sqrt)", d_rv < 1e-12,
      f"worst {d_rv:.2e}")

rho = e*R/math.sqrt(e*e - R*R)
d_on = e*e*R*R - (e*e - R*R)*rho*rho
check(1, "silhouette theorem: Delta = 0 exactly at r = rho",
      abs(d_on) < 1e-10, f"Delta(rho^2) = {d_on:.2e}, rho = {rho:.6f}")
Expect: four PASS lines, machine-class numbers.

🐍 Cell 2 — Chapter 2 battery: substitution, affinity, the constant-h special case

· · paste into colab.new → Run
banner(2, "half-vectors, three channels, two lights")

img_sub, _ = render(64, lights=[(SC['lights'][0][0], SC['lights'][0][1]),
                                (SC['lights'][1][0], np.zeros(3))])
img_omit, _ = render(64, lights=[SC['lights'][0]])
d_sub = np.abs(img_sub - img_omit).max()
check(2, "disable-by-substitution (c2->0) == omitting the light",
      d_sub < 1e-12, f"max {d_sub:.2e}")

def diffuse_field(c1r):
    e_, R_ = SC['e'], SC['R']
    W = 64; half = 1.25*e_*R_/math.sqrt(e_*e_ - R_*R_)
    xs = np.linspace(-half + half/W, half - half/W, W)
    ys = np.linspace(half - half/W, -half + half/W, W) - 0.15*half
    X, Y = np.meshgrid(xs, ys)
    Aq = X**2 + Y**2 + e_**2
    disc = e_**2*R_**2 - (e_**2 - R_**2)*(X**2 + Y**2)
    t = (e_**2 - np.sqrt(np.maximum(disc, 0)))/Aq
    P_ = np.stack([t*X, t*Y, e_*(1 - t)], -1)
    N_ = P_/R_
    L_ = SC['lights'][0][0]
    w_ = L_[None, None, :] - P_
    w_ /= np.linalg.norm(w_, axis=-1, keepdims=True)
    mu_ = np.einsum('ijk,ijk->ij', N_, w_)
    return c1r*np.maximum(mu_, 0)*(disc >= 0)

I1, I2, I3 = diffuse_field(0.6), diffuse_field(0.8), diffuse_field(1.0)
lin = np.abs(I1 - 2*I2 + I3).max()
check(2, "the image is exactly affine in light color (2nd difference = 0)",
      lin < 1e-12, f"max {lin:.2e}")

L_far = np.array([1.0, 1.0, 1.0])*1e6
def hvec(P_):
    e_ = SC['e']
    w_ = L_far - P_; w_ /= np.linalg.norm(w_)
    V_ = np.array([0., 0., e_]) - P_; V_ /= np.linalg.norm(V_)
    h_ = w_ + V_; return h_/np.linalg.norm(h_)
h1 = hvec(np.array([0.2, 0.1, math.sqrt(1 - 0.05)]))
h2 = hvec(np.array([-0.3, 0.2, math.sqrt(1 - 0.13)]))
dh = np.abs(h1 - h2).max()
check(2, "Blinn's constant-h special case (directional light limit)",
      dh < 1e-5, f"|dh| = {dh:.2e} across the sphere at |L| = 1e6")
Expect: three PASS lines — the substitution exact, the second difference zero, h constant to 1e-6.

🐍 Cell 3 — Chapter 3 battery: the Shadow Lemma, 200,000 strong

· · paste into colab.new → Run
banner(3, "shadows: the quantifier eliminated")
rng = np.random.default_rng(7)

M = 200_000
P = rng.normal(size=(M, 3))*2.0
L = rng.normal(size=(M, 3))*4.0 + [0, 0, 6.0]
Cb_ = rng.normal(size=(M, 3))*1.5
Rb = rng.uniform(0.2, 1.5, size=M)
seg = L - P
A2 = np.einsum('ij,ij->i', seg, seg)
B2 = np.einsum('ij,ij->i', P - Cb_, seg)
C2 = np.einsum('ij,ij->i', P - Cb_, P - Cb_) - Rb**2
DSC = B2**2 - A2*C2
keep = (C2 >= 0) & (A2 + 2*B2 + C2 >= 0) & (np.abs(DSC) > 1e-9)
sq = np.sqrt(np.maximum(DSC, 0))
s1 = (-B2 - sq)/A2; s2 = (-B2 + sq)/A2
marginal = (np.abs(s1) < 1e-7) | (np.abs(s1 - 1) < 1e-7) | \
           (np.abs(s2) < 1e-7) | (np.abs(s2 - 1) < 1e-7)
keep &= ~marginal
roots = (DSC >= 0) & (((s1 > 0) & (s1 < 1)) | ((s2 > 0) & (s2 < 1)))
signs = (DSC >= 0) & (B2 < 0) & (A2 + B2 > 0)
agree = (roots == signs)[keep]
check(3, "Shadow Lemma: sign logic == explicit roots",
      bool(agree.all()),
      f"{int(agree.sum())}/{int(keep.sum())} configurations "
      f"(skipped {M - int(keep.sum())} containment/tangent/marginal)")

img3, code3 = render(128)
inhabited = sorted(set(code3.ravel().tolist()))
check(3, "chamber partition: every pixel coded; all four regions inhabited",
      np.isin(code3, [0, 1, 2, 3]).all() and len(inhabited) == 4,
      f"regions present: {inhabited}")

edge = np.zeros(code3.shape, bool)
edge[1:, :] |= code3[1:, :] != code3[:-1, :]
edge[:, 1:] |= code3[:, 1:] != code3[:, :-1]
frac = edge.mean()
check(3, "the branch boundary is a curve, not a region",
      frac < 0.06, f"flip-eligible pixels: {100*frac:.2f}% of frame")
Expect: 200,000-configuration agreement with the skip policy printed; the partition and thinness checks pass.

🐍 Cell 4 — Chapter 4 battery: identities, the factoring root, the mirror flip, the ε-footprint

· · paste into colab.new → Run
banner(4, "reflection: recursion becomes composition")
rng = np.random.default_rng(42)

D = rng.normal(size=(100_000, 3))
N = rng.normal(size=(100_000, 3)); N /= np.linalg.norm(N, axis=1, keepdims=True)
dn = np.einsum('ij,ij->i', D, N)
Dp = D - 2*dn[:, None]*N
r1 = np.abs(np.linalg.norm(Dp, axis=1) - np.linalg.norm(D, axis=1)).max()
r2 = np.abs(np.einsum('ij,ij->i', Dp, N) + dn).max()
check(4, "reflection identities |D'| = |D|, D'.N = -D.N (100k trials)",
      r1 < 1e-12 and r2 < 1e-12, f"worst {max(r1, r2):.2e}")

e, R = SC['e'], SC['R']
x = rng.uniform(-.8, .8, 5000); y = rng.uniform(-.8, .8, 5000)
Aq = x*x + y*y + e*e
disc = e*e*R*R - (e*e - R*R)*(x*x + y*y)
m = disc >= 0
x, y, disc, Aq = x[m], y[m], disc[m], Aq[m]
t = (e*e - np.sqrt(disc))/Aq
P = np.stack([t*x, t*y, e*(1 - t)], -1)
c_self = np.abs(np.einsum('ij,ij->i', P, P) - R*R).max()
check(4, "the self-hit quadratic factors: c' == 0 at machine precision",
      c_self < 1e-12, f"max |c'| = {c_self:.2e}")

# the mirror flip, on floor pixels: t_mirrored_scene == 1 + t_reflected
W = 64
half = 1.25*e*R/math.sqrt(e*e - R*R)
xs = np.linspace(-half + half/W, half - half/W, W)
ys = np.linspace(half - half/W, -half + half/W, W) - 0.15*half
X, Y = np.meshgrid(xs, ys)
bx, by, bz = SC['Cb']; R2 = SC['R2']
aa = X**2 + Y**2 + e**2
bb = (X - bx)*X + (Y - by)*Y - bz*e
cc = (X - bx)**2 + (Y - by)**2 + bz**2 - R2**2
dd = bb**2 - aa*cc
t_rfl = np.where(dd >= 0, (-bb - np.sqrt(np.maximum(dd, 0)))/aa, np.inf)
t_rfl = np.where(t_rfl > 0, t_rfl, np.inf)
bm = -bx*X - by*Y + (bz - e)*e
cm = bx**2 + by**2 + (bz - e)**2 - R2**2
dm = bm**2 - aa*cm
t_mir = np.where(dm >= 0, (-bm - np.sqrt(np.maximum(dm, 0)))/aa, np.inf)
t_mir = np.where(t_mir > 1, t_mir, np.inf)      # beyond the mirror plane (s = 1)
both = np.isfinite(t_rfl) & np.isfinite(t_mir)
gap = np.abs((t_mir[both] - 1) - t_rfl[both]).max() if both.any() else np.nan
check(4, "mirror flip == bounce path (t_mir = 1 + t_refl)",
      both.sum() > 0 and gap < 1e-12,
      f"{int(both.sum())} floor px, max |Delta t| = {gap:.2e}")

# the epsilon-hack footprint: measurable, bounded
nrm = np.sqrt(aa)
Qx = X + 1e-4*X/nrm; Qy = Y + 1e-4*Y/nrm; Qz = 1e-4*e/nrm
bbh = (Qx - bx)*X + (Qy - by)*Y + (Qz - bz)*e
cch = (Qx - bx)**2 + (Qy - by)**2 + (Qz - bz)**2 - R2**2
ddh = bbh**2 - aa*cch
t_h = np.where(ddh >= 0, (-bbh - np.sqrt(np.maximum(ddh, 0)))/aa, np.inf)
t_h = np.where(t_h > 0, t_h, np.inf)
both2 = np.isfinite(t_rfl) & np.isfinite(t_h)
shift = np.abs(t_h[both2] - t_rfl[both2])
check(4, "the epsilon-hack footprint exists and is bounded",
      both2.sum() > 0 and 1e-9 < shift.mean() < 1e-2,
      f"mean |Delta t| = {shift.mean():.2e} on {int(both2.sum())} px")
Expect: identities and factorization at 1e-12; the mirror flip exact; the hack's footprint a small but nonzero number.

🐍 Cell 5 — Chapter 5 battery: the chart, the bitmap, the gamma theorem, the knee

· · paste into colab.new → Run
banner(5, "texture: albedo becomes a function")
rng = np.random.default_rng(4242)

N = rng.normal(size=(100_000, 3)); N /= np.linalg.norm(N, axis=1, keepdims=True)
u = (np.arctan2(N[:, 1], N[:, 0]) + np.pi)/(2*np.pi)
v = np.arccos(np.clip(N[:, 2], -1, 1))/np.pi
th = 2*np.pi*u - np.pi; sv = np.sin(np.pi*v)
N2 = np.stack([np.cos(th)*sv, np.sin(th)*sv, np.cos(np.pi*v)], -1)
r = np.abs(N - N2).max()
check(5, "chart roundtrip N -> (u,v) -> N (100k)", r < 1e-12, f"worst {r:.2e}")

e, R = SC['e'], SC['R']
x = rng.uniform(-.9, .9, 50000); y = rng.uniform(-.9, .9, 50000)
Aq = x*x + y*y + e*e
disc = e*e*R*R - (e*e - R*R)*(x*x + y*y)
m = disc >= 0
x, y, disc, Aq = x[m], y[m], disc[m], Aq[m]
t = (e*e - np.sqrt(disc))/Aq
lon = np.abs(np.arctan2(t*y/R, t*x/R) - np.arctan2(y, x)).max()
check(5, "longitude theorem: atan2(Ny,Nx) = atan2(y,x)", lon < 1e-12,
      f"worst {lon:.2e} over {m.sum()} hits")

Nz1 = e*(1 - t)/R
x2, y2 = -y, x                                    # rotate 90 deg: same r^2
Aq2 = x2*x2 + y2*y2 + e*e
disc2 = e*e*R*R - (e*e - R*R)*(x2*x2 + y2*y2)
t2 = (e*e - np.sqrt(disc2))/Aq2
rnz = np.abs(Nz1 - e*(1 - t2)/R).max()
check(5, "Nz is a function of r^2 alone (constant-v curves are circles)",
      rnz < 1e-12, f"worst {rnz:.2e}")

c = rng.random(100_000)
rr = np.abs(srgb_inv(srgb(c)) - c).max()
check(5, "sRGB encode/decode roundtrip (100k)", rr < 1e-9, f"worst {rr:.2e}")

NB = 16
gi, gj = np.meshgrid(np.arange(NB), np.arange(NB))
ring = (np.floor(np.hypot(gi - 7.5, gj - 7.5)) % 2) == 0
bmp = np.where(ring[..., None], [0.85, 0.62, 0.25], [0.16, 0.20, 0.34]).astype(float)
bmp[4, 4] = [0.80, 0.15, 0.15]
okm = tuple(bmp[4, 4]) == (0.80, 0.15, 0.15)
okr = tuple(bmp[0, 8]) == (0.16, 0.20, 0.34) and tuple(bmp[8, 8]) == (0.16, 0.20, 0.34)
check(5, "bitmap syncs with the simulator's JS (marker + ring corners)",
      okm and okr, "marker [0.80,0.15,0.15] at (4,4)")

d_noncomm = abs(float(srgb(1.0)) - 2*float(srgb(0.5)))
check(5, "E(a+b) != E(a)+E(b): gamma order is a theorem, not a taste",
      d_noncomm > 0.01, f"|E(1) - 2 E(0.5)| = {d_noncomm:.4f}")

uu, vv = np.meshgrid(np.linspace(0, 1, 512, endpoint=False),
                     np.linspace(0, 1, 512, endpoint=False))
def recon_err(k):
    f = np.sign(np.sin(np.pi*k*uu)*np.sin(np.pi*k*vv)) > 0
    fb = (np.sign(np.sin(np.pi*k*(np.arange(NB) + .5)/NB)[None, :]*
                  np.sin(np.pi*k*(np.arange(NB) + .5)/NB)[:, None]) > 0)
    rn = fb[np.floor(vv*NB).astype(int), np.floor(uu*NB).astype(int)]
    return np.abs(rn.astype(float) - f.astype(float)).mean()
e6, e8, e12 = recon_err(6), recon_err(8), recon_err(12)
check(5, "the bitmap Nyquist knee at k = NB/2 = 8",
      e8 < 0.06 and e12 > e8 + 0.05,
      f"reconstruction error: k=6 {e6:.3f}, k=8 {e8:.3f}, k=12 {e12:.3f}")
Expect: seven PASS lines — the chart trio at machine precision, the bitmap exact-class, gamma's constant printed, the knee bent at k=8.

🐍 Cell 6 — Chapter 6 battery: the box integral, the thin band, the order of operations, convergence

· · paste into colab.new → Run
banner(6, "antialiasing: the pixel is an integral")
rng = np.random.default_rng(11)

def Gn(v):
    v = float(v)
    m = math.floor(v); f = 1 - abs(1 - 2*(v - m))
    return m + (f*f if m % 2 == 0 else 1 - f*f)
def checker_box(x0, y0, h):
    return (Gn(x0 + h) - Gn(x0))*(Gn(y0 + h) - Gn(y0))/(h*h)

worst = 0.0
for _ in range(300):
    x0, y0 = rng.uniform(-20, 20, 2); h = float(rng.uniform(0.01, 2))
    exact = checker_box(x0, y0, h)
    K = 64
    xs = x0 + (np.arange(K) + .5)/K*h
    ys = y0 + (np.arange(K) + .5)/K*h
    quad = np.sign(np.sin(np.pi*xs)[:, None]*np.sin(np.pi*ys)[None, :]).mean()
    worst = max(worst, abs(exact - quad))
check(6, "the checker box integral == 64x64 quadrature (300 random pixels)",
      worst < 2e-3, f"worst {worst:.2e}  [statistical class: scan resolution]")

img6, code6 = render(128)
edge = np.zeros(code6.shape, bool)
edge[1:, :] |= code6[1:, :] != code6[:-1, :]
edge[:, 1:] |= code6[:, 1:] != code6[:, :-1]
frac = edge.mean()
check(6, "the coverage band is thin (the expensive machinery touches ~2%)",
      frac < 0.06, f"{100*frac:.2f}% of frame")

c = np.linspace(0, 1, 400)
order_err = np.abs(srgb(0.5*(c + c[::-1])) - 0.5*(srgb(c) + srgb(c[::-1]))).max()
check(6, "filter and encode do not commute", order_err > 0.01,
      f"max |encode(avg) - avg(encode)| = {order_err:.4f}")

half = 1.25*SC['e']*SC['R']/math.sqrt(SC['e']**2 - SC['R']**2)
h = 2*half/64
x0, y0 = -1.3, -0.75*half
def floor_cov_num(xa, ya, M):
    acc = 0.0
    for sj in range(M):
        for si in range(M):
            xs_ = xa + (si + .5)/M*h; ys_ = ya + (sj + .5)/M*h
            acc += 1.0 if math.sin(math.pi*xs_)*math.sin(math.pi*ys_) >= 0 else 0.0
    return acc/(M*M)
errs = []
for M in (1, 2, 4, 8):
    es = []
    for j in range(8):
        for i in range(8):
            xa, ya = x0 + i*h, y0 + j*h
            es.append(abs(floor_cov_num(xa, ya, M) - (1 + checker_box(xa, ya, h))/2))
    errs.append(float(np.mean(es)))
mono = all(errs[k] >= errs[k + 1] - 1e-12 for k in range(3)) and errs[0] > errs[3]
check(6, "supersampling converges to the analytic floor (M = 1, 2, 4, 8)",
      mono, "mean coverage error: " + ", ".join(f"{e_:.4f}" for e_ in errs))
Expect: the quadrature check with its confessed tolerance; the band thin; the order error nonzero; four decreasing means.

🐍 Cell 7 — Chapter 7 battery: the anchor, the regression, the conic, and the two currencies

· · paste into colab.new → Run
banner(7, "area lights: penumbra is an integral")

d, rho = 6.0, 0.45
exact = 2*d*(math.sqrt(d*d + rho*rho) - d)/rho**2
K = 160_000
rr = (np.arange(K) + .5)/K*rho
quad = float(np.mean(d/np.sqrt(d*d + rr*rr)))
check(7, "on-axis anchor: mu-bar in closed form (160k quadrature)",
      abs(exact - quad) < 1e-6, f"exact {exact:.8f}, quad {quad:.8f}")

lims = [2*d*(math.sqrt(d*d + r2*r2) - d)/r2**2 for r2 in (0.45, 0.05, 1e-3)]
check(7, "mu-bar -> 1 as rho -> 0 (the point-light regression)",
      abs(lims[-1] - 1) < 1e-3, ", ".join(f"{l_:.8f}" for l_ in lims))

def disk_basis(L):
    w = -np.array(L, float); w /= np.linalg.norm(w)
    up = np.array([0., 1., 0.]) if abs(w[1]) < 0.9 else np.array([1., 0., 0.])
    e1 = np.cross(w, up); e1 /= np.linalg.norm(e1)
    return e1, np.cross(w, e1)

def conic_det(Pv):
    C = SC['Cb']; Rb = SC['R2']; Q0 = SC['lights'][0][0]
    e1, e2 = disk_basis(Q0)
    def dsc(uu, vv):
        Q = Q0 + uu*e1 + vv*e2
        seg = Q - Pv
        a2 = float(seg @ seg); b2 = float((Pv - C) @ seg)
        c2 = float((Pv - C) @ (Pv - C)) - Rb**2
        return b2*b2 - a2*c2
    A11 = (dsc(1, 0) + dsc(-1, 0) - 2*dsc(0, 0))/2
    A22 = (dsc(0, 1) + dsc(0, -1) - 2*dsc(0, 0))/2
    A12 = (dsc(1, 1) - dsc(1, -1) - dsc(-1, 1) + dsc(-1, -1))/8
    return A11*A22 - A12*A12
dets = [conic_det(np.array([tt, 0.0, math.sqrt(max(1e-9, 1 - tt*tt))]))
        for tt in np.linspace(-.9, .9, 7)]
check(7, "the penumbra boundary on the source is a conic, classified per point",
      all(np.isfinite(dets)) and all(d_ != 0 for d_ in dets),
      "dets: " + " ".join(f"{d_:+.3f}" for d_ in dets))

def litfrac(Pv, K, blockers, seed=None):
    L = SC['lights'][0][0]
    e1, e2 = disk_basis(L)
    rngL = np.random.default_rng(seed) if seed is not None else None
    n = 0
    for sj in range(K):
        for si in range(K):
            uu = (si + (rngL.random() if rngL else .5))/K
            vv = (sj + (rngL.random() if rngL else .5))/K
            r = 0.45*math.sqrt(uu); th = 2*math.pi*vv
            Q = L + r*(math.cos(th)*e1 + math.sin(th)*e2)
            n += 0 if any(shadowed_roots(Pv, Q, C, Rr) for C, Rr in blockers) else 1
    return n/(K*K)

cands = []
for tt in np.linspace(-0.9, 0.9, 25):
    Pv = np.array([tt, 0.0, math.sqrt(max(1e-9, 1 - tt*tt))])
    f64 = litfrac(Pv, 32, [(SC['Cb'], SC['R2'])])
    if 0.02 < f64 < 0.98:
        cands.append((Pv, f64, [(SC['Cb'], SC['R2'])]))
    if len(cands) >= 5: break
if not cands:
    for xx in np.linspace(-1, 1, 25):
        Pv = np.array([xx, -0.2, 0.0])
        f64 = litfrac(Pv, 32, [(np.zeros(3), SC['R']), (SC['Cb'], SC['R2'])])
        if 0.02 < f64 < 0.98:
            cands.append((Pv, f64, [(np.zeros(3), SC['R']), (SC['Cb'], SC['R2'])]))
        if len(cands) >= 5: break
check(7, "penumbral test points exist on the frame", len(cands) > 0,
      f"{len(cands)} found")
if cands:
    Pv, _, bl = cands[0]
    n7 = litfrac(Pv, 4, bl, seed=7); n99 = litfrac(Pv, 4, bl, seed=99)
    s1 = litfrac(Pv, 4, bl); s2 = litfrac(Pv, 4, bl)
    check(7, "numeric estimate has a seed (variance)", n7 != n99,
          f"seed 7: {n7:.4f} vs seed 99: {n99:.4f}")
    check(7, "symbolic estimate is deterministic (bias only)", s1 == s2,
          f"{s1:.4f} == {s2:.4f}  [exact class]")
    e2 = np.mean([abs(litfrac(P_, 2, b_) - f_) for P_, f_, b_ in cands])
    e16 = np.mean([abs(litfrac(P_, 16, b_) - f_) for P_, f_, b_ in cands])
    check(7, "strata converge toward the fine truth (bias shrinks)",
          e16 < 0.05 and e16 <= e2 + 0.02,
          f"mean |err|: K=2 {e2:.4f}, K=16 {e16:.4f}")
Expect: the anchor to 1e-6, the regression approaching 1, the conic dets printed, and the two currencies demonstrated on a penumbral point.

🐍 Cell 8 — Chapter 8 battery: weights, normalizer, partition of unity, concentration

· · paste into colab.new → Run
banner(8, "glossy reflection: the lobe integral")

n_g, thm = 16, math.pi/3
worst = 0.0
for k in range(8):
    ta, tb = k/8*thm, (k + 1)/8*thm
    wE = (math.cos(ta)**(n_g + 1) - math.cos(tb)**(n_g + 1))/(n_g + 1)
    tt = np.linspace(ta, tb, 4000)
    wQ = np.trapezoid(np.cos(tt)**n_g*np.sin(tt), tt)
    worst = max(worst, abs(wE - wQ))
check(8, "exact Lambert stratum weights vs quadrature", worst < 1e-8,
      f"worst {worst:.2e}  [statistical class]")

Wex = (1 - math.cos(thm)**(n_g + 1))/(n_g + 1)
tt = np.linspace(0, thm, 200_000)
WQ = np.trapezoid(np.cos(tt)**n_g*np.sin(tt), tt)
check(8, "the lobe normalizer in closed form", abs(Wex - WQ) < 1e-8,
      f"exact {Wex:.10e}, quad {WQ:.10e}")

sw = sum((math.cos(k/8*thm)**(n_g + 1) - math.cos((k + 1)/8*thm)**(n_g + 1))/(n_g + 1)
         for k in range(8))
check(8, "constant integrand integrates exactly (partition of unity)",
      abs(sw - Wex) < 1e-14, f"|sum w_k - W| = {abs(sw - Wex):.2e}  [machine class]")

n_big = 4096
tot = sum((math.cos(k/8*thm)**(n_big + 1) - math.cos((k + 1)/8*thm)**(n_big + 1))/(n_big + 1)
          for k in range(8))
first = (1 - math.cos(thm/8)**(n_big + 1))/(n_big + 1)
check(8, "mirror-grade lobe concentrates at theta = 0 (n = 4096)",
      first/tot > 0.9, f"{100*first/tot:.2f}% of weight in the first ring")
Expect: four PASS lines — the weights and normalizer to 1e-8, the partition exact, and >90% first-ring concentration.

🐍 Cell 9 — Chapter 9 battery: the U-curve, the closed-form derivative, the regression

· · paste into colab.new → Run
banner(9, "bump mapping: the derivative is the texture")

ku, kv = SC['ku'], SC['kv']
def h_sine(u, v): return math.sin(math.pi*ku*u)*math.sin(math.pi*kv*v)
def grad_sine(u, v):
    return (math.pi*ku*math.cos(math.pi*ku*u)*math.sin(math.pi*kv*v),
            math.pi*kv*math.sin(math.pi*ku*u)*math.cos(math.pi*kv*v))
u0, v0 = 0.31, 0.47
hu_ex = math.pi*ku*math.cos(math.pi*ku*u0)*math.sin(math.pi*kv*v0)
errs = {}
for ep in (1e-1, 1e-2, 1e-4, 1e-6, 1e-8, 1e-9):
    errs[ep] = abs((h_sine(u0 + ep, v0) - h_sine(u0 - ep, v0))/(2*ep) - hu_ex)
slope = errs[1e-1] > errs[1e-2] > errs[1e-4]
cliff = errs[1e-9] > errs[1e-6]
check(9, "the FD U-curve: truncation slope AND roundoff cliff",
      slope and cliff,
      "  ".join(f"{e_:.0e}: {errs[e_]:.1e}" for e_ in errs))

rng = np.random.default_rng(5)
N = rng.normal(size=(5000, 3)); N /= np.linalg.norm(N, axis=1, keepdims=True)
g = rng.normal(size=(5000, 3))*4
b0 = 0.35
M = N + b0*g
ml = np.linalg.norm(M, axis=1, keepdims=True)
Np = M/ml
dot = np.einsum('ij,ij->i', Np, g)[:, None]
ex = (g - Np*dot)/ml
db = 1e-7
M2 = N + (b0 + db)*g
fd = (M2/np.linalg.norm(M2, axis=1, keepdims=True) - Np)/db
worst = np.abs(ex - fd).max()
check(9, "dN'/dbeta = (I - N'N')g/|M| vs finite differences (5000 trials)",
      worst < 1e-6, f"worst {worst:.2e}  [machine class]")

img_b0, _ = render(64, beta=0.0, grad=grad_sine)
img_plain, _ = render(64)
img_b15, _ = render(64, beta=1.5, grad=grad_sine)
d0 = np.abs(img_b0 - img_plain).max()
d15 = np.abs(img_b15 - img_plain).mean()
check(9, "beta -> 0 regression: the bumped render returns the plain render",
      d0 < 1e-12, f"max {d0:.2e}")
check(9, "sanity: the bump is not a no-op at beta = 1.5", d15 > 1e-3,
      f"mean |delta| {d15:.3f}")
Expect: the U-curve with both slopes in the printed numbers; the derivative to 1e-6; the regression to 1e-12; the sanity gap large.

🐍 Cell 10 — Chapter 10 battery: switching times, the bokeh law, regressions, and sympy's certificates

· · paste into colab.new → Run
banner(10, "motion & depth: the camera gains a clock and a pupil")

Cb = SC['Cb']; vel = np.array([-0.55, -0.35, 0.35]); R2 = SC['R2']; e = SC['e']
rng = np.random.default_rng(777)

def coverage_interval(x, y, vv):
    D = np.array([x, y, -e]); Aq = float(D @ D)
    w0 = np.array([0., 0., e]) - Cb
    b0 = float(w0 @ D); vd = float(-vv @ D)
    al = vd*vd - Aq*float(vv @ vv)
    be = -2*b0*vd + 2*Aq*float(w0 @ vv)
    ga = b0*b0 - Aq*(float(w0 @ w0) - R2**2)
    dsc = be*be - 4*al*ga
    if dsc <= 0 or abs(al) < 1e-14: return 0.0
    sq = math.sqrt(dsc)
    lo, hi = sorted(((-be - sq)/(2*al), (-be + sq)/(2*al)))
    if al < 0:
        return max(0.0, min(1.0, hi) - max(0.0, lo))
    return 0.0

worst = 0.0
for _ in range(100):
    x, y = rng.uniform(-1, 1, 2)
    exact = coverage_interval(x, y, vel)
    ts = np.linspace(0, 1, 400)
    D = np.array([x, y, -e]); Aq = float(D @ D)
    Wm = (np.array([0., 0., e]) - Cb)[None, :] - ts[:, None]*vel[None, :]
    bm = Wm @ D
    scan = float(np.mean(bm*bm - Aq*(np.einsum('ij,ij->i', Wm, Wm) - R2**2) >= 0))
    worst = max(worst, abs(exact - scan))
check(10, "switching-time coverage (quadratic in t) vs 400-pt scans",
      worst < 6e-3, f"worst {worst:.2e}  [statistical class: scan resolution]")

f0, a0 = 1.0, 0.2
worst_rel = 0.0
for z0 in (0.8, 2.0, 3.0):
    law = a0*e*abs(f0 - z0)/((e - z0)*(e - f0))
    ps = []
    for q in np.linspace(-a0, a0, 21):
        mu = (e - f0)/(e - z0)
        Fx = (1 - mu)*q + mu*0.3
        ps.append(Fx*e/(e - f0))
    ps = np.array(ps)
    meas = np.abs(ps - ps[len(ps)//2]).max()
    worst_rel = max(worst_rel, abs(meas - law)/law)
check(10, "bokeh law: footprint = exact scaled aperture",
      worst_rel < 1e-9, f"worst relative error {worst_rel:.2e}  [machine class]")

mism = 0
for _ in range(200):
    x, y = rng.uniform(-1, 1, 2)
    static = (np.array([0., 0., e]) - Cb) @ np.array([x, y, -e])
    D = np.array([x, y, -e])
    w0 = np.array([0., 0., e]) - Cb
    static_hit = float(w0 @ D)**2 - float(D @ D)*(float(w0 @ w0) - R2**2) >= 0
    cov = coverage_interval(x, y, vel*1e-6)
    if abs(cov - (1.0 if static_hit else 0.0)) > 1e-3: mism += 1
check(10, "v -> 0 regression: coverage degenerates to the static mask",
      mism == 0, f"{mism} mismatches in 200 pixels")

import sympy as sp
tt = sp.symbols('t', real=True)
Pvx, Pvy, Pvz = sp.symbols('P_x P_y P_z', real=True)
Llx, Lly, Llz = sp.symbols('L_x L_y L_z', real=True)
Pv = sp.Matrix([Pvx, Pvy, Pvz]); Ll = sp.Matrix([Llx, Lly, Llz])
Ct2 = sp.Matrix([Cb[0] + vel[0]*tt, Cb[1] + vel[1]*tt, Cb[2] + vel[2]*tt])
seg = Ll - Pv
a2e = seg.dot(seg); b2e = (Pv - Ct2).dot(seg)
c2e = (Pv - Ct2).dot(Pv - Ct2) - R2**2
sdisc = sp.expand(b2e**2 - a2e*c2e)
deg = sp.Poly(sdisc, tt).degree()
check(10, "the moving shadow mask is quartic in t (sympy certifies)",
      deg == 4, f"degree = {deg}")

xv, yv = sp.symbols('x y', real=True)
Dxy = sp.Matrix([xv, yv, -e])
E0v = sp.Matrix([0, 0, e])
Ctv = sp.Matrix([Cb[0] + vel[0]*tt, Cb[1] + vel[1]*tt, Cb[2] + vel[2]*tt])
Dt2 = sp.expand(((E0v - Ctv).dot(Dxy))**2
                - Dxy.dot(Dxy)*((E0v - Ctv).dot(E0v - Ctv) - R2**2))
pt = sp.Poly(Dt2, tt)
env = sp.expand(pt.nth(1)**2 - 4*pt.nth(2)*pt.nth(0))
degxy = sp.Poly(env, xv, yv).total_degree()
check(10, "the swept silhouette is algebraic, degree <= 4",
      2 <= degxy <= 4, f"total degree = {degxy}")
Expect: the scan check with its confessed tolerance; the bokeh law at 1e-9; the regression clean; sympy printing degree 4 twice.

🐍 Cell 11 — the tally: the archive counts itself

· · paste into colab.new → Run
print("\n" + "="*64)
print("THE MACHINE-CHECKED ARCHIVE — TALLY")
print("="*64)
n_pass = sum(1 for r in RESULTS if r[2])
by_ch = {}
for ch, name, ok, detail in RESULTS:
    by_ch.setdefault(ch, [0, 0])
    by_ch[ch][0] += ok; by_ch[ch][1] += 1
for ch in sorted(by_ch, key=str):
    p, t = by_ch[ch]
    print(f"  chapter {str(ch):>2}: {p}/{t}")
print(f"\n  ARCHIVE TOTAL: {n_pass}/{len(RESULTS)} checks pass")
fails = [r for r in RESULTS if not r[2]]
if fails:
    print("\n  FAILURES (investigate the environment or the edit, not the theorem):")
    for ch, name, _, detail in fails:
        print(f"    ch{ch}: {name}  ({detail})")
else:
    print("\n  The archive agrees with the book on every line it was asked to check.")
    print("  The proofs live in the chapters; the seeds are fixed; the tolerances")
    print("  are on the record. Rerun tomorrow: the contract is kept by re-running it.")
print("\nFINIS SERIEI.")
Expect: per-chapter tallies summing to 44/44, then FINIS SERIEI.

§3 · What the archive cannot certify

§4 · The book closes

Ten chapters and two appendices: from one sphere and one light to a five-dimensional product integral with a 44-check insurance policy. The discipline never changed — derive the anchor, run the regression, report the error in both currencies, show the work. What the expression could be, it was; where it could not be, the book said so at the boundary, with measurements. The renderer's decision tree became mathematics; the mathematics, in the end, became a suite that reruns itself. Both traditions got what they came for.

❦
Set by the Press · the archive is the book's answer to bit rot.
FINIS SERIEI — the checks, being re-runnable, continue without us.
Section ids are stable (s1…s4, lab, cell0…cell11) — cite the id when requesting revisions.