Symbolic Rendering

Chapter 1 — One Sphere, Two Computations

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–§3 set up the two pipelines and derive the math once, carefully. §4 is a live simulator: both renderers run in your browser, pixel by pixel, next to their difference image and the analytically derived silhouette circle. §5 is the laboratory: eight Python cells (sympy + numpy) to paste into Colab, which reproduce — and extend — everything the simulator shows. §6–§9 are the honest ledger, the pitfalls gallery, a little history, and the roadmap.
☰  contents
  1. Two ways to make a pixel
  2. The scene (Chapter-1 assumptions)
  3. The math, once, carefully
  4. Live simulator
  5. The Colab laboratory (Python / sympy)
  6. The honest ledger
  7. Pitfalls gallery
  8. A brief history (both of your worlds)
  9. Roadmap (measured steps)

§1 · Two ways to make a pixel

Numerical rendering computes a number per pixel by running an algorithm whose branches execute at runtime. Symbolic rendering computes a function once — a closed-form expression for the entire image, with the decision tree folded into piecewise mathematics — and then merely evaluates it per pixel. The image becomes a mathematical object you can inspect, differentiate, integrate, and re-parameterize.

scene θ = {R, e, L, k_a, k_d, k_s, n}→
NUMERIC: for each pixel: build ray → quadratic → if disc<0 → black else intersect → shade with max() clamps
→pixel value (float)
scene θ→
SYMBOLIC: derive once: I(x, y; θ) = Piecewise( I_Phong if Δ(x,y) ≥ 0, 0 ) → compile (lambdify → numpy) → for each pixel: evaluate the expression (branch = which sub-expression applies)
→pixel value (float)

Exactly as you framed it: the renderer's decision tree becomes standard mathematical logic — "if this situation, this expression applies; otherwise, that one" — carried by Piecewise, Max, and (seductively, dangerously) Heaviside.

§2 · The scene (Chapter-1 assumptions)

Grayscale intensity (one channel); one sphere of radius $R$ centered at the origin; pinhole camera at $E = (0,0,e)$ looking down $-z$; image plane $z = 0$, so pixel $(i,j)$ maps to world point $(x, y, 0)$ and ray direction $D = (x, y, -e)$ (unnormalized — a symbolic luxury); one point light at $L = (l_x, l_y, l_z)$; Phong 1975 illumination (ambient + diffuse + specular); background = 0 (black); no shadows — yet.

§3 · The math, once, carefully

Ray–sphere: $|E + tD|^2 = R^2$ is the quadratic $A t^2 + 2b t + c = 0$ with $A = x^2+y^2+e^2$, $b = -e^2$, $c = e^2-R^2$. The reduced discriminant (¼ of the textbook one) is

$$\Delta(x, y) = b^2 - A c = e^2R^2 - (e^2-R^2)(x^2+y^2)\qquad \text{hit} \iff \Delta \ge 0$$

Near root: $t = (e^2 - \sqrt{\Delta})/A$. Then the closed form the simulator and the Colab lab both use:

$$P = (t x,\ t y,\ e(1-t)), \quad N = P/R, \quad \hat w = \frac{L - P}{\|L - P\|}, \quad \hat V = \frac{(-x,\,-y,\,e)}{\sqrt A},$$ $$\mu = N\cdot\hat w, \qquad \nu = N\cdot\hat V = \frac{\sqrt{\Delta}}{R\sqrt A}, \qquad \hat r\cdot\hat V = 2\mu\nu - \hat w\cdot\hat V,$$ $$I(x, y; \theta) = \begin{cases} k_a + k_d\cdot\max(0, \mu) + k_s\cdot\max(0, \hat r\cdot\hat V)^n, & \text{if } \Delta(x,y) \ge 0\\[4pt] 0, & \text{otherwise}\end{cases}$$

Three gifts the symbolic route hands us for free — each invisible to a per-pixel number cruncher:

  1. The normal is exact and linear. Since $|P| = R$ by construction, $N = P/R$ with no normalization square root — the CAS can prove it; a float pipeline just hopes.
  2. The view direction never touches the hit point. $E - P = -t\cdot D$, so $\hat V = -D/\|D\|$: a function of the pixel alone. (The numeric/symbolic agreement in §5 validates this simplification.)
  3. The silhouette is a theorem. $\Delta$ depends only on $x^2+y^2$, so the sphere's image is exactly circularly symmetric, with outline $x^2+y^2 = e^2R^2/(e^2-R^2)$: $$\rho = \frac{eR}{\sqrt{e^2-R^2}} \qquad (e=4,\ R=1 \;\Rightarrow\; \rho = 4/\sqrt{15} \approx 1.03280).$$ Pixel rasterization approximates this curve; the expression is this curve.

