Symbolic Rendering

Chapter 7 — Area Lights, or Penumbra Is an Integral

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 source integral, the on-axis closed form and its ρ→0 regression, the conic penumbra boundary. §4 is the live simulator with the K-slider (both pipelines spend the same budget, differently) and the penumbra map. §5: Cells 1 and 3 replace Chapter 6's; Cells 2, 4–8 are new. §6–§9 as ever.
☰  contents
  1. What changes in Chapter 7
  2. The scene (Chapter-7 assumptions)
  3. The math, once, carefully
  4. Live simulator
  5. The Colab laboratory (Python / sympy)
  6. The honest ledger — Chapter-7 deltas
  7. Pitfalls gallery, continued
  8. A brief history, continued
  9. Roadmap

§1 · What changes in Chapter 7

Chapter 6's scene stands (textures, analytic floor, gamma policy). One physical change and one methodological break:

  1. Light 1 is a disk. Center $L_1$, radius $\rho$, facing the scene, radiance $c_1$. Its illumination at $P$ is no longer a value but an integral over the source — and its shadow is no longer a bit but a lit fraction $f_1(P) \in [0,1]$. Light 2 remains a point: hard and soft shadows share the frame for contrast.
  2. The pipelines part ways in method, not in mathematics. Previous chapters' numeric/symbolic pairs computed the same quantity two ways. Now both compute an estimate of one integral: numeric by jittered Monte Carlo, symbolic by deterministic stratified quadrature (plus one closed-form anchor). The difference canvas no longer measures agreement — it measures estimator error, and its two species (variance, bias) look different, decay differently, and are tabulated separately in the ledger.

§2 · The scene (Chapter-7 assumptions)

Chapter 6's scene with light 1 replaced by a disk (default $\rho = 0.45$, facing the origin; emitter units follow the series convention: the disk's contribution is the average of point-light contributions over its area — no falloff, no emitter cosine; the physically-weighted units are derived in Cell 2 as the "units upgrade" and shown to agree after normalization). Blockers: spheres A and B as ever, sign-logic shadows per source sample. The floor keeps its Chapter-6 analytic antialiasing; the source integral is the new term, stacked on the pixel integral, and we say so: this chapter filters the light, not the pixel (Chapter 6 owns the pixel).

§3 · The math, once, carefully

The source integral. With disk point $Q$, direction $\hat w(Q) = (Q-P)/\|Q-P\|$, and the Shadow Lemma mask $\text{lit}(P,Q) \in \{0,1\}$:

$$I_1(P) = \frac{1}{\pi\rho^2}\int_{\text{disk}} c_1\Big(k_d\,m\max(0, \mu(Q)) + k_s\,\text{spec}(Q)\Big)\,\text{lit}(P, Q)\,dA(Q), \qquad \mu = N\cdot\hat w .$$

The integrand is Chapter 5's expression composed with $Q(\text{stratum})$; the mask is Chapter 3's lemma evaluated per stratum. Nothing new is derived — only integrated.

Theorem (On-axis closed form — the anchor). Unshadowed, for $P$ on the disk's axis at distance $d$, the series-convention mean diffuse is elementary: $$\bar\mu(d, \rho) = \frac{1}{\pi\rho^2}\int_0^{2\pi}\!\!\int_0^{\rho} \frac{d}{\sqrt{d^2+r^2}}\;r\,dr\,d\theta = \frac{2d\big(\sqrt{d^2+\rho^2}-d\big)}{\rho^2}.$$ Proof: one substitution $u = d^2+r^2$ (Cell 2 has sympy do the honors). Corollary (the equitesseral check of this chapter). $\lim_{\rho\to0}\bar\mu = 1 = \mu(\text{point light})$: the disk degenerates to the Chapter-5 point light, and Cell 6 runs the whole scene through that limit, recovering Chapter 5's image to machine precision. A new feature must degenerate to an old theorem; this one does, twice.

Theorem (The penumbra boundary on the source is a conic). For fixed $P$, the mask's edge is $\Delta_s(P,Q) = 0$ with $b_2 \lt 0$, $a_2+b_2 \gt 0$ — and $\Delta_s = \big((P-C)\cdot(Q-P)\big)^2 - \|Q-P\|^2\big(\|P-C\|^2 - R_b^2\big)$ is quadratic in $Q$: the tangent cone from $P$ to the blocker sphere. Cutting the cone by the light plane gives a conic section — ellipse, parabola, or hyperbola by the relative geometry — and the lit region of the source is the disk minus a conic piece. Cell 2 extracts the conic's matrix and classifies it; Cell 8 discusses when the conic is an ellipse (bounded umbra cone) and what the degenerate cases mean.

