Chapter 2

Day 2: Fractals

2.1Overview

A fractal is a shape with structure at every scale—zoom in and you find more detail, forever. Today we make two of them:

The first is the Mandelbrot set, born from iterating $𝑧 \mapsto 𝑧 2 + 𝑐$ in the complex plane. The second is the Apollonian gasket, born from iterating circle inversions. Different mathematics, same principle: take a simple operation and repeat it.

Both fit naturally into shaders. Each pixel asks "what happens when I iterate from here?"—and since pixels don't depend on each other, the GPU can answer millions of these questions in parallel.

The Mandelbrot set lives in the complex plane, so we start there.

2.2Complex Numbers in GLSL

The complex numbers $ℂ$ are the plane $ℝ 2$ equipped with a multiplication. A point $(𝑎, 𝑏)$ represents $𝑎 + 𝑏 𝑖$ , and we define

(𝑎 + 𝑏 𝑖) (𝑐 + 𝑑 𝑖) = (𝑎 𝑐 - 𝑏 𝑑) + (𝑎 𝑑 + 𝑏 𝑐) 𝑖 .

This makes $ℂ$ a field, with $𝑖 2 = - 1$ encoding the rotation structure of complex multiplication.

2.2.1Representation

A complex number $𝑧 = 𝑎 + 𝑏 𝑖$ is naturally a vec2:

vec2 z = vec2(a, b);  // represents a + bi

We use the convention that z.x is the real part and z.y is the imaginary part.

2.2.2Arithmetic

Addition is componentwise, so GLSL's built-in + already does the right thing. No helper function needed.

Multiplication requires care. The * operator on vectors is componentwise—vec2(a,b) * vec2(c,d) gives vec2(a*c, b*d)—which is not complex multiplication. We implement the correct formula:

vec2 cmul(vec2 z, vec2 w) {
    return vec2(
        z.x * w.x - z.y * w.y,
        z.x * w.y + z.y * w.x
    );
}

The minus sign in the real part comes from $𝑖 2 = - 1$ .

2.2.3Magnitude

The magnitude $| 𝑧 | = \sqrt 𝑎 2 + 𝑏 2$ measures distance from the origin. GLSL provides a built-in length() function for this:

length(z)  // computes sqrt(z.x*z.x + z.y*z.y)

We'll use this to test conditions like $| 𝑧 | > 2$ .

2.3The Mandelbrot Set

Fix a complex number $𝑐$ . Starting from $𝑧 0 = 0$ , iterate:

𝑧 𝑛 + 1 = 𝑧 2 𝑛 + 𝑐

For some values of $𝑐$ , this sequence stays bounded. For others, it escapes to infinity. The Mandelbrot set $M$ is the set of all $𝑐$ for which the sequence remains bounded:

M = {𝑐 \in ℂ : t h e o r b i t o f 0 u n d e r 𝑧 \mapsto 𝑧 2 + 𝑐 i s b o u n d e d}

The boundary has intricate structure, but we never describe this geometry explicitly. Every pixel runs the same computation, asks "does my orbit escape?", and the structure emerges.

2.3.1The Escape Radius

We can't iterate forever, so we need a stopping criterion. Two facts make this practical:

Fact 1. If $| 𝑐 | > 2$ , then $𝑐 \notin M$ .

Fact 2. If $| 𝑧 𝑛 | > 2$ for some $𝑛$ , the orbit escapes to infinity, so $𝑐 \notin M$ .

Together, these justify the escape-time algorithm: iterate until $| 𝑧 | > 2$ (the orbit is escaping) or until we hit a maximum iteration count (the orbit is probably bounded). The number 2 here is the escape radius.

Proving these facts requires some careful estimates with the triangle inequality; we leave this as Challenge H4.

2.3.2Implementation

vec2 cmul(vec2 z, vec2 w) {
    return vec2(z.x * w.x - z.y * w.y, z.x * w.y + z.y * w.x);
}

vec2 normalize_coord(vec2 fragCoord) {
    vec2 uv = fragCoord / iResolution.xy;
    uv = uv - vec2(0.5, 0.5);
    uv.x *= iResolution.x / iResolution.y;
    return uv * 2.5;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 c = normalize_coord(fragCoord);
    c.x = c.x - 0.5;  // shift left to center the interesting part
    
    vec3 color = vec3(0.0, 0.0, 0.0);
    
    vec2 z = vec2(0.0, 0.0);
    
    for (int i = 0; i < 100; i++) {
        if (length(z) > 2.0) {
            color = vec3(1.0, 1.0, 1.0);
            break;
        }
        z = cmul(z, z) + c;
    }
    
    fragColor = vec4(color, 1.0);
}