Counterpoint: left alone, the CAS expands the Phong chain into a monster (measure it with sp.count_ops in Cell 2 — this is the classic expression swell of the Macsyma era). Human-guided simplification — the three gifts above — is what makes symbolic rendering practical. That division of labor is itself a finding.

§4 · Live simulator

One sphere, two computations: the algorithm run per pixel, vs. the derived expression evaluated per pixel — analytic silhouette ρ = eR/√(e²−R²) overlaid in red.

resolution 
initializing…
NUMERIC — run the algorithm per pixel (runtime ifs)
SYMBOLIC — evaluate I(x,y;θ) per pixel (piecewise math)
|difference|, amplified ×10¹⁵ — agreement is structural, not approximate
numeric: – symbolic: – max |ΔI|: – ρ = –
hover the canvases — the probe shows pixel → world (x,y), Δ, which branch of the piecewise applies, and I.

Try this: switch symbolic mode to Heaviside multiply (raw). The "one-formula" rendering I = H(Δ)·I_Phong fills the background with magenta NaNs. Symbolically, 0·anything = 0. In IEEE-754, 0·NaN = NaN — the miss branch computes √Δ of a negative number, and the zero no longer annihilates it. The regularized mode (√Δ → √max(0,Δ), a total expression) restores the identity. Symbolic and float semantics diverge at singularities — your first genuine symbolic-rendering pitfall, demonstrated live.

§5 · The Colab laboratory (Python / sympy)

Open a fresh notebook at colab.research.google.com, paste the cells below in order, run in order. Every cell carries three run paths on one line: 📋 Copy code to the clipboard, ⬇ notebook download (upload to Colab via File → Upload notebook), or paste into colab.new.

🐍 Cell 1 — the numerical baseline (per-pixel algorithm, branches at runtime)

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

def render_numeric(width, height, R=1.0, e=4.0, light=(5.0, 5.0, 10.0),
                   ka=0.1, kd=0.7, ks=0.5, shin=32.0):
    """Classic numerical renderer: per pixel, run the algorithm (with branches)."""
    img = np.zeros((height, width))
    E = np.array([0.0, 0.0, e])                       # pinhole camera
    L = np.array(light, dtype=float)                  # point light
    half = 1.25 * e * R / np.sqrt(e**2 - R**2)        # world half-extent of view
    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])                  # ray through pixel (unnormalized)
            A = float(D @ D)                          # quadratic coefficients
            b = float(E @ D)                          # = -e^2
            c = e * e - R * R
            disc = b * b - A * c                      # = e^2 R^2 - (e^2-R^2)(x^2+y^2)
            if disc < 0.0:                           # ---- decision tree: ray misses
                continue                              # background stays black
            t = (-b - np.sqrt(disc)) / A              # nearer intersection
            P = E + t * D                             # hit point
            N = P / R                                 # unit normal (|P| = R exactly)
            w = L - P
            w_hat = w / np.linalg.norm(w)             # toward light
            v_hat = (E - P) / np.linalg.norm(E - P)   # toward eye
            ndl = float(N @ w_hat)
            r = 2.0 * ndl * N - w_hat                 # mirror reflection of light
            rdv = float(r @ v_hat)
            img[j, i] = ka + kd * max(0.0, ndl) + ks * max(0.0, rdv) ** shin
    return np.clip(img, 0.0, 1.0)

t0 = time.perf_counter()
img_num = render_numeric(256, 256)
t1 = time.perf_counter()
print(f"numeric render: {t1 - t0:.3f} s")
plt.figure(figsize=(4, 4))
plt.imshow(img_num, cmap='gray', vmin=0, vmax=1)
plt.title("Numerical (per-pixel algorithm)"); plt.axis('off'); plt.show()
Expect: a lit sphere on black; the Python loop takes on the order of 0.3–1 s at 256² in Colab.

