Author

Steve Trettel

Published

January 2026

1 Day 2: Fractals — Homework

1.1 Checkpoints

These verify you’ve understood the core material.

C1. Julia Mouse. Modify the Mandelbrot shader to render a Julia set controlled by the mouse. Fix c to come from iMouse.xy (mapped to \([-2, 2] \times [-2, 2]\)) and initialize z from the pixel position. Drag around and watch the Julia set change—connected structures shatter into dust as you cross outside the Mandelbrot set.

C2. Higher Powers. The iteration \(z^2 + c\) generalizes to \(z^n + c\). For \(z^3\), you can use repeated multiplication:

vec2 ccube(vec2 z) {
    return cmul(cmul(z, z), z);
}

For arbitrary powers, use polar form:

vec2 cpow(vec2 z, float n) {
    float r = length(z);
    float theta = atan(z.y, z.x);
    return pow(r, n) * vec2(cos(n * theta), sin(n * theta));
}

Try \(n = 3, 4, 5\)—what rotational symmetry does each have? Why? Try \(n = 2.5\) and look for the discontinuity.

C3. Animated Zoom. Animate a zoom into the Mandelbrot set. To zoom smoothly, the magnification factor \(m(t)\) should satisfy \(\frac{dm}{dt} = k \cdot m\) for constant \(k\)—that is, the rate of magnification is proportional to current magnification. The solution is \(m(t) = a^t\) for some base \(a > 1\).

Modify normalize_coord to incorporate zoom:

vec2 normalize_coord(vec2 fragCoord) {
    vec2 uv = fragCoord / iResolution.xy;
    uv = uv - vec2(0.5, 0.5);
    uv.x *= iResolution.x / iResolution.y;
    
    float zoom = pow(2.0, iTime);  // doubles every second
    return uv * 2.5 / zoom;
}

Then zoom into a specific point by adding an offset:

vec2 center = vec2(-0.745, 0.186);  // seahorse valley
vec2 c = center + normalize_coord(fragCoord);

Try different zoom rates (change the base from 2.0) and different center points. What happens to the boundary as you zoom past \(10^6\times\) magnification?

C4. Apollonian Animation. Animate the Apollonian gasket by letting the iteration count grow with time:

int max_iter = int(mod(iTime * 5.0, 50.0)) + 1;

Watch each iteration reveal a new layer of circles.

1.2 Explorations

These extend the day’s topics in interesting directions.

E1. Julia Animation. Animate c along the boundary of the main cardioid:

vec2 cardioid(float t) {
    vec2 eit = vec2(cos(t), sin(t));
    return cmul(eit, (vec2(2.0, 0.0) - eit)) / 4.0;
}
vec2 c = cardioid(iTime * 0.5);

Every point on the cardioid gives a Julia set with a parabolic fixed point—right at the edge of connectivity. What do you observe about the Julia sets as \(c\) traces this curve?

E2. Other Escape-Time Fractals. The Mandelbrot iteration \(z \mapsto z^2 + c\) is one choice among many. Implement these variations:

Burning Ship: \(z_{n+1} = (|\text{Re}(z_n)| + i|\text{Im}(z_n)|)^2 + c\) — take absolute values before squaring. Center around \((-0.5, -0.5)\) to see the main structure. Then try zooming into the small blob in the lower-left (around \((-1.8, -0.05)\) at 20× zoom) to find a miniature “burning ship.”
Tricorn: \(z_{n+1} = \bar{z}_n^2 + c\) where \(\bar{z} = x - iy\) is the complex conjugate. Also called the “Mandelbar.”
Celtic: \(z_{n+1} = (|\text{Re}(z_n^2)| + i\,\text{Im}(z_n^2)) + c\) — absolute value of real part after squaring.

The escape condition \(|z| > 2\) still applies. Experiment with different palette parameters to make these look good.

E3. Sierpinski Carpet. Render the Sierpinski carpet using the same iterate-and-test approach as the Apollonian gasket.