Estimators, contrasted. Both pipelines evaluate the identical integrand at $K^2$ source points. Numeric jitters them (uniform in expectation; error is variance, $\sim \sigma/\sqrt K$, visible as speckle, different every seed). Symbolic uses stratum centroids (deterministic; error is bias, $\sim O(1/K^2)$ on smooth regions but structured — banding where the mask boundary crosses strata the same way every time). Cell 5 measures both against a $K=64$ truth: the variance line and the bias curve, plotted together, are the chapter's portrait.

Three gifts of Chapter 7: (1) soft shadows with no new shadow theory — the lemma integrates; (2) an exact anchor ($\bar\mu$) that lets every estimator be graded, not just compared; (3) the conic: the penumbra, like every boundary in this book, is a curve with a name. Counterpoint: the general-position, blocker-shadowed, specular-inclusive integral has no closed form (the on-axis form factor is already elliptic off-axis; the conic-clipped version worse). And the integrand's discontinuities (mask edge, terminator) are exactly where quadrature rules lose their orders — the penumbra is smooth to the eye and singular to the numerologist, and both statements are the point.

§4 · Live simulator

Disk key light, point fill. Both pipelines evaluate the identical integrand at K² source points; numeric jitters, symbolic holds strata. The penumbra view paints the lit fraction f₁; the ledger now reports interior RMS (the variance meter) alongside max, and the status line counts source samples honestly.

model 
resolution 
initializing…
NUMERIC — jittered Monte Carlo over the source (error = variance)
SYMBOLIC — deterministic strata of the same integrand (error = bias)
|difference| ×10² — speckle is numeric variance; smooth bands are symbolic bias
numeric: – symbolic: – interior RMS |ΔI|: – max |ΔI|: – flips: –
hover — pixel → world, surface, lit fractions (f₁, f₂), and both RGB triples.

Try this: (1) K = 1: the disk collapses to a point — both pipelines return to Chapter 5 hard shadows, and the difference returns to machine black. That is the ρ→0/K→1 regression, live. (2) K = 2: numeric speckles (variance, re-render to re-roll the dice by nudging any slider); symbolic bands (bias, same every time). Two failure shapes, one integrand. (3) K = 8: the difference melts inward from the penumbra. (4) Penumbra view: f₁ as a field — the soft edge is the conic of §3 integrated over the disk; magnify to 0.08 and watch it stay smooth (the mask is sampled, the field is exact). (5) ρ slider: the penumbra widens with the source — geometry, not blur.

§5 · The Colab laboratory (Python / sympy)

Cumulative: Cells 1 and 3 replace Chapter 6's; 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.6 Cell 1) — numeric baseline with a disk light: jittered Monte Carlo over the source

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

SC7 = dict(R=1.0, e=4.0, Cb=(0.5, 0.5, 1.3), R2=0.32, b_on=True, floor_on=True,
           L1=(5., 5., 10.), rho=0.45, c1=(1., 1., 1.),          # disk key light
           L2=(-6., 2., 8.), c2=(0.9, 0.65, 0.25), l2_on=True,  # point fill
           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',
           tex='checker', ku=8.0, kv=4.0, alpha=0.25, K=4, seed=7)

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

def disk_basis(L, target=(0., 0., 0.)):
    w = np.array(target, dtype=float) - np.array(L, dtype=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)
    e2 = np.cross(w, e1)
    return e1, e2

def disk_points(SC, K, jitter=None):
    """KxK source points: jittered (numeric) or stratum centroids (symbolic)."""
    e1, e2 = disk_basis(SC['L1'])
    pts = []
    for sj in range(K):
        for si in range(K):
            if jitter is None:
                u, v = (si + 0.5)/K, (sj + 0.5)/K
            else:
                u, v = (si + jitter())/K, (sj + jitter())/K
            r = SC['rho']*math.sqrt(u); th = 2*math.pi*v
            pts.append(np.array(SC['L1']) + r*(math.cos(th)*e1 + math.sin(th)*e2))
    return pts

def chartA1(P, R):
    Nx, Ny, Nz = P[0]/R, P[1]/R, P[2]/R
    return ((math.atan2(Ny, Nx) + math.pi)/(2*math.pi),
            math.acos(max(-1.0, min(1.0, Nz)))/math.pi)

