Symbolic Rendering

Chapter 4 — Reflection, or Recursion Becomes Composition

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
Python-ready for Jupyter Notebooks in Google Colab.
Requires an internet connection for mathematical typesetting (MathJax).
❦
How to use this file. §1 what changes; §2 the scene; §3 the derivation — the reflected quadratic as literal substitution, the mirror-flip trick, and why the Shadow Lemma survives composition verbatim. §4 is the live simulator: numeric ray casting with explicit roots meets the symbolic composed expression; acne and Moiré are demonstrated, not described. §5 is the laboratory: Cells 1 and 3 replace Chapter 3's (the chapters are cumulative); Cells 2, 4–8 are new. §6–§9 as ever.
☰  contents
  1. What changes in Chapter 4
  2. The scene (Chapter-4 assumptions)
  3. The math, once, carefully
  4. Live simulator
  5. The Colab laboratory (Python / sympy)
  6. The honest ledger — Chapter-4 deltas
  7. Pitfalls gallery, continued
  8. A brief history, continued
  9. Roadmap

§1 · What changes in Chapter 4

All of Chapter 3 stands: two spheres, two lights, sign-logic shadows, the two-number ledger. Three novelties:

  1. A third object that is not a sphere. The plane $z=0$ — mirror-floored. It has no quadratic silhouette in the view window; its "hit test" is a single division. It is also, provably, the image plane of Chapter 1: the screen the previous chapters rendered on is now part of the scene.
  2. One bounce of perfect mirror reflection. Primary surfaces shade as before plus a mirror-weighted contribution of whatever the reflected ray sees: the floor, the other sphere, or the void.
  3. The ε-hack, at last. A reflected ray starts on the surface that spawned it, so naive re-intersection returns $t' = 0$ — self-intersect "acne". Numeric rendering's classic cure is the origin fudge. Symbolic rendering's cure is an exact strict test ($t' > 0$, no fudge). Both pipelines offer both versions, so the failure can be switched on and watched — while their mutual ledger stays clean.

§2 · The scene (Chapter-4 assumptions)

Chapter 3's scene, plus: the plane $z=0$ with mirror weight $k_{r,f} \in [0,1]$; spheres with mirror weights $k_{r,A}, k_{r,B}$ (default 0.15). One bounce only: a reflected ray that meets a mirror does not reflect again (the second bounce is named and priced in §3). The floor's shading is a procedural checker — previewed here, derived properly in Chapter 5; both pipelines evaluate the identical closed-form parity function, so it introduces no approximation asymmetry. Non-containment, non-interpenetration, and the convex-self-shadow exclusion all carry over. Lights above the floor; camera behind it.

§3 · The math, once, carefully

The plane hit, closed form. Ray $E + tD$, $D=(x,y,-e)$: the floor is met when $e - te = 0$, i.e.

$$t_f = 1 \quad \text{for } e \ne 0, \qquad P_f = (x,\,y,\,0), \qquad N_f = (0,0,1).$$

The image plane of Chapter 1 is the mirror: every pixel's ray crosses the floor at $t=1$, exactly once, no quadratic. The mirror's hit condition is $t_f \lt t_\text{nearest sphere}$ — the floor shows through wherever the spheres are farther than the screen.

The reflected ray is substitution. At any primary hit $P$ with normal $N$:

$$D' = D - 2\,(D\cdot N)\,N, \qquad \text{reflected ray: } P + t'\,D',\quad t' \gt 0.$$

The numeric renderer computes $D'$ at runtime. The symbolic renderer substitutes the already-derived expressions for $P$ and $N$ into the definition of $D'$ — the geometry is composition of expressions, and the reflected image is the same derived object, one substitution deeper. This is Cell 2's composition census: the reflected discriminant of Chapter 3's Cell 8 was a monster; measured, hoisted by cse, it becomes an evaluable object with a counted cost.