2.4Coloring Escape-Time Fractals

Black and white shows the set, but we're throwing away information. The iteration count tells us how quickly a point escapes—points that escape after 5 iterations are different from points that escape after 50.

2.4.1Grayscale

Map the iteration count to brightness. To use the iteration count after the loop, declare it beforehand:

vec3 color = vec3(0.0, 0.0, 0.0);

vec2 z = vec2(0.0, 0.0);
int i;
for (i = 0; i < 100; i++) {
    if (length(z) > 2.0) break;
    z = cmul(z, z) + c;
}

if (i < 100) {
    float t = float(i) / 100.0;
    float gray = 1.0 - t;
    color = vec3(gray, gray, gray);
}

Structure appears: tendrils, spirals, bulbs around the boundary. Points far from the Mandelbrot set escape quickly (white); points near the boundary take many iterations to escape (gray); points in the set never escape (black).

2.4.2Color

Grayscale reveals structure, but color can reveal more. We use a cosine palette (see Appendix: Color) to map iteration counts to smooth color gradients:

vec3 palette(float t) {
    vec3 a = vec3(0.5, 0.5, 0.5);
    vec3 b = vec3(0.5, 0.5, 0.5);
    vec3 c = vec3(1.0, 1.0, 1.0);
    vec3 d = vec3(0.00, 0.33, 0.67);
    return a + b * cos(6.28318 * (c * t + d));
}

// in mainImage, after the loop:
if (i < 100) {
    float t = float(i) / 100.0;
    color = palette(t);
}

The color bands correspond to iteration counts. There's another fractal hiding in the same iteration.

2.5Julia Sets

The Mandelbrot set asks: for which $𝑐$ does the orbit of $0$ stay bounded? We can ask a different question: for a fixed $𝑐$ , which starting points $𝑧 0$ have bounded orbits?

Fix $𝑐 \in ℂ$ . The filled Julia set $𝐾 𝑐$ is

𝐾 𝑐 = {𝑧 0 \in ℂ : t h e o r b i t o f 𝑧 0 u n d e r 𝑧 \mapsto 𝑧 2 + 𝑐 i s b o u n d e d}

Same iteration, different question. For the Mandelbrot set, we vary $𝑐$ and always start at $𝑧 0 = 0$ . For a Julia set, we fix $𝑐$ and vary the starting point.

The code change is minimal—just swap which variable comes from the pixel and which is fixed. The iteration and coloring logic stay the same:

// Mandelbrot: c comes from pixel position, z starts at 0
vec2 c = normalize_coord(fragCoord);
c.x = c.x - 0.5;
vec2 z = vec2(0.0, 0.0);

// Julia: c is fixed, z comes from pixel position
vec2 c = vec2(-0.7, 0.27015);  // a fixed parameter
vec2 z = normalize_coord(fragCoord);

The escape condition $| 𝑧 | > 2$ still works for most Julia sets; the precise escape radius depends on $| 𝑐 |$ (see Exercise H3).

Different values of $𝑐$ produce very different Julia sets. Some are connected (one piece), others are totally disconnected (Cantor dust). The relationship is governed by a precise theorem.

2.5.1The Mandelbrot-Julia Correspondence

Each point $𝑐$ in the complex plane has an associated Julia set $𝐾 𝑐$ , and the Mandelbrot set tells you what kind:

Theorem (Douady-Hubbard). The filled Julia set $𝐾 𝑐$ is connected if and only if $𝑐 \in M$ .

If $𝑐$ lies inside the Mandelbrot set, the Julia set is a single connected piece—perhaps with a complicated boundary, but topologically one blob. If $𝑐$ lies outside the Mandelbrot set, the Julia set shatters into uncountably many pieces, a Cantor set of dust. The Mandelbrot set is precisely the set of parameters for which the Julia set holds together.

Points on the boundary of $M$ give the most intricate Julia sets—neither solid nor completely fractured, but balanced at the edge of connectivity.

Drag from inside the Mandelbrot set to outside and watch the Julia set transform. Connected structures with complicated boundaries give way to scattered dust as you cross into the exterior.

2.6Circle Inversion

We've been iterating a polynomial. But polynomials aren't the only thing we can iterate.