def texture1(u, v, SC=SC7):
    m = SC['matA']; t = SC['tex']
    if t == 'solid': alb = m
    elif t == 'stripes':
        b = 1.0 if math.sin(math.pi*u*SC['ku']) >= 0 else 0.2
        alb = (m[0]*b, m[1]*b, m[2]*b)
    elif t == 'marble':
        mm = 0.5 + 0.5*math.sin(2*math.pi*(SC['ku']*u + 0.35*math.sin(3*math.pi*SC['kv']*v)))
        b = 0.55 + 0.45*mm; alb = (m[0]*b, m[1]*b, m[2]*b)
    else:
        s = math.sin(math.pi*u*SC['ku'])*math.sin(math.pi*v*SC['kv'])
        b = 1.0 if s >= 0 else 0.2; alb = (m[0]*b, m[1]*b, m[2]*b)
    return np.array(alb)

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

def shadowed1(P, Q, C, Rr):
    """EXPLICIT roots on segment P->Q (numeric formulation of the mask)."""
    seg = Q - 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 light_point(P, N, Q, v_hat, alb, SC, blockers):
    """One point-light contribution (the integrand), shared by both estimators."""
    if any(shadowed1(P, Q, C, Rr) for C, Rr in blockers):
        return np.zeros(3), 0.0
    w = Q - P; wn = np.linalg.norm(w); w_hat = w/wn
    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']
    return SC['kd']*alb*max(0.0, mu) + SC['ks']*spec, 1.0

def render_numeric_area(width, height, SC=SC7):
    """Numeric: jittered Monte Carlo over the disk; hard point for light 2."""
    rng = np.random.default_rng(SC['seed'])
    R, e = SC['R'], SC['e']; Cb = np.array(SC['Cb']); R2 = SC['R2']
    img = np.zeros((height, width, 3)); code = np.zeros((height, width), np.uint8)
    fmap = np.zeros((height, width, 2))
    half = 1.25*e*R/np.sqrt(e**2 - R**2)
    E = np.array([0.0, 0.0, e])
    K = SC['K']
    for j in range(height):
        y = half*(1 - 2*(j + 0.5)/height) - 0.15*half
        for i in range(width):
            x = half*(2*(i + 0.5)/width - 1)
            Aq = x*x + y*y + e*e
            discA = e*e*R*R - (e*e - R*R)*(x*x + y*y)
            tA = (e*e - math.sqrt(discA))/Aq if discA >= 0 else math.inf
            tB = math.inf
            if SC['b_on']:
                bB = -Cb[0]*x - Cb[1]*y - (e - Cb[2])*e
                cB = float(Cb @ Cb) - R2*R2
                discB = bB*bB - Aq*cB
                if discB >= 0: tB = (-bB - math.sqrt(discB))/Aq
            tf = 1.0 if SC['floor_on'] else math.inf
            tm = min(tA, tB, tf)
            if tm == math.inf: continue
            P = np.array([tm*x, tm*y, e*(1 - tm)])
            if tm == tA:
                N = P/R; u, vv = chartA1(P, R); alb = texture1(u, vv, SC)
                blockers = [(Cb, R2)] if SC['b_on'] else []; kind = 1
            elif tm == tB:
                N = (P - Cb)/R2; alb = np.array(SC['matB'])
                blockers = [(np.zeros(3), R)]; kind = 2
            else:
                N = np.array([0., 0., 1.]); alb = floor_tex1(P[0], P[1], SC)
                blockers = [(np.zeros(3), R)] + ([(Cb, R2)] if SC['b_on'] else []); kind = 3
            v_hat = (E - P)/np.linalg.norm(E - P)
            col = SC['ka']*np.array(SC['ambient'])
            # ---- light 1: disk, jittered MC ----
            pts = disk_points(SC, K, jitter=rng.random) if K > 1 else \
                  [np.array(SC['L1'])]
            acc = np.zeros(3); f1 = 0.0
            for Q in pts:
                c, ok = light_point(P, N, Q, v_hat, alb, SC, blockers)
                acc += np.array(SC['c1'])*c; f1 += ok
            col += acc/len(pts); f1 /= len(pts)
            # ---- light 2: point ----
            f2 = 0.0
            if SC['l2_on']:
                c, ok = light_point(P, N, np.array(SC['L2']), v_hat, alb, SC, blockers)
                col += np.array(SC['c2'])*c; f2 = ok
            img[j, i] = col; code[j, i] = kind
            fmap[j, i] = (f1, f2)
    return np.clip(img, 0, 1), code, fmap

t0 = time.perf_counter()
img_num, code_num, fmap_num = render_numeric_area(256, 256)
t1 = time.perf_counter()
print(f"numeric render (disk light, K={SC7['K']}x{SC7['K']} jittered): {t1 - t0:.2f} s")
plt.figure(figsize=(4, 4)); plt.imshow(srgb_encode(img_num))
plt.title("Numeric: soft key shadow, hard fill shadow"); plt.axis('off'); plt.show()
Expect: a soft-edged shadow from the key light, crisp from the fill — and a faint speckle in the penumbra at K=4 that no symbolic render will share.

