Mandelbrot Set Algorithms: Escape Time, Distance Estimation & GPU Rendering

1. The Definition: One Equation, Infinite Complexity

The Mandelbrot set M is the set of all complex numbers c for which the orbit of 0 under repeated application of f(z) = z² + c stays bounded forever:

f_c(z) = z² + c Starting value: z₀ = 0 Orbit: z₁ = 0² + c = c z₂ = c² + c z₃ = (c²+c)² + c ... c \in M ⟺ |z_n| does not \to \infty as n \to \infty Key fact: if |z_n| ever exceeds 2, the orbit escapes to infinity. Proof: if |z| > 2 and |c| < 2, then |z²+c| \geq |z|² - |c| > 2|z| - |z| = |z|. The modulus grows by at least |z| each step \to diverges.

The boundary of M — the infinitely intricate coastline visible in fractal zoom videos — is where the interesting algorithmic challenges live. Points strictly inside M iterate forever; points outside escape at some finite step. The boundary itself has Hausdorff dimension 2 (it is area-zero but so complex it fills the plane in a topological sense).

The set lives in the complex plane with real part (Re c) on the x-axis and imaginary part (Im c) on the y-axis. The full set fits within a circle of radius 2 centred at the origin, so we only need to sample c values with |c| ≤ 2.

2. Escape-Time Algorithm

The classic algorithm iterates each pixel's corresponding c value until |z|² > 4 (equivalent to |z| > 2, but avoids a square root) or until a maximum iteration count MAX_ITER is reached:

// Pseudocode — naive escape time function mandelbrot(cx, cy, max_iter): zx = 0.0, zy = 0.0 for n in 0 .. max_iter: if (zx*zx + zy*zy) > 4.0: return n // escaped — return iteration count zx_new = zx*zx - zy*zy + cx zy = 2.0*zx*zy + cy zx = zx_new return max_iter // did not escape \to "inside" the set

The returned value n is then mapped to a colour. The simplest mapping assigns a hue proportional to n / max_iter. But this produces harsh visible bands at every integer iteration count — which leads to the smooth colouring trick described in the next section.

Choosing MAX_ITER

At the default view (zoom ≈ 1×), 256–512 iterations suffice for recognisable detail. Deep zoom images require millions of iterations for points near the boundary. A rule of thumb: MAX_ITER ≈ 200 × (zoom level)^0.4. Too few iterations makes the set look "fat" — exterior points near the boundary incorrectly classified as interior.

Performance note: At 1080p resolution (1920×1080 ≈ 2 million pixels) with 1000 iterations each, a naive single-threaded CPU renderer performs ~2 billion floating-point multiplications per frame. A modern CPU takes ~2–5 seconds per frame. A GPU shader runs this in under 16 ms (60 fps) because each pixel is an independent computation — the definition of embarrassingly parallel workloads.

Optimisation: Skip the Square Root

Checking zx*zx + zy*zy > 4.0 avoids the expensive sqrt() call. Another early optimisation is to periodically check whether z has converged to a cycle — indicating the point is inside M — saving time on "definitely interior" points.

3. Smooth Colouring & Normalised Iteration Count

Integer escape counts produce colour bands — concentric rings clearly separated at every integer n. The normalised iteration count (also called the Büldt or Linas formula) interpolates between integers using the logarithm of the modulus at escape:

// After escape at iteration n, z has modulus |z| = r: nu = n - log2(log2(|z|)) = n - log2( log(|z|) / log(2) ) Interpretation: The "true" fractional escape time is the value of n at which |z| would have equalled exactly 2, had we allowed non-integer steps. Using log₂(log₂(|z|)) corrects for the quadratic growth of |z|. Smooth colour: colour(nu) — interpolate a colour palette by nu mod 1

This removes all visible banding, producing the smooth gradient colouring seen in professional fractal renders. Palette design then determines the final appearance: cycling through hue, using sinusoidal RGB functions, or mapping to hand-crafted gradient stops.

Orbit Traps

A more artistic colouring method records the minimum distance from the orbit to a geometric shape (a point, line, circle, or cross). The colour depends on this minimum distance. Orbit traps produce dramatically different visuals — repeating patterns of the trap shape across the fractal boundary — and can generate textures resembling real objects.