Circle inversion is to circles what reflection is to lines: it swaps inside and outside, preserves angles, and applying it twice returns you to where you started. We'll build another fractal the same way we built the Mandelbrot set: iterate a map, check a condition on where the point ends up, and color accordingly. For Mandelbrot, we checked whether the orbit escaped past radius 2. For our next fractal, we'll check which region the orbit lands in after bouncing between inversions.

2.6.1Definition

Inversion in the unit circle sends a point $𝐩$ to:

i n v (𝐩) = 𝐩 | 𝐩 | 2

The inverted point lies on the same ray from the origin, but at reciprocal distance: if $𝐩$ is at distance $𝑟$ , its image is at distance $1 / 𝑟$ . Points inside the circle map outside; points outside map inside; points on the circle stay fixed. (The origin maps to "infinity"—inversion extends naturally to the Riemann sphere.)

vec2 invert(vec2 p) {
    return p / dot(p, p);
}

2.6.2Visualizing Inversion

To see what inversion does, let's draw some shapes and watch them transform. The shader below toggles between original and inverted coordinates:

vec2 p_inv = invert(p);

// Toggle every second
vec2 q;
if (fract(iTime * 0.5) < 0.5) {
    q = p;
} else {
    q = p_inv;
}

// Draw shapes using q

Lines not through the origin become circles through the origin. Circles map to circles, or to lines if they pass through the center of inversion.

GLSL Shortcuts: mix and step

The toggle logic can be written more compactly:

step(edge, x) returns 0 if x < edge, otherwise 1
mix(a, b, t) linearly interpolates: returns a when t = 0, b when t = 1

float t = step(0.5, fract(iTime * 0.5));
vec2 q = mix(p, p_inv, t);

Invert a grid. The function mod(q, 0.5) gives the position of q within a repeating $0.5 \times 0.5$ cell. When either component is near zero, we're on a grid line:

vec2 grid = mod(q, 0.5);
if (grid.x < 0.02 || grid.y < 0.02) color = vec3(1.0, 1.0, 0.0);

The rectilinear grid becomes a web of circles, all passing through the origin.

2.6.3General Circle Inversion

So far we've inverted through the unit circle at the origin. For a circle with center $𝐜$ and radius $𝑅$ :

i n v (𝐩) = 𝐜 + 𝑅 2 𝐩 - 𝐜 | 𝐩 - 𝐜 | 2

The point is reflected through the circle: same ray from the center, reciprocal distance (scaled by $𝑅 2$ ).

2.7Structs

To invert through an arbitrary circle, we need to pass both a center and a radius. We could write:

vec2 invert(vec2 p, vec2 center, float radius) { ... }

But when working with multiple circles, this gets unwieldy. GLSL lets us bundle related data into a struct:

struct Circle {
    vec2 center;
    float radius;
};

Now Circle is a type. We can create instances and access their fields:

Circle c = Circle(vec2(1.0, 0.5), 0.7);
// c.center is vec2(1.0, 0.5)
// c.radius is 0.7

Our inversion function becomes cleaner:

vec2 invert(vec2 p, Circle c) {
    vec2 d = p - c.center;
    return c.center + c.radius * c.radius * d / dot(d, d);
}

And we can write helper functions that take circles as arguments:

float distToCircle(vec2 p, Circle c) {
    return abs(length(p - c.center) - c.radius);
}

bool isInside(vec2 p, Circle c) {
    return length(p - c.center) < c.radius;
}

2.8The Apollonian Gasket

The Apollonian gasket is a fractal circle packing, named for Apollonius of Perga who studied tangent circles in the 3rd century BCE. Start with four mutually tangent circles, three inside one, and fill each curved gap with a circle tangent to its three neighbors. Repeat.

2.8.1Setup

Place three circles of radius $𝑟$ with centers at the vertices of an equilateral triangle, all tangent to each other and enclosed by an outer circle tangent to all three:

float r = 1.0;
float circumradius = 2.0 * r / sqrt(3.0);  // center-to-vertex distance

Circle c1 = Circle(vec2(0.0, circumradius), r);
Circle c2 = Circle(vec2(-circumradius * sqrt(3.0)/2.0, -circumradius * 0.5), r);
Circle c3 = Circle(vec2(circumradius * sqrt(3.0)/2.0, -circumradius * 0.5), r);
Circle outer = Circle(vec2(0.0, 0.0), circumradius + r);

The gaps between circles are curvilinear triangles.