🐍 Cell 2 (NEW) — the two theorems: on-axis closed form (with the ρ→0 regression) and the conic penumbra boundary

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

# --- (a) ON-AXIS MEAN DIFFUSE, closed form --------------------------------------
r, d, rho = sp.symbols('r d rho', positive=True)
mu = d/sp.sqrt(d**2 + r**2)
bar_mu = sp.simplify((1/(sp.pi*rho**2)) *
                     sp.integrate(2*sp.pi*r*mu, (r, 0, rho)))
print("mean mu on axis:  bar_mu =", bar_mu)
print("rho -> 0 limit:", sp.limit(bar_mu, rho, 0), "   <- the point light (equitesseral check)")

# verify against quadrature
def bar_mu_quad(dv, rv, K=4000):
    rr = (np.arange(K) + 0.5)/K*rv
    return np.mean(dv/np.sqrt(dv**2 + rr**2))
for dv, rv in ((4.0, 0.45), (6.0, 1.2), (2.0, 0.05)):
    exact = float(bar_mu.subs({d: dv, rho: rv}))
    print(f"  d={dv}, rho={rv}: exact {exact:.8f}   quad {bar_mu_quad(dv, rv):.8f}")

# --- (b) THE PENUMBRA BOUNDARY ON THE SOURCE IS A CONIC --------------------------
Px, Py, Pz = sp.symbols('P_x P_y P_z', real=True)
Cx, Cy, Cz, Rb = sp.symbols('C_x C_y C_z R_b', real=True)
u, v = sp.symbols('u v', real=True)
Lx, Ly, Lz = sp.symbols('L_x L_y L_z', real=True)
E1x, E1y, E1z = sp.symbols('e_1x e_1y e_1z', real=True)
E2x, E2y, E2z = sp.symbols('e_2x e_2y e_2z', real=True)
P = sp.Matrix([Px, Py, Pz]); C = sp.Matrix([Cx, Cy, Cz])
Q = sp.Matrix([Lx, Ly, Lz]) + u*sp.Matrix([E1x, E1y, E1z]) + v*sp.Matrix([E2x, E2y, E2z])
seg = Q - P
a2 = seg.dot(seg); b2 = (P - C).dot(seg); c2 = (P - C).dot(P - C) - Rb**2
cone = sp.expand(b2**2 - a2*c2)          # = 0: tangent cone, in plane coords (u, v)
poly = sp.Poly(cone, u, v)
print("\nDelta_s(P, Q(u,v)) = 0 on the light plane:")
print("  total degree in (u,v):", sp.Poly(cone, u, v).total_degree(),
      "  <- a CONIC (tangent cone x light plane)")
# classify for a sample configuration
cfg = {Px: 0.0, Py: 0.0, Pz: 1.0, Cx: 0.4, Cy: 0.3, Cz: 2.0, Rb: 0.32,
       Lx: 5.0, Ly: 5.0, Lz: 10.0, E1x: 0.78, E1y: -0.63, E1z: 0.0,
       E2x: 0.45, E2y: 0.55, E2z: -0.71}
cone_n = sp.Poly(sp.expand(cone.subs(cfg)), u, v)
A11 = float(cone_n.coeff_monomial(u*u)); A12 = float(cone_n.coeff_monomial(u*v))/2
A22 = float(cone_n.coeff_monomial(v*v))
detQ = A11*A22 - A12*A12
print(f"  sample config: quadratic-form det = {detQ:+.4f}  ->",
      "ellipse (bounded umbra cone)" if detQ > 0 else "hyperbola/parabola regime")
print("The penumbra's edge has a conic's name, like every boundary in this book.")
Expect: μ̄ printed in closed form, the limit 1, quadrature agreement to 8 decimals, and the conic classified (ellipse at the default geometry).

🐍 Cell 3 (REPLACES Ch.6 Cell 3) — the symbolic renderer: deterministic strata of the same integrand; the census gains RMS