// Orbit trap: point trap at (tx, ty) min_dist = +\infty for each iteration (zx, zy): d = sqrt((zx-tx)² + (zy-ty)²) min_dist = min(min_dist, d) // Colour based on min_dist (e.g., log(min_dist) \to hue)

4. Interior Detection: Period Checking & Cardioid Test

Points inside M iterate all the way to MAX_ITER, making interior pixels the most expensive per-pixel computation. Two techniques dramatically reduce interior cost.

Cardioid and Period-2 Bulb Test

The main cardioid and the largest circular bulb (centred at c = −1) together contain a large fraction of all interior points. They can be checked analytically before iterating:

// Cardioid test q = (cx - 0.25)² + cy² in_cardioid = q * (q + (cx - 0.25)) < 0.25 * cy² // Period-2 bulb test in_p2_bulb = (cx + 1)² + cy² < 0.0625 if (in_cardioid || in_p2_bulb): return max_iter immediately

Periodicity / Cycle Checking

Interior orbits eventually enter a periodic cycle — the orbit repeats with some period p. Brent's or Floyd's cycle detection can identify this: save a "slow" copy of z every k iterations, compare with the current z. If |z_current - z_saved| < ε, the point is in M and iteration stops early.

// Brent-style period checking (simplified) zx_old = 0, zy_old = 0 period = 0, check = 0 for n in 0 .. max_iter: iterate z if |z - z_old| < 1e-10: return max_iter // cycle detected \to inside M check++ if check == period: zx_old = zx; zy_old = zy check = 0; period += period // double check interval

Combined, cardioid test + period checking reduces interior computation by 40–80% at typical zoom levels, roughly doubling average rendering speed.

5. Distance Estimation for Boundary Detail

Computing the exterior distance estimator (DEM) gives the approximate distance from a point c to the nearest boundary of M, without knowing exactly where the boundary is. This is achieved by simultaneously tracking the derivative of the iteration function:

// Distance estimation: track z and dz/dc simultaneously z = 0 + 0i, dz = 1 + 0i // dz is derivative of z w.r.t. c for n in 0 .. max_iter: dz = 2 * z * dz + 1 // chain rule: d/dc [z²+c] = 2z\cdot(dz/dc) + 1 z = z*z + c if |z|² > bailout²: break // Exterior distance estimate: d \approx |z| * log(|z|) / |dz| If d < pixel_size: the pixel lies on or near the boundary \to colour it.

Applications of DEM:

Anti-aliasing: Pixels closer to the boundary than one pixel width receive sub-pixel sampling (supersampling or Monte Carlo).
Edge rendering: Assign a distinct colour to all pixels where DEM < threshold — draws the set's outline without iterating to high depths.
3D relief rendering: Use DEM as a height field; shade with a virtual directional light to produce a "embossed" fractal relief map.

Accuracy: DEM is only an approximation valid for the exterior of M. Near Misiurewicz points (where the critical orbit is pre-periodic rather than periodic), the formula can over-estimate distances. In practice the estimate is close enough for rendering use.

6. Deep Zoom & Perturbation Theory

Zooming in to depths beyond ~10¹⁴ exhausts the precision of 64-bit double floating-point numbers (which have ~15–16 decimal digits of mantissa). At these depths, rounding errors accumulate and the image degenerates into incorrect solid-colour regions called glitches.

Arbitrary-Precision Arithmetic

The straightforward fix replaces doubles with arbitrary-precision floating-point (e.g., 128-bit, 256-bit, or even 1000+ digit libraries). But cost scales roughly quadratically with precision — 10× more digits → 100× slower. This makes interactive deep-zoom renders impossible without additional cleverness.

Perturbation Theory

The key insight: most pixels in a deep-zoom view are very close to each other in the complex plane. Compute one reference orbit at the viewport centre using high-precision arithmetic. Then represent every other pixel as a small perturbation δ from the reference:

Z = reference orbit (computed in high precision) z = Z + δ (pixel orbit, small perturbation) z_{n+1} = z_n² + c (Z + δ)_{n+1} = (Z + δ)_n² + (C + Δc) ≈ Z_n² + 2·Z_n·δ_n + δ_n² + C + Δc δ_{n+1} = 2·Z_n·δ_n + δ_n² + Δc Because |δ| ≪ 1, compute δ iterations in double precision. Only the reference Z needs high precision.