2.8.2Iteration

If a point is inside one of the inner circles, invert it through that circle, pushing it out. If it's outside the outer circle, invert through the outer circle—pulling it in. Repeat until the point lands in a gap (inside the outer circle but outside all inner circles) or we hit a maximum iteration count.

int i;
for (i = 0; i < 50; i++) {
    if (isInside(p, c1)) {
        p = invert(p, c1);
    } else if (isInside(p, c2)) {
        p = invert(p, c2);
    } else if (isInside(p, c3)) {
        p = invert(p, c3);
    } else if (!isInside(p, outer)) {
        p = invert(p, outer);
    } else {
        break;  // in a gap—done
    }
}

Color by iteration count, just like escape-time fractals:

float t = float(i) / 50.0;
vec3 color = palette(t);

Each inversion maps the configuration into a smaller copy of itself, creating self-similar structure at every scale.

2.8.3Full Implementation

struct Circle {
    vec2 center;
    float radius;
};

vec2 invert(vec2 p, Circle c) {
    vec2 d = p - c.center;
    return c.center + c.radius * c.radius * d / dot(d, d);
}

bool isInside(vec2 p, Circle c) {
    return length(p - c.center) < c.radius;
}

vec3 palette(float t) {
    vec3 a = vec3(0.5, 0.5, 0.5);
    vec3 b = vec3(0.5, 0.5, 0.5);
    vec3 c = vec3(1.0, 1.0, 1.0);
    vec3 d = vec3(0.00, 0.33, 0.67);
    return a + b * cos(6.28318 * (c * t + d));
}

void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 uv = fragCoord / iResolution.xy;
    uv = uv - vec2(0.5, 0.5);
    uv.x *= iResolution.x / iResolution.y;
    vec2 p = uv * 6.0;
    
    // Setup circles
    float r = 1.0;
    float circumradius = 2.0 * r / sqrt(3.0);
    
    Circle c1 = Circle(vec2(0.0, circumradius), r);
    Circle c2 = Circle(vec2(-circumradius * sqrt(3.0)/2.0, -circumradius * 0.5), r);
    Circle c3 = Circle(vec2(circumradius * sqrt(3.0)/2.0, -circumradius * 0.5), r);
    Circle outer = Circle(vec2(0.0, 0.0), circumradius + r);
    
    // Iterate inversions
    int i;
    for (i = 0; i < 50; i++) {
        if (isInside(p, c1)) {
            p = invert(p, c1);
        } else if (isInside(p, c2)) {
            p = invert(p, c2);
        } else if (isInside(p, c3)) {
            p = invert(p, c3);
        } else if (!isInside(p, outer)) {
            p = invert(p, outer);
        } else {
            break;
        }
    }
    
    // Color by iteration count
    float t = float(i) / 50.0;
    vec3 color = palette(t);
    
    fragColor = vec4(color, 1.0);
}

2.8.4The Limit Set

Points that land in a gap quickly are dark. Points near the fractal boundary, the limit set, take many iterations to settle and appear bright. We can emphasize the limit set with nonlinear coloring:

float t = float(i) / 100.0;
float t2 = pow(t, 2.0);
vec3 color = 30.0 * vec3(t2, t2, t2);

The squaring suppresses low iteration counts while the factor of 30 boosts high ones:

All of this from iterating four circle inversions.

2.9Summary

Today we built fractals from two different iterations: complex quadratic polynomials (Mandelbrot and Julia sets) and circle inversions (Apollonian gasket). Despite their different origins, both share the same algorithmic structure:

Iterate a transformation
Test a stopping condition (escape, or landing in a fundamental region)
Color based on iteration count

This pattern—per-pixel iteration with no communication between pixels—is ideal for GPUs.

2.9.1GLSL Skills

Complex arithmetic: Representing $ℂ$ as vec2, implementing cmul
Structs: Bundling related data (Circle with center and radius)
Control flow: for loops with early break, nested if/else
Built-in functions: dot, length, mod, step, mix
Cosine palettes: Mapping scalar values to smooth color gradients

2.9.2Key Concepts

The escape radius lets us stop iterating early—once $| 𝑧 |$ exceeds the threshold, we know the orbit escapes
The Mandelbrot set indexes Julia sets: $𝑐 \in M$ if and only if $𝐾 𝑐$ is connected
Circle inversion generalizes reflection and preserves angles (it's a conformal map)
Iteration count encodes geometric information—how close a point is to the fractal boundary