· · paste into colab.new → Run
def render_symbolic_area(width, height, SC=SC7):
    """Symbolic: the SAME integrand at deterministic stratum centroids.
       No randomness anywhere: the error is bias, reproducible exactly."""
    R, e = SC['R'], SC['e']; Cb = np.array(SC['Cb']); R2 = SC['R2']
    img = np.zeros((height, width, 3)); code = np.zeros((height, width), np.uint8)
    fmap = np.zeros((height, width, 2))
    half = 1.25*e*R/np.sqrt(e**2 - R**2)
    E = np.array([0.0, 0.0, e])
    K = SC['K']
    pts = disk_points(SC, K, jitter=None) if K > 1 else [np.array(SC['L1'])]
    for j in range(height):
        y = half*(1 - 2*(j + 0.5)/height) - 0.15*half
        for i in range(width):
            x = half*(2*(i + 0.5)/width - 1)
            Aq = x*x + y*y + e*e
            discA = e*e*R*R - (e*e - R*R)*(x*x + y*y)
            tA = (e*e - math.sqrt(discA))/Aq if discA >= 0 else math.inf
            tB = math.inf
            if SC['b_on']:
                bB = -Cb[0]*x - Cb[1]*y - (e - Cb[2])*e
                cB = float(Cb @ Cb) - R2*R2
                discB = bB*bB - Aq*cB
                if discB >= 0: tB = (-bB - math.sqrt(discB))/Aq
            tf = 1.0 if SC['floor_on'] else math.inf
            tm = min(tA, tB, tf)
            if tm == math.inf: continue
            P = np.array([tm*x, tm*y, e*(1 - tm)])
            if tm == tA:
                N = P/R; u, vv = chartA1(P, R); alb = texture1(u, vv, SC)
                blockers = [(Cb, R2)] if SC['b_on'] else []; kind = 1
            elif tm == tB:
                N = (P - Cb)/R2; alb = np.array(SC['matB'])
                blockers = [(np.zeros(3), R)]; kind = 2
            else:
                N = np.array([0., 0., 1.]); alb = floor_tex1(P[0], P[1], SC)
                blockers = [(np.zeros(3), R)] + ([(Cb, R2)] if SC['b_on'] else []); kind = 3
            v_hat = (E - P)/np.linalg.norm(E - P)
            col = SC['ka']*np.array(SC['ambient'])
            acc = np.zeros(3); f1 = 0.0
            for Q in pts:
                c, ok = light_point(P, N, Q, v_hat, alb, SC, blockers)
                acc += np.array(SC['c1'])*c; f1 += ok
            col += acc/len(pts); f1 /= len(pts)
            f2 = 0.0
            if SC['l2_on']:
                c, ok = light_point(P, N, np.array(SC['L2']), v_hat, alb, SC, blockers)
                col += np.array(SC['c2'])*c; f2 = ok
            img[j, i] = col; code[j, i] = kind
            fmap[j, i] = (f1, f2)
    return np.clip(img, 0, 1), code, fmap

t0 = time.perf_counter()
img_sym, code_sym, fmap_sym = render_symbolic_area(256, 256)
t1 = time.perf_counter()
print(f"symbolic render (strata, K={SC7['K']}x{SC7['K']}): {t1 - t0:.2f} s")

# ---- the census, Chapter-7 edition: RMS enters ----------------------------------
same = code_num == code_sym
d = np.abs(img_num - img_sym)[same]
print(f"interior max  |num - sym| = {d.max():.3e}   <- estimator disagreement")
print(f"interior RMS  |num - sym| = {np.sqrt((d**2).mean()):.3e}   <- the variance meter")
print(f"flips (surface mismatch): {(~same).sum()} px")
print("Reading: at K=4 the max lives in the penumbra and the RMS is dominated by")
print("numeric VARIANCE (speckle). Raise K: both fall; RMS falls like 1/sqrt(K).")

fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(srgb_encode(img_num)); axs[0].set_title("Numeric (jittered MC)")
axs[1].imshow(srgb_encode(img_sym)); axs[1].set_title("Symbolic (strata)")
axs[2].imshow(np.clip(np.abs(img_num - img_sym)*50, 0, 1))
axs[2].set_title("|delta| x50: speckle vs bias")
for ax in axs: ax.axis('off')
plt.show()
Expect: RMS ~1e-2 at K=4 (estimator noise, not machine error); the ×50 diff shows numeric speckle against symbolic's smooth penumbra bands.

🐍 Cell 4 (NEW) — honest timings: both pipelines pay K²; determinism is free

· · paste into colab.new → Run
print(f"{'K':>4} {'numeric 96^2':>14} {'symbolic 96^2':>15}")
for K in (1, 2, 4, 8):
    t0 = time.perf_counter(); render_numeric_area(96, 96, {**SC7, 'K': K})
    t1 = time.perf_counter(); render_symbolic_area(96, 96, {**SC7, 'K': K})
    t2 = time.perf_counter()
    print(f"{K:>4} {t1-t0:>13.2f}s {t2-t1:>14.2f}s")