This reduces high-precision arithmetic to a single orbit per frame (the reference), while all other pixels use cheap double arithmetic. Glitches from approximation failure are detected and corrected by picking additional reference orbits near problematic pixels.

State-of-the-art Mandelbrot renderers (such as Kalles Fraktaler and Mandel Machine) combine perturbation theory with Series Approximation (SA), which skips entire blocks of early iterations analytically — enabling interactive zoom to depths of 10¹⁰⁰⁰ and beyond.

7. GPU Rendering with GLSL Shaders

Modern GPUs execute thousands of shader invocations in parallel. A fragment shader running one pixel per invocation maps perfectly onto the Mandelbrot iteration loop. Each pixel needs no data from any other pixel — no synchronisation required.

// GLSL fragment shader skeleton uniform vec2 u_center; // view centre in complex plane uniform float u_scale; // units per pixel uniform int u_max_iter; void main() { vec2 c = u_center + (gl_FragCoord.xy - u_resolution*0.5) * u_scale; vec2 z = vec2(0.0); float n = 0.0; for (int i = 0; i < MAX_ITER_CONST; i++) { if (dot(z, z) > 4.0) break; z = vec2(z.x*z.x - z.y*z.y + c.x, 2.0*z.x*z.y + c.y); n++; } // Smooth colouring float nu = n - log2(log2(length(z))); float t = nu / float(MAX_ITER_CONST); gl_FragColor = palette(t); }

The main constraint on GPU shaders is that the loop bound must be a compile-time constant in older GLSL versions. Packing too many iterations overflows the GPU's instruction cache or exceeds driver watchdog timeouts, causing screen resets. Practical GPU limits are 10 000–50 000 iterations per pixel at interactive frame rates.

Double Precision on the GPU

Consumer graphics cards offer hardware double-precision (64-bit) arithmetic at roughly 1/32 the throughput of single precision. For deep zoom rendering where doubles are required, the GPU advantage shrinks. A common workaround is double-float emulation — representing each double as a pair of floats using the Dekker–Knuth double-float algorithm. This regains double-precision accuracy at ~3–4× the cost of single precision, still much faster than CPU-side arbitrary precision.

WebGL Implementation Considerations

WebGL 1 requires extension OES_element_index_uint for 32-bit indices; no native double support.
WebGL 2 supports EXT_color_buffer_float for 32-bit float render targets (useful for accumulation passes).
For interactive demos, render at half resolution and upscale — halves total fragment invocations at modest quality cost.
Progressive refinement: render a coarse pass first, then refine. Users perceive the coarse preview as instant response.

8. Julia Sets: A Family of Relatives

For every complex number c, the corresponding Julia set J(c) is computed by the same iteration — but now c is fixed and z varies across the image:

// Julia set: c is a constant parameter, z is the variable function julia(zx, zy, cx, cy, max_iter): for n in 0 .. max_iter: if (zx² + zy²) > 4.0: return n zx_new = zx² - zy² + cx zy = 2\cdotzx\cdotzy + cy zx = zx_new return max_iter

The Mandelbrot and Julia sets are intimately linked: if c is inside M, J(c) is a connected fractal; if c is outside M, J(c) is a disconnected "dust" (Cantor set topology). The Mandelbrot set can literally be seen as a map of all possible Julia shapes — each point in M corresponds to a connected Julia set, and the shape of M at small scales mirrors nearby Julia sets.

This connection is exploited in interactive explorers: clicking on any point c in a Mandelbrot viewer instantly updates a side panel showing J(c). For aesthetic exploration, the c values along the boundary of M — the filaments and antennas — produce the most intricate Julia shapes.

Mandelbrot as a parameter space: Technically, the Mandelbrot set is the connectedness locus of the quadratic family f_c(z) = z²+c. Higher-degree families (z³+c, z⁴+c, …) produce Multibrot sets — similar in spirit, with (d−1)-fold rotational symmetry and correspondingly more intricate structure. The same escape-time and distance-estimation algorithms apply with minor modifications.