🐍 Cell 2 — the symbolic derivation (sympy does Chapter 1's algebra; watch expression swell)

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

# --- Symbols -----------------------------------------------------------------
x, y = sp.symbols('x y', real=True)                 # image-plane coordinates
R, e = sp.symbols('R e', positive=True)             # sphere radius, eye distance
lx, ly, lz = sp.symbols('l_x l_y l_z', real=True)   # point-light position
ka, kd, ks, n = sp.symbols('k_a k_d k_s n', positive=True)

# --- Ray and intersection ----------------------------------------------------
E = sp.Matrix([0, 0, e])          # pinhole
D = sp.Matrix([x, y, -e])         # ray through pixel (x, y, 0), unnormalized

A = D.dot(D)                      # x^2 + y^2 + e^2
b = E.dot(D)                      # -e^2
c = e**2 - R**2
disc = sp.expand(b**2 - A*c)      # hit test: disc >= 0  (reduced discriminant)
t_hit = sp.simplify((-b - sp.sqrt(disc)) / A)   # nearer root

# --- Hit point, normal, lighting vectors -------------------------------------
P = (E + t_hit * D).applyfunc(sp.simplify)
N = P / R                                     # exact unit normal, |P| = R
L = sp.Matrix([lx, ly, lz])
w = L - P
w_hat = w / sp.sqrt(w.dot(w))                 # toward light
v_hat = -D / sp.sqrt(A)                       # toward eye: E - P = -t*D  (gift #2)

ndl = sp.simplify(N.dot(w_hat))               # diffuse geometry term mu
r_vec = 2*ndl*N - w_hat
rdv = sp.simplify(r_vec.dot(v_hat))           # specular geometry term

# --- The rendered image as a symbolic object: the decision tree becomes math --
I_shade = ka + kd*sp.Max(0, ndl) + ks*sp.Max(0, rdv)**n
I_img = sp.Piecewise((I_shade, disc >= 0), (0, True))

print("disc   =", disc)
print("t_hit  =", t_hit)
print("N.V    =", sp.simplify(N.dot(v_hat)), "   (expect sqrt(disc)/(R*sqrt(A)))")
print("\n--- expression swell meter ---")
for name, expr in [("t_hit", t_hit), ("ndl", ndl), ("rdv", rdv),
                   ("I_shade", I_shade), ("I_img", I_img)]:
    print(f"  {name:8s}: {sp.count_ops(expr):4d} ops, {len(str(expr)):5d} chars")
Expect: a few seconds of one-time derivation cost (the simplify calls), then the swell table. This cost is paid ONCE — the compiled expression is then evaluated forever.

🐍 Cell 3 — compile the image to numpy, evaluate everywhere, verify against numeric

· · paste into colab.new → Run
params = (x, y, R, e, lx, ly, lz, ka, kd, ks, n)
f_img = sp.lambdify(params, I_img, modules='numpy')   # the image, compiled

SC = dict(R=1.0, e=4.0, light=(5.0, 5.0, 10.0), ka=0.1, kd=0.7, ks=0.5, shin=32.0)

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)

X, Y = grid(256, 256)
with np.errstate(invalid='ignore', divide='ignore'):
    t0 = time.perf_counter()
    img_sym = f_img(X, Y, SC['R'], SC['e'], *SC['light'],
                    SC['ka'], SC['kd'], SC['ks'], SC['shin'])
    t1 = time.perf_counter()
img_sym = np.clip(np.nan_to_num(img_sym), 0, 1)
print(f"symbolic render (derived once, evaluated per pixel): {t1 - t0:.4f} s")

img_num = render_numeric(256, 256, **SC)
diff = np.abs(img_num - img_sym)
print(f"max |numeric - symbolic| = {diff.max():.3e}   (structural agreement, not approximation)")

fig, axs = plt.subplots(1, 3, figsize=(12, 4))
for ax, im, ttl in zip(axs, [img_num, img_sym, diff],
                       ["Numerical", "Symbolic (compiled expression)", "|difference| (auto-scaled)"]):
    ax.imshow(im, cmap='gray'); ax.set_title(ttl); ax.axis('off')
plt.show()
Expect: identical images; max difference ~1e-15 (float operation ordering), and the symbolic eval in milliseconds. sympy's Piecewise lambdifies to np.select; Max to np.maximum — the decision tree is now numpy code.

🐍 Cell 4 — honest timings: "numeric vs symbolic" is orthogonal to "loop vs vectorized"