The carpet lives in the unit square \([0,1]^2\). At each iteration: 1. Scale coordinates by 3 to see which of 9 sub-cells you’re in 2. If you’re in the center cell, stop (this region is removed) 3. Otherwise, use fract to zoom into that cell and repeat

Color the carpet black and removed regions white (or vice versa). The key insight: floor(p * 3.0) gives the cell indices (0, 1, or 2 in each dimension), and fract(p * 3.0) gives position within that cell, ready for the next iteration.

E4. Sierpinski Triangle. Apply the same subdivision technique to the Sierpinski triangle.

We use a right triangle with vertices at \((0,0)\), \((1,0)\), and \((0,1)\) because the region tests are simple: - Bottom-left sub-triangle: \(x + y < 0.5\) - Bottom-right sub-triangle: \(x > 0.5\)
- Top sub-triangle: \(y > 0.5\) - Middle (removed): everything else

Each sub-triangle maps back to the full triangle by scaling by 2 and translating. Work out the three transformations and implement the shader.

Extension: The right triangle looks tilted. Apply an affine transformation at the end to map it to an equilateral triangle centered at the origin. The transformation is a shear followed by a scale.

E5. Apollonian Domain Coloring. Modify the Apollonian gasket to track which circle you inverted through most recently. After iteration, color by this: red for c1, green for c2, blue for c3, yellow for outer.

What structure do you see? Where are the boundaries between colored regions? How does this relate to the limit set?

1.3 Challenges

These require ideas beyond the lecture.

H1. Split-Screen Julia Explorer. Build a side-by-side display with the Mandelbrot set on the left and a Julia set on the right. The mouse position in the left half selects the parameter \(c\), and the Julia set updates in real time.

Requirements: 1. Left half shows Mandelbrot, right half shows Julia—each properly centered and aspect-corrected within its half 2. Mouse position determines \(c\), but only when clicking/dragging on the left half 3. A small colored dot marks the current \(c\) value on the Mandelbrot side 4. Both halves should use the same coloring scheme

The coordinate setup is the tricky part: each half of the screen needs its own independent mapping to \([-2, 2] \times [-2, 2]\) (or similar). Think carefully about aspect ratios for half-width rectangles.

Watch the Julia set transform as you drag across the Mandelbrot boundary—connected sets shatter into dust.