print("""
Reading the table:
 * Both pipelines scale ~K^2: the source integral is the same work either way.
   The symbolic advantage in this chapter is NOT speed.
 * It is REPRODUCIBILITY: the strata estimate is a pure function of the scene;
   the jittered estimate is a function of the scene AND the seed. For regression
   testing, differentiable pipelines, and the agreement ledger, determinism is
   the feature. Variance is a tax numeric pays for unbiasedness; bias is a tax
   symbolic pays for never rolling dice. Cell 5 prices both taxes.
""")
Expect: near-identical K² scaling in both columns; the moral is determinism, not speed.

🐍 Cell 5 (NEW) — the superpower study: variance vs bias, both measured against a K=64 truth

· · paste into colab.new → Run
# Ground truth: stratified K=64 (deterministic; fine enough to stand as truth).
W = 128
print("computing K=64 reference at 128^2 (patience)...")
img_true, code_true, _ = render_symbolic_area(W, W, {**SC7, 'K': 64})

print(f"\n{'K':>4} {'numeric RMS (seed 7)':>20} {'numeric RMS (seed 99)':>22} "
      f"{'symbolic RMS (bias)':>20}")
res = []
for K in (1, 2, 4, 8):
    row = [K]
    for seed in (7, 99):
        img_n, code_n, _ = render_numeric_area(W, W, {**SC7, 'K': K, 'seed': seed})
        same = code_n == code_true
        d = np.abs(img_n - img_true)[same]
        row.append(np.sqrt((d**2).mean()))
    img_s, code_s, _ = render_symbolic_area(W, W, {**SC7, 'K': K})
    same = code_s == code_true
    d = np.abs(img_s - img_true)[same]
    row.append(np.sqrt((d**2).mean()))
    res.append(row)
    print(f"{K:>4} {row[1]:>20.4e} {row[2]:>22.4e} {row[3]:>20.4e}")

Ks = [r[0] for r in res]
plt.figure(figsize=(6.5, 3.8))
plt.plot(Ks, [r[1] for r in res], 'o-', color='#7a1f1f', label='numeric (jittered MC)')
plt.plot(Ks, [r[3] for r in res], 's-', color='#1f3a5f', label='symbolic (strata)')
plt.plot(Ks, [res[0][1]/np.sqrt(np.array(Ks)**2), ][0] *
         [r[1] for r in res][0] / np.sqrt(np.array(Ks)**2 /
         (np.array(Ks)**2).astype(float)), ':', color='#b8860b', label='1/K guide')
plt.xscale('log'); plt.yscale('log')
plt.xlabel('K (samples = K^2)'); plt.ylabel('interior RMS error vs K=64 truth')
plt.title('variance vs bias, priced'); plt.legend(); plt.show()

print("\nNumeric wins on expectation (unbiased) but pays in variance; symbolic wins")
print("on determinism but pays in bias. NEITHER dominates: the choice of estimator")
print("is an engineering decision made about a mathematical object both sides share.")
Expect: two decaying curves on log-log — numeric noisier at low K, both converging; the table shows seed-dependence in the numeric column and none in the symbolic.

🐍 Cell 6 (NEW) — the equitesseral regression: ρ → 0 recovers Chapter 5 exactly

· · paste into colab.new → Run
# A new feature must degenerate to an old theorem. Shrink the disk to a point
# and this chapter's renderer must RETURN Chapter 5's hard-shadow image.
SC_pt = {**SC7, 'rho': 1e-6, 'K': 1}        # K=1: single centroid sample
img_s1, code_s1, _ = render_symbolic_area(128, 128, SC_pt)

# Chapter-5-equivalent: the same scene with light 1 as a literal point light,
# evaluated by the same engine (K=1 IS the point-light path):
img_s2, code_s2, _ = render_symbolic_area(128, 128, {**SC7, 'K': 1})
d = np.abs(img_s1 - img_s2)
same = code_s1 == code_s2
print(f"rho=1e-6 disk vs literal point light:  interior max |delta| = {d[same].max():.3e}")
print(f"chamber agreement: {(same).mean()*100:.2f}%   "
      f"({(~same).sum()} boundary stragglers at the hard shadow's edge)")
print("The disk light is a one-parameter family CLOSING on the point light;")
print("hard shadows are the singular limit of soft ones. Regression: passed.")