· · paste into colab.new → Run
def render_numeric_vec(width, height, R=1.0, e=4.0, light=(5., 5., 10.),
                       ka=.1, kd=.7, ks=.5, shin=32.):
    """Same math, numpy-vectorized: masks instead of ifs. Still 'numerical':
       no symbolic object is ever constructed or reused."""
    half = 1.25 * e * R / np.sqrt(e**2 - R**2)
    xs = np.linspace(-half + half/width,  half - half/width,  width)
    ys = np.linspace( half - half/height, -half + half/height, height)
    X, Y = np.meshgrid(xs, ys)
    A = X**2 + Y**2 + e**2
    disc = e**2 * R**2 - (e**2 - R**2) * (X**2 + Y**2)
    hit = disc >= 0
    sq = np.sqrt(np.maximum(disc, 0))
    t = (e**2 - sq) / A
    Px, Py, Pz = t*X, t*Y, e*(1 - t)
    Nx, Ny, Nz = Px/R, Py/R, Pz/R
    Lx, Ly, Lz = light
    wx, wy, wz = Lx - Px, Ly - Py, Lz - Pz
    wn = np.sqrt(wx**2 + wy**2 + wz**2)
    ndl = (Nx*wx + Ny*wy + Nz*wz) / wn
    va = np.sqrt(A)
    vx, vy, vz = -X/va, -Y/va, e/va
    rdv = 2*ndl*(Nx*vx + Ny*vy + Nz*vz) - (wx*vx + wy*vy + wz*vz)/wn
    img = np.where(hit, ka + kd*np.maximum(0, ndl) + ks*np.maximum(0, rdv)**shin, 0.0)
    return np.clip(img, 0, 1)

print(f"{'grid':>9} {'numeric loop':>14} {'numeric vec':>13} {'symbolic compiled':>19}")
for W in (64, 128, 256, 512):
    H = W; X, Y = grid(W, H)
    t0 = time.perf_counter(); render_numeric(W, H, **SC);     t1 = time.perf_counter()
    t2 = time.perf_counter(); render_numeric_vec(W, H, **SC); t3 = time.perf_counter()
    with np.errstate(invalid='ignore', divide='ignore'):
        t4 = time.perf_counter()
        f_img(X, Y, SC['R'], SC['e'], *SC['light'], SC['ka'], SC['kd'], SC['ks'], SC['shin'])
        t5 = time.perf_counter()
    print(f"{W:>4}x{H:<4} {t1-t0:>13.3f}s {t3-t2:>12.3f}s {t5-t4:>18.3f}s")

print("\nMoral: the symbolic advantage is NOT per-pixel speed (both compiled forms")
print("are vectorized float evaluation). The advantage is the OBJECT you keep:")
print("inspectable, differentiable, integrable, re-parameterizable.")
Expect: loop ≫ vectorized ≈ compiled-symbolic. The payoff of symbolic is elsewhere — next cells.

🐍 Cell 5 — symbolic superpower #1: the silhouette is a theorem, overlaid on the pixels

· · paste into colab.new → Run
# Solve disc = 0 symbolically. disc depends only on s = x^2 + y^2, which itself
# proves the image of a centered sphere is exactly circularly symmetric.
s = sp.symbols('s', positive=True)
disc_s = sp.expand(disc).xreplace({x**2: s - y**2})   # rewrite disc in terms of s
rho2 = sp.solve(sp.Eq(disc_s, 0), s)[0]
rho = sp.sqrt(sp.simplify(rho2))
print("silhouette: x^2 + y^2 =", rho2, "   =>   rho =", rho)
print("numeric check (e=4, R=1):", float(rho.subs({R: 1.0, e: 4.0})))  # 4/sqrt(15)

half = 1.25 * SC['e'] * SC['R'] / np.sqrt(SC['e']**2 - SC['R']**2)
fig, ax = plt.subplots(figsize=(4.5, 4.5))
ax.imshow(img_sym, cmap='gray', extent=[-half, half, -half, half])
ax.add_patch(plt.Circle((0, 0), float(rho.subs({R: SC['R'], e: SC['e']})),
                        color='red', fill=False, lw=1))
ax.set_title("exact analytic silhouette vs. rasterized pixels")
plt.show()
Expect: the red circle threading the jagged pixel boundary — the continuous truth the raster approximates.

🐍 Cell 6 — superpowers #2 and #3: exact image derivatives; partial evaluation (baking θ)