H2. Newton Fractal. Newton’s method finds roots of \(f(z) = 0\) by iterating \(z \mapsto z - f(z)/f'(z)\). For \(f(z) = z^3 - 1\): \[z_{n+1} = \frac{2z_n^3 + 1}{3z_n^2}\]

This converges to one of the three cube roots of unity: \(1\), \(e^{2\pi i/3}\), \(e^{4\pi i/3}\). Color each pixel by which root its orbit approaches.

You’ll need complex division:

vec2 cdiv(vec2 z, vec2 w) {
    float denom = dot(w, w);
    return vec2(z.x*w.x + z.y*w.y, z.y*w.x - z.x*w.y) / denom;
}

How do you determine which root the orbit converged to? The basin boundaries form a fractal—every point on the boundary has all three basins arbitrarily close by.

H3. Escape Radius Proofs. Prove the facts that justify the escape-time algorithm.

For the Mandelbrot set: 1. If \(|c| > 2\), then \(c \notin \mathcal{M}\). 2. If \(|z_n| > 2\) for some \(n\) (with \(|c| \leq 2\)), the orbit escapes to infinity.

For Julia sets: If \(|z| > R\) where \(R = \frac{1 + \sqrt{1 + 4|c|}}{2}\), then the orbit escapes.

Hint: The key inequality is \(|z^2 + c| \geq |z|^2 - |c|\) (triangle inequality). For what values of \(|z|\) is \(|z|^2 - |c| > |z|\)?

H4. Period Detection. Some orbits in the Mandelbrot set are periodic—they settle into cycles. A period-1 orbit converges to a fixed point; period-2 alternates between two values (\(z_{n+2} = z_n\)); and so on.

To detect periodicity, track the last few iterates and check if the current value matches one of them. You’ll need separate variables for \(z_{n-1}\), \(z_{n-2}\), etc. Check for matches only after enough iterations have passed for the orbit to settle (say, \(n > 50\)).

Implement a shader that colors the Mandelbrot set by period: one color for period-1, another for period-2, another for period-3. Points that escape get colored by iteration count as usual; points that don’t converge to a detected period stay black.

What structure do you find?

H5. Orbit Traps. The escape-time algorithm colors points by when the orbit escapes. An orbit trap colors points by where the orbit goes—specifically, how close it passes to a geometric “trap” shape.

During iteration, track the minimum distance from \(z\) to your trap. For points inside the Mandelbrot set (bounded orbits), color by this minimum distance. For points outside, just color black—the interesting structure is inside, where orbits pass near the trap repeatedly.

Good traps to try: - Point trap: length(z) — distance to origin - Circle trap: abs(length(z) - 1.0) — distance to unit circle - Cross trap: min(abs(z.x), abs(z.y)) — distance to axes

Experiment with coloring: palette(trap_dist * 5.0) or glow effects like exp(-trap_dist * 20.0).

The results reveal structure invisible in standard escape-time coloring.

1.4 Projects

Project 1: Grid of Julia Sets

Create a grid where each cell shows the Julia set for a different parameter \(c\). Use the grid formula from Day 1:

// Grid setup
float L = 0.5;                            // cell size

// Grid formula
vec2 cell_id = floor(p / L + 0.5);        // which cell (round to nearest)
vec2 cell_center = cell_id * L;           // center of this cell
vec2 local = p - cell_center;             // offset from center, in [-L/2, L/2]

Map cell_center to a parameter \(c\) in an interesting region of the complex plane (try the window \([-2, 0.5] \times [-1.25, 1.25]\)). Use local scaled to \([-2, 2]\) as the starting point \(z\) for the Julia iteration.

Color the Julia sets dark on a light background (inverted from the usual convention). When you zoom out, the grid reveals the Mandelbrot set: cells with \(c \in \mathcal{M}\) show connected Julia sets, cells with \(c \notin \mathcal{M}\) show scattered dust. This is the Mandelbrot-Julia correspondence made visible.

Project 2: Deep Zoom with Distance Estimation

The escape-time algorithm breaks down at high magnification: boundaries become blocky staircases of iteration bands. Increasing max_iter helps but never fully solves it. Distance estimation produces smooth, resolution-independent boundaries by computing the actual distance to the Mandelbrot set.

Part A: The Problem

Implement an animated zoom (C3) targeting the seahorse valley at \(c = -0.745 + 0.186i\). Watch the boundary as magnification passes \(10^3\), \(10^6\), \(10^9\). The curves degrade into chunky steps—this is fundamental to the escape-time method, not a bug you can fix by iterating more.

Part B: Distance to an Implicit Curve

For any implicit curve \(f(\mathbf{x}) = 0\), the distance from a nearby point to the curve is approximately \[d \approx \frac{|f(\mathbf{x})|}{|\nabla f(\mathbf{x})|}\]

This is the first-order Taylor approximation: \(f(\mathbf{x} + \boldsymbol{\epsilon}) \approx f(\mathbf{x}) + \nabla f \cdot \boldsymbol{\epsilon}\). Setting this to zero and solving for \(|\boldsymbol{\epsilon}|\) gives the formula above.

For the Mandelbrot set, we need the gradient of the escape-time function with respect to \(c\). Let \(z'_n = \frac{\partial z_n}{\partial c}\). Starting from \(z_0 = 0\) and \(z'_0 = 0\), differentiate \(z_{n+1} = z_n^2 + c\): \[z'_{n+1} = 2 z_n z'_n + 1\]

This is just the chain rule. Track \(z'\) alongside \(z\) during iteration.

When the orbit escapes, the distance to the Mandelbrot boundary is: \[d \approx \frac{|z_n| \log|z_n|}{|z'_n|}\]

The \(\log|z|\) factor comes from the potential theory of the Mandelbrot set (the Green’s function of the complement).

To read more, see

Distance Estimation Method for Fractals, Lindsay Robert Wilson

Distance to Fractals, Inigo Quilez

Part C: Implementation

Track z and dz as two vec2s. Important: update dz before z in the loop, since the derivative formula uses the current z. Use a large escape radius (256) for numerical stability in the log formula:

vec2 z = vec2(0.0, 0.0);
vec2 dz = vec2(0.0, 0.0);

for (int i = 0; i < 256; i++) {
    if (length(z) > 256.0) break;
    
    // Update dz BEFORE z (uses current z)
    dz = 2.0 * cmul(z, dz) + vec2(1.0, 0.0);
    z = cmul(z, z) + c;
}

float r = length(z);
float dr = length(dz);
float d = 0.5 * r * log(r) / dr;

Build a complete shader that renders the Mandelbrot set using distance estimation with a glow effect: exp(-d * scale). Start with the standard view (no zoom) and adjust scale until you see a bright glow near the boundary that fades to dark farther away. Protect against division by zero when dr is very small.

Part D: Adding Zoom

Add animated zoom to your Part C shader using the exponential zoom from C3. The challenge: as you zoom in, the glow becomes invisibly thin (the distance shrinks in world coordinates).

The fix: multiply d by the zoom factor before computing the glow. This keeps the visual thickness constant:

float zoom = pow(2.0, iTime);
// ... compute d as before ...
float scaled_d = d * zoom;
float glow = exp(-scaled_d * 200.0);

Target the seahorse valley at \(c = -0.745 + 0.186i\). Watch the boundary stay crisp as you zoom past \(10^6\times\) magnification—compare this to your C3 result.

Part E: Smooth Coloring for the Exterior

Distance estimation handles the boundary, but the exterior still shows discrete iteration bands. A fractional escape time removes them.

When the orbit escapes at iteration \(n\), the “true” escape time is slightly less than \(n\). The correction comes from the potential function \(G(z) = \lim_{n\to\infty} \frac{1}{2^n}\log|z_n|\):

// After the loop, if escaped:
float log_zn = log(length(z));
float nu = log(log_zn / log(2.0)) / log(2.0);
float smooth_iter = float(i) + 1.0 - nu;

Use smooth_iter instead of i for coloring. The bands vanish.

Part F: Putting It Together

Build a polished deep zoom combining: 1. Distance estimation for crisp boundaries at any magnification 2. Smooth coloring for the exterior 3. Animated exponential zoom (C3) 4. A visually interesting target (minibrots, spirals, dendrites)

Compare your result to pure escape-time rendering at the same magnification. The difference is dramatic.

Extension: Dual Complex Numbers

Instead of tracking z and dz separately, pack them into a vec4 where .xy is the value and .zw is the derivative. Define arithmetic so derivatives propagate automatically:

vec4 dcSqr(vec4 a) {
    return vec4(
        a.x*a.x - a.y*a.y,            // z²: real
        2.0*a.x*a.y,                   // z²: imag
        2.0*(a.x*a.z - a.y*a.w),      // (z²)': real
        2.0*(a.x*a.w + a.y*a.z)       // (z²)': imag
    );
}

vec4 dcAdd(vec4 a, vec4 b) {
    return a + b;
}

Then the iteration is just z = dcAdd(dcSqr(z), c) where c = vec4(c.xy, 1.0, 0.0).

Project 3: Orbit Visualization

Instead of coloring by iteration count, visualize the actual orbit. Let the mouse control the starting point \(z_0\):

Compute iterates \(z_0, z_1, \ldots, z_N\)
Draw a small circle at each iterate
Connect consecutive iterates with line segments (use the segment SDF from Day 1)
Color by iteration index—early iterates blue, later iterates red

Start with the Mandelbrot/Julia iteration: for a fixed \(c\) inside \(\mathcal{M}\), drag \(z_0\) around. Bounded orbits spiral toward an attractor; escaping orbits fly outward. You’ll see sensitive dependence—nearby starting points can have very different fates.

Once it works for complex dynamics, try the same visualization for the Apollonian gasket or Sierpinski carpet. Watch the point bounce between inversions or zoom into sub-cells. The iteration logic differs but the visualization technique is identical.