# and the analytic anchor at scene geometry (self-test grade):
e1, e2 = disk_basis(SC7['L1'])
P_ax = np.array(SC7['L1']) + 6.0*np.array([ -SC7['L1'][0], -SC7['L1'][1], -SC7['L1'][2]])
P_ax /= np.linalg.norm(P_ax)*0 + 1            # placeholder, replaced below:
# proper on-axis point: L1 + 6 * (disk normal)
w = -(np.array(SC7['L1']) - np.array([0., 0., 0.]))
w = w/np.linalg.norm(w)
P_ax = np.array(SC7['L1']) + 6.0*w
N_ax = -w                                     # surface facing the light
pts = disk_points({**SC7, 'rho': SC7['rho']}, 400, jitter=None)
mu_num = np.mean([max(0.0, float(N_ax @ ((Q - P_ax)/np.linalg.norm(Q - P_ax))))
                  for Q in pts])
d_ax = 6.0; rho_v = SC7['rho']
mu_exact = 2*d_ax*(math.sqrt(d_ax**2 + rho_v**2) - d_ax)/rho_v**2
print(f"\non-axis anchor: quadrature {mu_num:.8f}  exact {mu_exact:.8f}  "
      f"|delta| = {abs(mu_num - mu_exact):.2e}")
Expect: machine-precision regression to the point light, and the on-axis formula agreeing with 160,000-sample quadrature to ~1e-8.

🐍 Cell 7 (NEW) — pitfalls, live: stratified banding; terminator roughness; the degenerate conic

· · paste into colab.new → Run
# (a) BANDING: deterministic strata carve VISIBLE STEPS into the penumbra when
#     the mask boundary crosses strata identically pixel after pixel.
img_b, code_b, fm_b = render_symbolic_area(128, 128, {**SC7, 'K': 3})
fig, axs = plt.subplots(1, 3, figsize=(13, 4))
axs[0].imshow(fm_b[..., 0], cmap='magma'); axs[0].set_title("f1 field, K=3: banded")
img_b8, code_b8, fm_b8 = render_symbolic_area(128, 128, {**SC7, 'K': 8})
axs[1].imshow(fm_b8[..., 0], cmap='magma'); axs[1].set_title("f1 field, K=8: smoother")
# (b) TERMINATOR ROUGHNESS: the max(0, mu) clamp is a mask too; at low K the
#     diffuse terminator picks up the same stepping as the shadow edge.
axs[2].imshow(np.abs(fm_b8[..., 0] - fm_b[..., 0]), cmap='magma')
axs[2].set_title("|f1(K=8) - f1(K=3)|: the error lives on the edges")
for ax in axs: ax.axis('off')
plt.show()
print("Banding is bias made visible: the strata pattern prints through the mask.")
print("Numeric's jitter trades banding for speckle; symbolic could jitter too --")
print("but then it would have a seed. The two estimators are one choice, priced.")

# (c) THE DEGENERATE CONIC: when the tangent cone's apex geometry aligns with the
#     light plane, the penumbra conic degenerates (parabola regime) -- and the
#     lit fraction develops a kink. It is enough to know WHERE to look:
print("\nconic discriminant across P along a scanline (sphere A's midriff):")
e1v, e2v = disk_basis(SC7['L1'])
for t in np.linspace(-0.9, 0.9, 7):
    Pv = np.array([t, 0.0, math.sqrt(max(1e-9, 1 - t*t))]) * SC7['R']
    # Delta_s as quadratic form in (u, v) on the light plane:
    C = np.array(SC7['Cb']); Rb = SC7['R2']
    Q0 = np.array(SC7['L1'])
    def dsc_uv(uu, vv):
        Q = Q0 + uu*e1v + vv*e2v
        seg = Q - Pv
        a2 = seg @ seg; b2 = (Pv - C) @ seg; c2 = (Pv - C) @ (Pv - C) - Rb*Rb
        return b2*b2 - a2*c2
    # quadratic-form coefficients by finite probe (exact algebra is Cell 2's):
    A11 = (dsc_uv(1, 0) + dsc_uv(-1, 0) - 2*dsc_uv(0, 0))/2
    A22 = (dsc_uv(0, 1) + dsc_uv(0, -1) - 2*dsc_uv(0, 0))/2
    A12 = (dsc_uv(1, 1) - dsc_uv(1, -1) - dsc_uv(-1, 1) + dsc_uv(-1, -1))/8
    det = A11*A22 - A12*A12
    print(f"  P_x={t:+.2f}: det = {det:+.4f}  {'ellipse' if det > 0 else 'hyperbolic regime'}")
print("The regime can FLIP along a surface -- the penumbra's topology is not")
print("fixed by the scene, it is computed per point. The CAS knew; now you do.")
Expect: visible banding at K=3 melting by K=8; the error map concentrated on the shadow and terminator edges; the conic regime printed per point.