· · paste into colab.new → Run
# (a) Differentiate the IMAGE with respect to a scene parameter.
#     This is the seed of modern differentiable / inverse rendering,
#     obtained here for free, exactly. (We differentiate the unclamped field;
#     the clamps Max(0,.) contribute kinks, and the silhouette contributes a
#     Dirac-delta "visibility gradient" -- see pitfalls section.)
dI = sp.diff(ka + kd*ndl + ks*rdv**n, lx)
print("dI/dl_x :", sp.count_ops(dI), "ops (swell again -- differentiation amplifies)")
g = sp.lambdify(params, dI, modules='numpy')
X, Y = grid(256, 256)
discN = SC['e']**2 * SC['R']**2 - (SC['e']**2 - SC['R']**2) * (X**2 + Y**2)
with np.errstate(invalid='ignore', divide='ignore'):
    G = g(X, Y, SC['R'], SC['e'], *SC['light'], SC['ka'], SC['kd'], SC['ks'], SC['shin'])
G = np.where(discN >= 0, np.nan_to_num(G), 0.0)
plt.figure(figsize=(4.5, 4)); plt.imshow(G, cmap='RdBu'); plt.colorbar()
plt.title("∂I/∂l_x : brightness sensitivity to light x"); plt.axis('off'); plt.show()

# (b) Partial evaluation: bake the scene constants into the expression.
#     Moving the light = new constants, SAME derivation -- no re-derivation ever.
subs_scene = {R: SC['R'], e: SC['e'], lx: SC['light'][0], ly: SC['light'][1],
              lz: SC['light'][2], ka: SC['ka'], kd: SC['kd'], ks: SC['ks'], n: SC['shin']}
I_baked = I_img.subs(subs_scene)
print("ops: free-params:", sp.count_ops(I_img), " baked-scene:", sp.count_ops(I_baked))
f_baked = sp.lambdify((x, y), I_baked, modules='numpy')
with np.errstate(invalid='ignore'):
    img_baked = np.clip(np.nan_to_num(f_baked(X, Y)), 0, 1)
print("baked matches compiled render:", np.allclose(img_baked, img_sym, atol=1e-12))
Expect: a signed sensitivity field (red/blue), and a baked expression with slightly fewer ops. Inverse rendering = gradient descent on images like G.

🐍 Cell 7 — pitfall, live: the Heaviside temptation and IEEE 0·NaN

· · paste into colab.new → Run
# "One formula, no Piecewise": I = H(disc) * I_shade. Symbolically impeccable.
Hfun = sp.Heaviside(disc, 0)
I_heaviside = I_shade * Hfun
heavi = {'Heaviside': lambda u, h0: (u >= 0).astype(float)}
fH = sp.lambdify(params, I_heaviside, modules=['numpy', heavi])
X, Y = grid(256, 256)
with np.errstate(invalid='ignore'):
    bad = fH(X, Y, SC['R'], SC['e'], *SC['light'], SC['ka'], SC['kd'], SC['ks'], SC['shin'])
print("NaNs produced (miss branch: 0 * sqrt(negative)):", np.isnan(bad).sum(), "of", bad.size)
# Symbolically 0*anything = 0. IEEE-754: 0*NaN = NaN. The zero cannot annihilate
# a subexpression that is singular off-branch. The algebra and the floats diverge.

# Regularization: make every subexpression TOTAL, so the 0 can safely annihilate.
I_reg = I_heaviside.subs(sp.sqrt(disc), sp.sqrt(sp.Max(0, disc)))
fR = sp.lambdify(params, I_reg, modules=['numpy', heavi])
with np.errstate(invalid='ignore'):
    good = fR(X, Y, SC['R'], SC['e'], *SC['light'], SC['ka'], SC['kd'], SC['ks'], SC['shin'])
print("NaNs after regularization:", np.isnan(good).sum())
print("matches Piecewise render:", np.allclose(np.nan_to_num(good), img_sym, atol=1e-12))
Expect: tens of thousands of NaNs, then zero after regularization. (The §4 simulator demonstrates the same thing live in JavaScript.)

🐍 Cell 8 — where it hurts: branch-space growth (teaser for Chapter 2)