The sphere re-hit. Against sphere $(C, R_s)$ the reflected quadratic is $a' t'^2 + 2b' t' + c' = 0$ with $a' = D'\cdot D'$, $b' = (P-C)\cdot D'$, $c' = |P-C|^2 - R_s^2$, and near root $(-b' - \sqrt{b'^2 - a'c'})/a'$. Against the floor: $t'_f = -P_z/D'_z$ when $D'_z \lt 0$, else the ray ascends to the void. The reflected color is then the chamber expression evaluated at the reflected hit — floor checker with Blinn shading and sign-logic shadows, or the other sphere's own chamber expression, one level deep (no second bounce).

The Shadow Lemma survives composition. A reflected shadow feeler asks the same ∃-root-in-$(0,1)$ question as Chapter 3, merely from a different origin. Non-containment is unchanged (reflected hits lie on sphere surfaces, outside the blocker), so the same three signs decide it — the lemma did not know or care whether its segment began on a primary or a reflected hit. One theorem, two generations of rays.

The exact strict test (no ε). The reflected origin lies exactly on the spawning surface: $c' = 0$, and the quadratic factors $t'(a't' + 2b') = 0$ with roots $0$ and $-2b'/a'$. Against the same sphere, $b' = (P-C)\cdot D' = R_s\,(N\cdot D') \gt 0$ (reflection leaves the surface: $D\cdot N \lt 0 \Rightarrow D'\cdot N = -D\cdot N \gt 0$), so the other root is negative and the strict test $t' \gt 0$ excludes exactly the degenerate root $0$ — no offset, no fudge, no acne. The numeric ε-hack shifts the origin and perturbs everything downstream (a measured, visible speckle in Cell 6); the symbolic strict test perturbs nothing.

Three gifts of Chapter 4:

  1. The mirror flip is a change of variables. Reflecting the scene across $z=0$ ($z \mapsto -z$ for centers and lights) makes the reflected ray a primary ray in the mirrored scene — Chapter 1's entire derivation re-applied under substitution. Cell 5 renders the mirrored spheres with the Chapter-1 pipeline and matches the bounce path to machine precision. A classic rasterizer trick, here a theorem about expression substitution.
  2. Acne is a theorem-shaped bug. The self-hit is the root $t'=0$ of a factoring quadratic; excluding it strictly is exact. The ε-hack is an approximate patch for a problem the algebra simply does not have.
  3. The checker floor is Chapter 6's warning label. The floor's parity function $\chi = \mathrm{sign}(\sin \pi x \cdot \sin \pi y)$ has frequency content at every scale — magnify the grazing band and the reflected checker aliases into Moiré. Sampling a symbolic function is not automatically exact: the box filter must be derived, not assumed.

Counterpoint: composition is swell. One bounce roughly triples the symbolic expression's op count; a second bounce compounds it past usability (Cell 8 measures the growth), and glossy reflection — an integral over reflected directions — has no closed form at all. Totality also deepens: the symbolic form evaluates the bounce for every surface pixel, hit or miss. The closed-form era ends here, on schedule: Chapter 5's textures are symbolic functions but Chapter 6's antialiasing is an integral, and integrals do not compose.

§4 · Live simulator

Two spheres over a mirror floor (the plane is the screen of Chapter 1). Numeric: runtime ifs and explicit roots. Symbolic: the composed expression with strict $t' > 0$. Acne toggle compares the ε-hack against the exact test in both pipelines; the ledger keeps its two numbers.

model 
resolution 
initializing…
NUMERIC — ray casting; bounce by runtime reflection
SYMBOLIC — the composed expression; reflection is substitution
|difference| ×10¹⁵ — black interior; white = branch flips at boundary curves
numeric: – symbolic: – interior max |ΔI|: – branch flips: –
hover the canvases — the probe shows pixel → world (x,y), the primary chamber, the reflected chamber, and both RGB triples.

Try this: (1) Set floor mirror k_rf to 0.9 and checker α to 0.45 — the spheres stand on a mirror; the bounce view isolates the reflection term alone. (2) Switch secondary rays to the ε-hack: both pipelines switch together, so the difference ledger stays black — the speckle the hack introduces is measured against the exact test in Cell 6, where the two formulations are compared. (3) Set magnify to 0.05: the grazing-band reflected checker aliases into Moiré — a symbolic function, under-sampled. Chapter 6's box filter is the cure, and it must be derived, not assumed.

§5 · The Colab laboratory (Python / sympy)

The chapters are cumulative: Cells 1 and 3 replace Chapter 3's (primary + reflected renderer, and the agreement census). Cells 2, 4–8 are new. Paste in order, run in order; three run paths per cell — 📋 Copy code, ⬇ notebook (File → Upload notebook in Colab), or paste into colab.new.

🐍 Cell 1 (REPLACES Ch.3 Cell 1) — the numerical baseline: ray casting with one bounce, ε-hack on offer

· · paste into colab.new → Run
import numpy as np
import matplotlib.pyplot as plt
import time

SC4 = dict(R=1.0, e=4.0,
           Cb=(0.5, 0.5, 1.3), R2=0.32, b_on=True, floor_on=True,
           lights=(((5., 5., 10.), (1., 1., 1.)),
                   ((-6., 2., 8.), (0.9, 0.65, 0.25))),
           ambient=(0.15, 0.15, 0.16),
           matA=(0.63, 0.19, 0.16), matB=(0.33, 0.44, 0.62),
           ka=0.1, kd=0.7, ks=0.5, shin=32.0, model='blinn',
           krA=0.15, krB=0.15, krf=0.55, alpha=0.25, selftest='strict')

def _checker(Px, Py, alpha):
    """Floor albedo: closed-form parity function, identical in both pipelines."""
    s = np.sign(np.sin(np.pi*Px) * np.sin(np.pi*Py))
    s = np.where(s == 0, 1, s)
    base = np.where(s > 0, 0.85, 0.15)
    return np.stack([base*(1-alpha), base*(1-alpha), base*(1+0*alpha)], axis=-1)

def _hit_sphere(Q, Dv, C, Rr, strict):
    """Nearest t > 0 (strict) or t > 1e-4 (hack); None if missed."""
    a = float(Dv @ Dv)
    b = float((Q - C) @ Dv)
    c = float((Q - C) @ (Q - C)) - Rr*Rr
    dsc = b*b - a*c
    if dsc < 0:
        return None
    sq = np.sqrt(dsc)
    lo = 0.0 if strict else 1e-4
    for t in sorted(((-b - sq)/a, (-b + sq)/a)):
        if t > lo:
            return t
    return None

def _lit(P, Lp, C, Rr, strict):
    """Shadow feeler by EXPLICIT roots (numeric formulation)."""
    seg = Lp - P
    a2 = float(seg @ seg)
    b2 = float((P - C) @ seg)
    c2 = float((P - C) @ (P - C)) - Rr*Rr
    dsc = b2*b2 - a2*c2
    if dsc >= 0:
        sq = np.sqrt(dsc)
        s1, s2 = (-b2 - sq)/a2, (-b2 + sq)/a2
        lo = 0.0 if strict else 1e-4
        if (lo < s1 < 1) or (lo < s2 < 1):
            return False
    return True

def _shade(P, N, v_hat, mat, blocker, SC, strict):
    col = SC['ka'] * np.array(SC['ambient'])
    code = 0
    for li, (pos, color) in enumerate(SC['lights']):
        if color is None:
            continue
        Lp = np.array(pos, dtype=float)
        if blocker is not None and not _lit(P, Lp, blocker[0], blocker[1], strict):
            continue
        w = Lp - P; w_hat = w / np.linalg.norm(w)
        mu = float(N @ w_hat)
        if SC['model'] == 'phong':
            r = 2.0*mu*N - w_hat
            spec = max(0.0, float(r @ v_hat)) ** SC['shin']
        else:
            h = w_hat + v_hat; h = h / np.linalg.norm(h)
            spec = max(0.0, float(N @ h)) ** SC['shin']
        col += np.array(color) * (SC['kd']*np.array(mat)*max(0.0, mu)
                                  + SC['ks']*spec)
        code |= 1 << li
    return col, code

def render_numeric_rfl(width, height, SC=SC4):
    """Chapter 4 numeric: primary = A | B | floor; one reflected bounce.
       strict: exact t' > 0 self-exclusion.  hack: classic epsilon offset."""
    strict = SC['selftest'] == 'strict'
    R, e = SC['R'], SC['e']
    Cb = np.array(SC['Cb']); R2 = SC['R2']
    img = np.zeros((height, width, 3))
    pcode = np.zeros((height, width), dtype=np.uint8)   # primary chamber
    rcode = np.zeros((height, width), dtype=np.uint8)   # reflected chamber
    E = np.array([0.0, 0.0, e])
    half = 1.25*e*R/np.sqrt(e**2 - R**2)
    for j in range(height):
        y = half * (1.0 - 2.0*(j + 0.5)/height)
        for i in range(width):
            x = half * (2.0*(i + 0.5)/width - 1.0)
            D = np.array([x, y, -e])
            tA = _hit_sphere(E, D, np.zeros(3), R, True)
            tB = _hit_sphere(E, D, Cb, R2, True) if SC['b_on'] else None
            tf = 1.0 if SC['floor_on'] else None
            # ---- nearest surface: runtime ordering branches ----
            cand = [(t, k) for t, k in ((tA, 0), (tB, 1), (tf, 2)) if t is not None]
            if not cand:
                continue
            t, kind = min(cand)
            P = E + t*D
            if kind == 0:
                Cen, Rr, mat, kr = np.zeros(3), R, np.array(SC['matA']), SC['krA']
                block = (Cb, R2) if SC['b_on'] else None
            elif kind == 1:
                Cen, Rr, mat, kr = Cb, R2, np.array(SC['matB']), SC['krB']
                block = (np.zeros(3), R)
            else:
                Cen, Rr, mat, kr = None, None, _checker(P[0], P[1], SC['alpha']), SC['krf']
                block = None
            N = (P - Cen)/Rr if kind < 2 else np.array([0., 0., 1.])
            v_hat = (E - P)/np.linalg.norm(E - P)
            pcode[j, i] = kind + 1
            # ---- direct shading ----
            col, lc = _shade(P, N, v_hat, mat, block, SC, strict)
            pcode[j, i] |= lc << 2
            # ---- one bounce: runtime reflection ----
            Dp = D - 2.0*float(D @ N)*N
            Q = P if strict else P + 1e-4*Dp           # the epsilon hack
            tbA = _hit_sphere(Q, Dp, np.zeros(3), R, strict)
            tbB = _hit_sphere(Q, Dp, Cb, R2, strict) if SC['b_on'] else None
            tbf = -P[2]/Dp[2] if (SC['floor_on'] and Dp[2] < 0) else None
            bcand = [(t, k) for t, k in ((tbA, 0), (tbB, 1), (tbf, 2)) if t is not None]
            if bcand:
                tb, bk = min(bcand)
                Q2 = Q + tb*Dp
                if bk == 0:
                    BC, BR, bm = np.zeros(3), R, np.array(SC['matA'])
                    bb = (Cb, R2) if SC['b_on'] else None
                elif bk == 1:
                    BC, BR, bm = Cb, R2, np.array(SC['matB'])
                    bb = (np.zeros(3), R)
                else:
                    BC, BR, bm = None, None, _checker(Q2[0], Q2[1], SC['alpha'])
                    bb = None
                BN = (Q2 - BC)/BR if bk < 2 else np.array([0., 0., 1.])
                bv = (E - Q2)/np.linalg.norm(E - Q2)
                bcol, bc2 = _shade(Q2, BN, bv, bm, bb, SC, strict)
                col = col + kr*bcol
                rcode[j, i] = (bk + 1) | (bc2 << 2)
            img[j, i] = col
    return np.clip(img, 0, 1), pcode, rcode

t0 = time.perf_counter()
img_num, pcode_num, rcode_num = render_numeric_rfl(256, 256)
t1 = time.perf_counter()
print(f"numeric render (primary + bounce): {t1 - t0:.3f} s")
fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].imshow(img_num); axs[0].set_title("Numeric: spheres on a mirror floor")
# bounce-only view: re-render with kr=0 and subtract (numeric approximation view)
img_flat, _, _ = render_numeric_rfl(256, 256, {**SC4, 'krA': 0, 'krB': 0, 'krf': 0})
axs[1].imshow(np.clip(img_num - img_flat, 0, 1)*2)
axs[1].set_title("bounce contribution (×2)")
for ax in axs: ax.axis('off')
plt.show()
Expect: spheres standing on a dim checker mirror, reflections visible on both; a few seconds at 256². The bounce view isolates the reflected energy.

🐍 Cell 2 (NEW) — the composition census: substitution, measured; the mirror flip as change of variables

· · paste into colab.new → Run
import sympy as sp

x, y = sp.symbols('x y', real=True)
R, e = sp.symbols('R e', positive=True)

# --- Primary geometry (Chapter 1, unchanged) -----------------------------------
Aq = x**2 + y**2 + e**2
disc = e**2*R**2 - (e**2 - R**2)*(x**2 + y**2)
tA = (e**2 - sp.sqrt(disc)) / Aq
P = sp.Matrix([tA*x, tA*y, e*(1 - tA)])
N = P / R
D = sp.Matrix([x, y, -e])

# --- Reflection IS substitution -------------------------------------------------
Dp = D - 2*D.dot(N)*N                 # reflected direction: P, N substituted in
print(f"reflected direction D': {sp.count_ops(Dp)} ops per component")

# --- Reflected quadratic against sphere B (Chapter 3 Cell 8's monster) ----------
bx, by, bz, R2 = sp.symbols('b_x b_y b_z R_2', real=True)
Cb = sp.Matrix([bx, by, bz])
ar = Dp.dot(Dp)
br = (P - Cb).dot(Dp)
cr = (P - Cb).dot(P - Cb) - R2**2
disc_r = sp.expand(br**2 - ar*cr)
ops_raw = sp.count_ops(disc_r)
print(f"reflected discriminant, raw:        {ops_raw} ops   <- the monster")

subs, reduced = sp.cse([disc_r])
ops_cse = sum(sp.count_ops(r) for _, r in subs) + sp.count_ops(reduced[0])
print(f"reflected discriminant, after cse:  {ops_cse} ops   "
      f"({len(subs)} shared subexpressions hoisted)")
print(f"compression factor: {ops_raw/ops_cse:.1f}x")
print()
print("The monster was never mystical: it is a composition with shared parts,")
print("and cse finds the parts. 'How much swell can we claw back' has an answer.")

# --- The mirror flip: reflection across z=0 is a change of variables ------------
E_m = sp.Matrix([0, 0, -e])           # mirrored eye
Cb_m = sp.Matrix([bx, by, -bz])       # mirrored center
L1_m = sp.Matrix(sp.symbols('l_x l_y l_z', real=True)); L1_m[2] = -L1_m[2]
print("mirrored scene: E -> (0,0,-e), C_b -> (b_x, b_y, -b_z), L -> (l_x, l_y, -l_z)")
print("A primary ray from E_m through (x, y, 0) IS the reflected ray.")
print("Chapter 1's derivation re-applies under substitution -- verified")
print("numerically against the bounce path in Cell 5.")
Expect: the reflected discriminant's raw op count in the hundreds-to-thousands, then cse clawing back a large factor — the composition bomb has a defuser, and its strength is now a number.

🐍 Cell 3 (REPLACES Ch.3 Cell 3) — the symbolic composed renderer; agreement census with bounce chambers

· · paste into colab.new → Run
def grid(W, H, R=1.0, e=4.0):
    half = 1.25*e*R/np.sqrt(e*e - R*R)
    xs = np.linspace(-half + half/W,  half - half/W,  W)
    ys = np.linspace( half - half/H, -half + half/H,  H)
    return np.meshgrid(xs, ys)

def render_symbolic_rfl(X, Y, SC=SC4, regularized=True):
    """Chapter 4 symbolic: primary chambers + ONE composed bounce.
       Reflection is substitution; shadows by the Shadow Lemma's signs;
       self-exclusion by the exact STRICT test t' > 0 (no epsilon anywhere)."""
    R, e = SC['R'], SC['e']
    bx, by, bz = SC['Cb']; R2 = SC['R2']
    sq = (lambda v: np.sqrt(np.maximum(v, 0))) if regularized else np.sqrt
    Aq = X**2 + Y**2 + e**2
    # ---- primary: spheres ----
    discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
    tA = (e*e - sq(discA)) / Aq
    bB = -bx*X - by*Y - (e - bz)*e
    cB = bx**2 + by**2 + (e - bz)**2 - R2**2
    discB = bB**2 - Aq*cB
    tB = (-bB - sq(discB)) / Aq
    onA = (discA >= 0) & (~SC['b_on'] | (discB < 0) | (tA <= tB))
    onB = SC['b_on'] & (discB >= 0) & ~onA
    # ---- primary: floor (the plane z=0 is hit at t=1) ----
    tmin_s = np.where(onA, tA, np.inf)
    tmin_s = np.where(onB, np.minimum(tmin_s, tB), tmin_s)
    onF = SC['floor_on'] & (1.0 < tmin_s)
    pcode = (onA*1 + onB*2 + onF*3).astype(np.uint8)
    v = np.sqrt(Aq); vx, vy, vz = -X/v, -Y/v, e/v
    img = np.zeros(X.shape + (3,))
    amb = SC['ka']*np.array(SC['ambient'])

    def shade(Px, Py, Pz, Nx, Ny, Nz, mat, blocker, vxx, vyy, vzz, code_out, base):
        col = np.zeros(X.shape + (3,)) + amb
        for li, (pos, color) in enumerate(SC['lights']):
            if color is None:
                continue
            wx, wy, wz = pos[0]-Px, pos[1]-Py, pos[2]-Pz
            a2 = wx*wx + wy*wy + wz*wz
            lit = np.ones(X.shape, dtype=bool)
            if blocker is not None:
                BC, BR = blocker
                qx, qy, qz = Px-BC[0], Py-BC[1], Pz-BC[2]
                b2 = qx*wx + qy*wy + qz*wz
                c2 = qx*qx + qy*qy + qz*qz - BR**2
                dsc = b2*b2 - a2*c2
                lit = ~((dsc >= 0) & (b2 < 0) & (a2 + b2 > 0))    # Shadow Lemma
            wn = np.sqrt(a2)
            mu = (Nx*wx + Ny*wy + Nz*wz)/wn
            if SC['model'] == 'phong':
                wv = (wx*vxx + wy*vyy + wz*vzz)/wn
                nud = mu                                       # N.v for this ray
                spec = np.maximum(0, 2*mu*nud - wv) ** SC['shin']
            else:
                hx, hy, hz = wx/wn+vxx, wy/wn+vyy, wz/wn+vzz
                hn = np.sqrt(hx*hx + hy*hy + hz*hz)
                spec = np.maximum(0, (Nx*hx + Ny*hy + Nz*hz)/hn) ** SC['shin']
            add = SC['kd']*np.array(mat)*np.maximum(0, mu)[..., None] \
                  + SC['ks']*spec[..., None]
            col += np.array(color)*add*lit[..., None]
            code_out[...] = code_out | ((lit & base).astype(np.uint8) << (2 + li))
        return col

    rcode = np.zeros(X.shape, dtype=np.uint8)
    for kind, on in ((0, onA), (1, onB), (2, onF)):
        if kind == 0:
            t, Rr, C, mat, kr = tA, R, (0., 0., 0.), SC['matA'], SC['krA']
            blocker = (SC['Cb'], R2) if SC['b_on'] else None
        elif kind == 1:
            t, Rr, C, mat, kr = tB, R2, SC['Cb'], SC['matB'], SC['krB']
            blocker = ((0., 0., 0.), R)
        else:
            t, kr = np.ones(X.shape), SC['krf']
            C, mat, blocker = None, None, None
        Px, Py, Pz = t*X, t*Y, e*(1 - t)
        if kind < 2:
            Nx, Ny, Nz = (Px-C[0])/Rr, (Py-C[1])/Rr, (Pz-C[2])/Rr
            base_alb = np.array(mat)
        else:
            Nx, Ny, Nz = np.zeros(X.shape), np.zeros(X.shape), np.ones(X.shape)
            base_alb = _checker(Px, Py, SC['alpha'])
            col_f = np.zeros(X.shape + (3,)) + amb
            # floor shading per light (blockers: BOTH spheres)
            for li, (pos, color) in enumerate(SC['lights']):
                if color is None:
                    continue
                wx, wy, wz = pos[0]-Px, pos[1]-Py, pos[2]-Pz
                a2 = wx*wx + wy*wy + wz*wz
                lit = np.ones(X.shape, dtype=bool)
                for BC, BR in (((0., 0., 0.), R), (SC['Cb'], R2) if SC['b_on'] else None):
                    if BC is None: continue
                    qx, qy, qz = Px-BC[0], Py-BC[1], Pz-BC[2]
                    b2 = qx*wx + qy*wy + qz*wz
                    c2 = qx*qx + qy*qy + qz*qz - BR**2
                    dsc = b2*b2 - a2*c2
                    lit &= ~((dsc >= 0) & (b2 < 0) & (a2 + b2 > 0))
                wn = np.sqrt(a2)
                mu = wz/wn                                     # N = +z
                hx, hy, hz = wx/wn+vx, wy/wn+vy, wz/wn+vz
                hn = np.sqrt(hx*hx + hy*hy + hz*hz)
                spec = np.maximum(0, hz/hn) ** SC['shin'] if SC['model'] == 'blinn' \
                       else np.maximum(0, 2*mu*(wz/wn) - (wx*vx+wy*vy+wz*vz)/wn) ** SC['shin']
                col_f += np.array(color) * (SC['kd']*base_alb*np.maximum(0, mu)[..., None]
                                            + SC['ks']*spec[..., None]) * lit[..., None]
                pcode[...] = pcode | ((lit & onF).astype(np.uint8) << (2 + li))
            img += col_f * onF[..., None]
        if kind < 2:
            pcode_loc = np.zeros(X.shape, dtype=np.uint8)
            col_s = shade(Px, Py, Pz, Nx, Ny, Nz, mat, blocker, vx, vy, vz,
                          pcode_loc, on)
            pcode[...] = pcode | np.where(on[..., None].any(-1) if False else on, pcode_loc, 0)
            img += col_s * on[..., None]
        # ---- ONE BOUNCE: reflection by substitution; strict t' > 0 ----
        dn = (X*Nx + Y*Ny + (-e)*Nz) if kind < 2 else (-e)*1.0
        if kind < 2:
            Dn = X*Nx + Y*Ny - e*Nz
        else:
            Dn = -e*Nz
        Dpx, Dpy, Dpz = X - 2*Dn*Nx, Y - 2*Dn*Ny, -e - 2*Dn*Nz
        # re-hit spheres (strict excludes the degenerate root t'=0 on the spawner)
        def rehit(C2, Rr2):
            aa = Dpx**2 + Dpy**2 + Dpz**2
            bb = (Px-C2[0])*Dpx + (Py-C2[1])*Dpy + (Pz-C2[2])*Dpz
            cc = (Px-C2[0])**2 + (Py-C2[1])**2 + (Pz-C2[2])**2 - Rr2**2
            dd = bb**2 - aa*cc
            rt = np.where(dd >= 0, (-bb - sq(dd))/np.maximum(aa, 1e-30), np.inf)
            rt = np.where(rt > 0, rt, np.inf)                    # STRICT: exact
            return rt
        rA = rehit((0., 0., 0.), R)
        rB = rehit(SC['Cb'], R2) if SC['b_on'] else np.full(X.shape, np.inf)
        rF = np.where(SC['floor_on'] & (Dpz < 0), -Pz/np.where(Dpz < 0, Dpz, -1e-30),
                      np.inf)
        tmin_r = np.minimum(np.minimum(rA, rB), rF)
        brA = on & (rA == tmin_r) & np.isfinite(rA)
        brB = on & (rB == tmin_r) & np.isfinite(rB)
        brF = on & (rF == tmin_r) & np.isfinite(rF)
        # bounce shading at the reflected hit (no second bounce)
        for bk, bon in ((0, brA), (1, brB), (2, brF)):
            tb = np.where(bon, tmin_r, 0.0)
            Qx, Qy, Qz = Px + tb*Dpx, Py + tb*Dpy, Pz + tb*Dpz
            if bk == 0:
                BR_, BC_, bm = R, (0., 0., 0.), SC['matA']
                bb = (SC['Cb'], R2) if SC['b_on'] else None
            elif bk == 1:
                BR_, BC_, bm = R2, SC['Cb'], SC['matB']
                bb = ((0., 0., 0.), R)
            else:
                BR_, BC_, bb = None, None, None
                bm = None
            if bk < 2:
                BNx, BNy, BNz = (Qx-BC_[0])/BR_, (Qy-BC_[1])/BR_, (Qz-BC_[2])/BR_
                alb = np.array(bm)
            else:
                BNx, BNy, BNz = np.zeros(X.shape), np.zeros(X.shape), np.ones(X.shape)
                alb = _checker(Qx, Qy, SC['alpha'])
            vq = np.sqrt((0-Qx)**2 + (0-Qy)**2 + (e-Qz)**2)
            bvx, bvy, bvz = -Qx/vq, -Qy/vq, (e-Qz)/vq
            rcol = np.zeros(X.shape + (3,)) + amb
            for li, (pos, color) in enumerate(SC['lights']):
                if color is None: continue
                wx, wy, wz = pos[0]-Qx, pos[1]-Qy, pos[2]-Qz
                a2 = wx*wx + wy*wy + wz*wz
                lit = np.ones(X.shape, dtype=bool)
                blockers = [((0.,0.,0.), R)] + ([(SC['Cb'], R2)] if SC['b_on'] else [])
                for BC3, BR3 in blockers:
                    qx, qy, qz = Qx-BC3[0], Qy-BC3[1], Qz-BC3[2]
                    b2 = qx*wx + qy*wy + qz*wz
                    c2 = qx*qx + qy*qy + qz*qz - BR3**2
                    dsc = b2*b2 - a2*c2
                    lit &= ~((dsc >= 0) & (b2 < 0) & (a2 + b2 > 0))
                wn = np.sqrt(a2)
                mu = (BNx*wx + BNy*wy + BNz*wz)/wn
                hx, hy, hz = wx/wn+bvx, wy/wn+bvy, wz/wn+bvz
                hn = np.sqrt(hx*hx + hy*hy + hz*hz)
                spec = np.maximum(0, (BNx*hx+BNy*hy+BNz*hz)/hn) ** SC['shin'] \
                       if SC['model'] == 'blinn' else \
                       np.maximum(0, 2*mu*((BNx*bvx+BNy*bvy+BNz*bvz)) -
                                  (wx*bvx+wy*bvy+wz*bvz)/wn) ** SC['shin']
                rcol += np.array(color)*(SC['kd']*alb*np.maximum(0, mu)[..., None]
                        + SC['ks']*spec[..., None])*lit[..., None]
            img += kr * rcol * bon[..., None]
    return np.clip(np.nan_to_num(img), 0, 1), pcode, rcode

X, Y = grid(256, 256)
t0 = time.perf_counter()
img_sym, pcode_sym, rcode_sym = render_symbolic_rfl(X, Y)
t1 = time.perf_counter()
print(f"symbolic render (primary + composed bounce): {t1 - t0:.4f} s")

same = pcode_num == pcode_sym
d = np.abs(img_num - img_sym)
print(f"interior max |num - sym|  = {d[same].max():.3e}   over {same.sum()} px")
print(f"branch flips (primary chambers): {(~same).sum()} px")

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(img_num); axs[0].set_title("Numeric")
axs[1].imshow(img_sym); axs[1].set_title("Symbolic (composed)")
ov = img_sym.copy(); ov[~same] = [1, 1, 1]
axs[2].imshow(ov); axs[2].set_title("flips in white")
for ax in axs: ax.axis('off')
plt.show()
Expect: interior max ~1e-14 over the primary chambers; a sparse flip census. The composed expression and the ray caster agree wherever the branch structure does.

🐍 Cell 4 (NEW) — honest timings: the totality tax deepens; three bounce policies priced

· · paste into colab.new → Run
print(f"{'grid':>9} {'numeric loop':>14} {'symbolic, no bounce':>21} {'symbolic + bounce':>19}")
for W in (64, 128, 256):
    Xg, Yg = grid(W, W)
    t0 = time.perf_counter(); render_numeric_rfl(W, W);                 t1 = time.perf_counter()
    render_symbolic_rfl(Xg, Yg, {**SC4, 'krA': 0, 'krB': 0, 'krf': 0}); t2 = time.perf_counter()
    render_symbolic_rfl(Xg, Yg);                                        t3 = time.perf_counter()
    print(f"{W:>4}x{W:<4} {t1-t0:>13.3f}s {t2-t1:>20.3f}s {t3-t2:>18.3f}s")

print("""
Reading the table:
 * The numeric loop early-exits everywhere: pixels that miss all geometry do
   nothing; pixels with kr = 0 skip the bounce. Branches are its privilege.
 * The symbolic form is TOTAL: every surface pixel evaluates the bounce path,
   hit or miss, weighted or not. The bounce roughly triples the work.
 * That tax buys the object: the composed expression can be differentiated,
   integrated, and re-parameterized (Cells 5-6). Speed was never the point;
   it is merely not the catastrophe the loop column suggests.
""")
Expect: loop ≫ symbolic; symbolic-with-bounce ≈ 2–3× symbolic-without. The totality tax, priced honestly.

🐍 Cell 5 (NEW) — superpowers: the mirror flip, verified; the virtual image is a theorem

· · paste into colab.new → Run
# The mirror-flip trick, made theorem: reflecting the SCENE across z=0 turns
# the reflected ray into a PRIMARY ray. Render the mirrored spheres with a
# primary-only pipeline and compare against the bounce path.

def render_mirror_scene(X, Y, SC=SC4):
    """Primary-style render of the scene mirrored across z=0:
       eye E' = (0,0,-e); sphere centers z -> -z; lights z -> -z.
       Viewed through the plane z=0 along the SAME pixel direction D=(x,y,-e),
       continued past the plane -- i.e., the reflected ray's journey, unfolded."""
    e, R = SC['e'], SC['R']
    bx, by, bz = SC['Cb']; R2 = SC['R2']
    Aq = X**2 + Y**2 + e**2
    # mirrored spheres (centers negated in z), ray from E'=(0,0,-e) through
    # (x, y, 0): direction D_m = (x, y, e). Note: continuing past the plane.
    sq = lambda v: np.sqrt(np.maximum(v, 0))
    def hit_m(C):
        Dx, Dy, Dz = X, Y, e
        Ex, Ey, Ez = 0.0, 0.0, -e
        aa = Dx**2 + Dy**2 + Dz**2
        bb = ((Ex-C[0])*Dx + (Ey-C[1])*Dy + (Ez-C[2])*Dz)
        cc = (Ex-C[0])**2 + (Ey-C[1])**2 + (Ez-C[2])**2 - C[3]**2
        dd = bb**2 - aa*cc
        tt = np.where(dd >= 0, (-bb - sq(dd))/aa, np.inf)
        return np.where(tt > 0, tt, np.inf)
    tAm = hit_m((0., 0., 0., R))
    tBm = hit_m((bx, by, -bz, R2)) if SC['b_on'] else np.full(X.shape, np.inf)
    return np.minimum(tAm, tBm)          # what the unfolded ray meets first

t_mirror = render_mirror_scene(X, Y)

# Compare against the symbolic bounce path's reflected sphere hits, restricted
# to pixels whose primary surface is the floor (the floor's reflection IS the
# mirror scene; spheres' reflections include self/other terms).
sq = lambda v: np.sqrt(np.maximum(v, 0))
e, R = SC4['e'], SC4['R']
bx, by, bz = SC4['Cb']; R2 = SC4['R2']
Aq = X**2 + Y**2 + e**2
onF = (e*0 + 1) > 0                      # floor pixels: primary t_f = 1 wins
discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
tA = (e*e - sq(discA))/Aq
bB = -bx*X - by*Y - (e - bz)*e
cB = bx**2 + by**2 + (e - bz)**2 - R2**2
discB = bB**2 - Aq*cB
tB = (-bB - sq(discB))/Aq
onA = (discA >= 0) & ((discB < 0) | (tA <= tB))
onB = (discB >= 0) & ~onA
floor_px = SC4['floor_on'] & ~onA & ~onB
# floor primary P = (x, y, 0), N = +z: reflected direction D' = (x, y, +e)
# = exactly the unfolded ray above. Sphere re-hit from P along D':
def rehit(C, Rr):
    aa = X**2 + Y**2 + e**2
    bb = (X-C[0])*X + (Y-C[1])*Y + (0-C[2])*e
    cc = (X-C[0])**2 + (Y-C[1])**2 + C[2]**2 - Rr**2
    dd = bb**2 - aa*cc
    tt = np.where(dd >= 0, (-bb - sq(dd))/aa, np.inf)
    return np.where(tt > 0, tt, np.inf)
t_refl = np.minimum(rehit((0., 0., 0.), R),
                    rehit((bx, by, bz), R2) if SC4['b_on'] else np.inf)

both = floor_px & np.isfinite(t_mirror) & np.isfinite(t_refl)
gap = np.abs(t_mirror - t_refl)[both]
print(f"floor pixels where both paths hit a sphere: {both.sum()}")
print(f"max |t_mirror - t_reflected| = {gap.max():.3e}   <- the flip is exact")
print("The mirrored scene and the reflected ray are one object, measured twice.")

vis = np.zeros(X.shape + (3,))
vis[floor_px & np.isfinite(t_refl)] = [0.55, 0.7, 0.95]
vis[floor_px & ~np.isfinite(t_refl)] = [0.12, 0.14, 0.2]
vis[~floor_px] = [0.05, 0.05, 0.06]
plt.figure(figsize=(4.5, 4)); plt.imshow(vis)
plt.title("what the floor mirror sees (reflected sphere footprints)")
plt.axis('off'); plt.show()
Expect: max |Δt| ~1e-15 between the mirrored-scene primary render and the composed bounce — the rasterizer's mirror trick, promoted to a theorem and verified numerically.

🐍 Cell 6 (NEW) — the ε-hack, measured: acne, and the perturbation the fudge buys

· · paste into colab.new → Run
# (a) ACNE CENSUS: with NO self-exclusion at all, reflected rays from sphere
#     pixels re-hit their spawning sphere at t' = 0 (exactly, in exact
#     arithmetic; at t' ~ 1e-16 in floats). Count the poisoned pixels:
sq = lambda v: np.sqrt(np.maximum(v, 0))
e, R = SC4['e'], SC4['R']
Aq = X**2 + Y**2 + e**2
discA = e*e*R*R - (e*e - R*R)*(X**2 + Y**2)
tA = (e*e - sq(discA))/Aq
onA = discA >= 0
Px, Py, Pz = tA*X, tA*Y, e*(1 - tA)
Nx, Ny, Nz = Px/R, Py/R, Pz/R
Dn = X*Nx + Y*Ny - e*Nz
Dpx, Dpy, Dpz = X - 2*Dn*Nx, Y - 2*Dn*Ny, -e - 2*Dn*Nz
aa = Dpx**2 + Dpy**2 + Dpz**2
bb = Px*Dpx + Py*Dpy + Pz*Dpz                      # blocker = spawner: c' = 0
dd = bb**2  # c' = 0 exactly: disc = b'^2
t_self = np.where(aa > 0, (-bb - np.sqrt(dd))/aa, np.inf)
self_hits = onA & np.isfinite(t_self) & (np.abs(t_self) < 1e-12)
print(f"sphere-A pixels whose reflected ray re-hits A at t' ≈ 0: "
      f"{self_hits.sum()} of {onA.sum()}  <- acne, without exclusion")
print("(the factorization t'(a't' + 2b') = 0 is exact; the near-root is 0)")

# (b) STRICT vs HACK: the epsilon offset perturbs every downstream quantity.
img_strict, _, _ = render_symbolic_rfl(X, Y)
img_hack_num, _, _ = render_numeric_rfl(256, 256, {**SC4, 'selftest': 'hack'})
d_hack = np.abs(img_strict - img_hack_num)
same = pcode_num == pcode_sym
interior = d_hack[same]
print(f"\nstrict-symbolic vs hack-numeric, interior pixels: "
      f"max |ΔI| = {interior.max():.3e}")
print(f"pixels differing by > 1e-6: {(interior > 1e-6).sum()} of {same.sum()}")
print("The hack shifts every secondary origin by 1e-4: hits, shading, shadows,")
print("all slightly wrong, everywhere the bounce matters. The strict test")
print("excludes exactly the degenerate root and perturbs nothing.")

fig, axs = plt.subplots(1, 2, figsize=(9, 4))
axs[0].imshow(np.where(self_hits[..., None], [1.0, 0.3, 0.2], img_sym))
axs[0].set_title("acne sites (red): the t′=0 root, excluded strictly")
axs[1].imshow(np.clip(d_hack*1e4, 0, 1))
axs[1].set_title("|strict − hack| ×10⁴: the ε-fudge's footprint")
for ax in axs: ax.axis('off')
plt.show()
Expect: hundreds of acne sites concentrated at sphere rims (where D′ grazes); the hack's error footprint faint but global — an approximate patch for an exactly-solvable problem.

🐍 Cell 7 (NEW) — pitfalls, live: raw Heaviside at the third level; the Moiré warning

· · paste into colab.new → Run
# (a) RAW Heaviside assembly, now with a composed bounce: NaNs from every
#     unselected chamber at BOTH levels. The census is per-level:
img_raw, _, _ = render_symbolic_rfl(X, Y, SC4, regularized=False)
nan_px = np.isnan(img_raw).any(axis=-1)
print(f"raw assembly NaN pixels: {nan_px.sum()} of {nan_px.size}")
img_reg, _, _ = render_symbolic_rfl(X, Y, SC4, regularized=True)
print(f"after regularization:    {np.isnan(img_reg).sum()}")
print("Three levels of chambers (primary x bounce x shadow masks): the 0·NaN")
print("trap waits at every multiply. Regularize every sqrt, or select, always.")

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(nan_px, cmap='autumn'); axs[0].set_title("raw NaN support")

# (b) MOIRE: the checker is chi = sign(sin(pi x) sin(pi y)) -- symbolic, and
#     NOT band-limited. Magnify the grazing band and the mirror checker aliases.
SC_mag = {**SC4, 'krf': 0.9}
Xz, Yz = grid(256, 256)
# zoomed window: emulate by shrinking the view half-extent
half_full = 1.25*SC4['e']*SC4['R']/np.sqrt(SC4['e']**2 - SC4['R']**2)
def grid_zoom(W, H, halfw, cy=0.0):
    xs = np.linspace(-halfw, halfw, W)
    ys = np.linspace(cy + halfw, cy - halfw, H)
    return np.meshgrid(xs, ys)
Xz2, Yz2 = grid_zoom(256, 256, half_full*0.04, cy=-0.55*half_full)
img_zoom, _, _ = render_symbolic_rfl(Xz2, Yz2, SC_mag)
axs[1].imshow(img_zoom); axs[1].set_title("grazing band, magnified: Moiré")
# full view for reference
axs[2].imshow(render_symbolic_rfl(Xz, Yz, SC_mag)[0])
axs[2].set_title("k_rf = 0.9 full view")
for ax in axs: ax.axis('off')
plt.show()
print("\nThe Moiré is the honest warning: a symbolic function evaluated at pixel")
print("CENTERS is a sampler, not an integrator. Chapter 6 must derive the box")
print("filter -- integrate the checker (and the chambers) over the pixel square --")
print("or accept aliasing as the price of point sampling, symbolic or not.")
Expect: thousands of raw NaNs, zero after regularization; and in the magnified grazing band, checkerboard Moiré — aliasing of an exact function, the case for Chapter 6's derived filter.

🐍 Cell 8 (NEW) — the frontier, updated: bounce 2 compounds; glossy reflection is an integral (teaser for Ch. 5–6)

· · paste into colab.new → Run
# (a) THE SECOND BOUNCE, PRICED. Reflect TWICE: D'' = D' - 2(D'.N')N', with the
#     second hit's normal N' itself a composed expression. Count, don't trust:
Dp_ops = sp.count_ops(Dp)
print(f"bounce-1 direction D':        {Dp_ops} ops/component (measured, Cell 2)")
# A second bounce substitutes the FIRST bounce's hit point (itself containing
# sqrt(disc_r)) into the reflection formula again:
print("bounce-2 requires N' at the reflected hit: P' contains sqrt(disc_r),")
print("so N'' = (P' - C)/R and D'' = D' - 2(D'.N')N' compose sqrt-of-composed-")
print("expressions. Estimated op growth per bounce: ~3-4x (measured at bounce 1")
print("in Cell 2: raw ~%.0f ops; bounce 2 lands in the thousands even after cse)."
      % 400)
print("Exact, in principle. Usable, no. The closed-form arc ends by attrition.")

# (b) GLOSSY REFLECTION IS AN INTEGRAL -- the real frontier:
print("""
Perfect mirror:        I = I_reflected (one composed expression)      [this chapter]
Glossy (blurred) reflection:
    I(P) = integral over the lobe of  I_reflected(D') * W(D'.mirror) dD'
The integrand is a composed chamber expression; the domain is a solid angle.
No closed form exists even for one sphere: the lobe crosses visibility
chambers, and the integral of a Piecewise over its own boundary curves is
a semi-algebraic quadrature problem. Options, honestly labeled:
  * Monte Carlo over the lobe (numerical -- the symbolic form still helps:
    the integrand is exact and differentiable),
  * quadrature with chamber-aware splitting (semi-symbolic),
  * give up and call it roughness (empirical -- the industry's answer).
Chapter 5 (textures as symbolic functions) and Chapter 6 (the derived box
filter) pick up the two remaining tractable threads of the classic arc.""")
Expect: the second bounce priced as compounding swell, and the glossy frontier named precisely — an integral over a chambered integrand, where the symbolic object remains useful even when the closed form dies.

§6 · The honest ledger — Chapter-4 deltas

DimensionChapter 3Chapter 4
Objectstwo spheres+ the plane z=0 (no silhouette, hit at t=1 — the Chapter-1 screen, in the scene)
Reflection—one bounce; symbolic side is substitution/composition of expressions; mirror flip = change of variables, verified (Cell 5)
Self-intersect—numeric ε-hack vs. exact strict t′ > 0 (roots factor: 0 and −2b′/a′); both toggleable, the hack's footprint measured (Cell 6)
ShadowsShadow Lemma, primary feelerssame lemma on reflected feelers — composition-invariant (§3)
Chambers9 primary regionsprimary × reflected chamber codes; flip census unchanged in method (Cell 3)
Totality taxboth spheres, every pixel+ the bounce path for every surface pixel; ~3× and honestly tabled (Cell 4)
Floor shading—symbolic checker χ = sign(sin πx · sin πy); Moiré under magnification = Chapter 6's warning (Cell 7b)
Failure modescross-chamber NaNsthe same, at three chamber levels (Cell 7a)
Frontier∀ over blockers = productbounce 2 = compounding swell; glossy = an integral, not an expression (Cell 8)

§7 · Pitfalls gallery, continued

  1. The t′ = 0 root is exact. The reflected ray's self-intersection is not a rounding accident; the quadratic factors. Strict t′ > 0 excludes exactly it; the ε-hack perturbs everything else to avoid it. Prefer the theorem to the fudge — and when you must fudge, measure the footprint (Cell 6).
  2. Acne clusters at silhouettes. Where the view grazes a sphere, N·D → 0 and the reflected direction barely leaves the surface — the acne sites concentrate at rims (visible in Cell 6's census map). Boundary pixels are where all of this book's approximations go to be tested.
  3. Cross-chamber NaNs, level 3. Raw Heaviside with a composed bounce: every unselected chamber at both levels is a NaN source. Regularize or select; there is no third option.
  4. The checker aliases. A symbolic function is exact; sampling it is not. Magnify the grazing band and the mirror floor Moirés. The box filter must be derived — Chapter 6.
  5. Ordering ties, again. Floor at t=1 vs. sphere grazing the plane; reflected hit grazing the floor. The ≤ convention must match across pipelines or the flip census lights up the tie curve.

§8 · A brief history, continued

Turner Whitted, 1980 — "An Improved Illumination Model for Shaded Display" — put reflection, refraction, and shadow feelers into one recursive ray tracer and gave rendering its canonical image: spheres on a checkered plane, bouncing light at each other. (This chapter is, deliberately, Whitted's scene with two lights and honest weights.) The ε-offset hack is likewise as old as the recursive bounce itself — every ray tracer carries one; few can say exactly what it costs. The mirror-flip trick is older than ray tracing: planar reflection by rendering the mirrored scene is the standard rasterizer idiom (the stencil-reflection era of the late 1990s ran on it), and it survives here as a change of variables in a composed expression. In the other world: composition has always been the CAS's double-edged primitive — Macsyma users of the 1970s watched substituting an expression into itself produce exactly this chapter's swell, and common subexpression elimination (Cocke 1970; Macsyma's optimize()) was invented to claw it back. The plane itself is the quiet joke of the chapter: the screen these chapters have been rendering onto since Chapter 1 was a mirror all along.

§9 · Roadmap

❦
Set by the Press · single HTML5 file · the expression is the ground truth.
Section ids are stable (s1…s9, sim, lab, cell1…cell8) — cite the id when requesting revisions.