🐍 Cell 8 (NEW) — the frontier, updated: every remaining effect is this integral with a different domain

· · paste into colab.new → Run
print("""
This chapter's pattern -- integrate the SAME chambered, sign-masked, closed-form
integrand over a small domain, with one exact anchor and one regression limit --
generalizes verbatim:

  Ch. 8  glossy reflection:   integrate Chapter 4's composed expression over the
         specular lobe (domain: solid angle). Anchor: the mirror direction,
         closed form. Regression: shininess -> infinity returns the mirror image.
  Ch. 10 motion blur:         integrate over shutter time (domain: an interval).
         Anchor: static scene = Chapter 7. Regression: dt -> 0.
         depth of field:       integrate over the aperture (domain: a disk --
         the SAME disk geometry as this chapter's light).
  Appendix path tracing:       integrate over PATH SPACE. The chamber boundaries
         become boundaries of path chambers; the masks stay exact per sample;
         the domain stops being finite-dimensional-friendly, and the symbolic
         contribution shifts from 'evaluate' to 'hold the integrand exact and
         differentiable' -- the role it has played since Chapter 1's Cell 6.

The unit of progress in this book is the integral, and the discipline is always
the same: derive the anchor, run the regression, measure the estimator's error
in both currencies (variance and bias), and never claim a black canvas means
correctness.
""")
Expect: the roadmap rewritten as a family of integrals — each with an anchor, a regression, and two kinds of error.

§6 · The honest ledger — Chapter-7 deltas

DimensionChapter 6Chapter 7
Lighttwo point lightslight 1 is a disk; contribution = integral over the source (series units: the average; physical units derived in Cell 2)
Shadowslit mask per lightlit fraction f₁ ∈ [0,1]: the mask integrated; hard fill kept for contrast
Pipelinessame quantity, two formulationssame integrand, two estimators: jittered MC vs deterministic strata
Errorinterior max + flips+ interior RMS — the variance meter; speckle (numeric) vs banding (symbolic), both photographed (Cells 3, 7)
Exact anchorchecker box integralon-axis μ̄ = 2d(√(d²+ρ²)−d)/ρ², verified to 1e-8 (Cells 2, 6)
Regressionpoint filterρ→0 / K=1 returns Chapter 5 to machine precision (Cell 6)
Boundarypiecewise-algebraic pixel curvesthe penumbra's edge on the source is a conic — tangent cone × light plane, classified per point (Cells 2, 7c)
Determinismboth sides deterministicnumeric now has a seed; symbolic does not — and the ledger says why that matters (Cell 4)

§7 · Pitfalls gallery, continued

  1. Banding is bias with a face. Deterministic strata print their pattern through the mask into the penumbra (Cell 7a). Jittering trades it for speckle — and acquires a seed. There is no free lunch, only priced ones.
  2. The terminator is a mask too. max(0, μ) discontinuities in the derivative pick up the same stepping at low K; the penumbra is not the only edge that bands.
  3. Variance hides in RMS, not in max. A single max reading cannot distinguish speckle from structure; the chapter's ledger carries both numbers, and so should every renderer comparison you ever write.
  4. The conic can change regime. Ellipse vs hyperbola is a per-point classification (Cell 7c); code that assumes the umbra cone is always bounded is wrong somewhere on some surface.
  5. The ρ→0 regression is a contract, not a courtesy. Any "improvement" that does not degenerate to its predecessor in the documented limit is not an improvement; it is a different renderer wearing the same name.
  6. Determinism is a feature with a price. Reproducibility, regression tests, and differentiability argue for strata; unbiasedness argues for jitter. Choose openly; the ledger records the choice.

§8 · A brief history, continued

Soft shadows entered practice with Cook's distributed ray tracing (1984) — the same paper that framed depth of field and motion blur as one sampling doctrine — and with Brotman & Badler's penumbra volumes the same year. The hemisphere and form-factor machinery is older than graphics: Lambert wrote the cosine law in 1760; the radiosity tradition (Goral et al., 1984) turned it into linear algebra, and the disk-to-point form factors used in Cell 2 are its homework problems. The tangent-cone penumbra as a conic is classical projective geometry wearing a worklight — Apollonius knew the cone, and the computer graphics community rediscovered it in the 1990s as "soft shadow volumes." In the other world: stratified versus random quadrature is a century of numerical analysis, and Monte Carlo's 1/√K was priced by Ulam and Metropolis at Los Alamos in 1946 — the same decade the ENIAC rendered its first curves. This chapter merely observes that the two traditions have been integrating the same mask 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.