· · paste into colab.new → Run
# Add a second sphere on the z-axis. Even before shading, the CONDITION space
# multiplies: visibility = hit(A) and not occluded-by(B), etc.
R2, cz = sp.symbols('R_2 c_z', positive=True)
C2 = sp.Matrix([0, 0, cz])
b2 = (E - C2).dot(D)
c2 = (E - C2).dot(E - C2) - R2**2
disc2 = sp.expand(b2**2 - A*c2)

I_two = sp.Piecewise(
    (sp.Symbol('I_A'), (disc >= 0) & (disc2 < 0)),   # A visible, B missed
    (sp.Symbol('I_B'), disc2 >= 0),                    # B hit (naive: ignores occlusion order)
    (0, True))
print(I_two)
print("\nThe decision tree is now PART of the algebra. With k objects, visibility is a")
print("quantified statement (exists-t / for-all-blockers) -- not a closed form.")
print("Shadows are the same quantifier one level deeper. This is the frontier for")
print("symbolic rendering, and where measured escalation (Chapter 2+) begins.")
Expect: a nested Piecewise. Discussion: expression growth, occlusion ordering, and why shadows resist closed forms.

§6 · The honest ledger

DimensionNumerical renderingSymbolic rendering
Artifact per framean algorithm run W×H timesan expression derived once, evaluated W×H times
Branchingruntime if / continuePiecewise, Max, (carefully) Heaviside — logic inside the math
Scene parameters θbaked into the run; change θ → re-runremain symbolic; change θ → substitute into the same derivation (Cell 6b)
Exactnessrounding at every step; structure invisibleexact algebra until final float evaluation; structure provable (silhouette, symmetry)
Derivatives ∂I/∂θfinite differences or autodiff machineryexact, immediate (Cell 6a) — the seed of inverse/differentiable rendering; kinks at Max, Dirac deltas at silhouettes
Antialiasingsupersample and hopein principle, integrate I(x,y) over the pixel square symbolically — an exact box filter (Chapter 6)
Performancedecades of engineering, GPUs; loop slow, vectorized fastderivation slow (once, seconds); compiled eval ≈ vectorized numeric (Cell 4); raw sympy evaluation unusably slow — always lambdify
Scaling to scenesmore code, same asymptoticsbranch-space and expression-swell explosion; shadows/reflections are quantified/recursive — resist closed forms (Cell 8)
Verificationground truth hard to come bythe expression is an inspectable ground truth; numeric/symbolic agreement cross-validates both (Cell 3)
Failure modesepsilon hacks, z-fighting0·NaN at singularities, Max non-differentiability, assumption management in the CAS

§7 · Pitfalls gallery

  1. 0 · NaN ≠ 0. Symbolic annihilation fails on singular off-branch subexpressions. Cure: regularize to total subexpressions (√max(0,Δ)), or stay with Piecewise/np.select semantics. (§4 toggle, Cell 7.)
  2. Clamps are kinks. max(0, μ) is C⁰ but not C¹: ∂/∂θ exists a.e. but the terminator curve carries a kink; the silhouette carries a Dirac delta (a visibility gradient). Differentiate the unclamped field, then reason about the boundary — exactly what modern differentiable renderers formalize.
  3. Expression swell. Raw expansion of the Phong chain grows fast; differentiation amplifies it (measure, don't trust — count_ops). The three gifts of §3 are the human-guided simplifications that keep Chapter 1 tractable.
  4. Quantifiers are not closed forms. "Ray hits A before any blocker" is ∃/∀ logic. Symbolic rendering must either enumerate orderings (branch explosion) or change representation. This is the frontier, not a bug.

§8 · A brief history (both of your worlds)

Utah, 1968–77: Warnock's hidden-surface algorithm ('69), Gouraud shading ('71), Bui Tuong Phong's illumination model ('73 thesis / '75 CACM paper — the one used here), Newell's teapot ('75), Blinn's bump and reflection models ('77), Catmull's z-buffer and subdivision — the cradle of algorithmic rendering. Meanwhile, a few years earlier and a continent away: Project MAC's Macsyma (Moses, Martin et al., ~1968–82) put "symbolic manipulation" itself into the machine — the lineage through DOE Macsyma to today's Maxima, and culturally to sympy. This project simply reunites the two traditions: let the CAS hold the mathematics, let the renderer evaluate it. The modern echo is differentiable/inverse rendering (∂image/∂scene), which you can now see was implicit in the symbolic view from the start.

§9 · Roadmap (measured steps)

❦
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.