Computing fractals

An escape-time fractal (Chapter 15) is computed one pixel at a time, and each pixel is a loop: take z, square it, add c, check whether it has escaped, repeat. The total work is roughly the number of pixels times the average number of trips around that loop.

Both numbers are bigger than they look. A 1920 × 1080 frame has 2,073,600 pixels. Points inside the set never escape, so every one of them runs all the way to the iteration limit; with a limit of 1,000, a frame that is mostly inside costs about two billion iterations. And the limit cannot stay at 1,000. The deeper you zoom, the longer points near the boundary take to make up their minds, so deep zooms routinely need tens of thousands of iterations, sometimes millions.

One iteration of z² + c is cheap: with x and y the two parts of z, it is x² − y² + a and 2xy + b, which is three multiplications and a handful of additions if you keep x² and y² around for the escape test too. A modern processor core does a few billion simple operations a second. So a single deep frame is seconds of work on one core, and an animation at 60 frames a second is impossible by brute force. Everything that follows is a way of not doing the brute force.

The most expensive pixels are the ones inside the set, because they always run to the limit. Two tricks deal with most of them.

The cardioid and bulb check. The big heart-shaped body of the Mandelbrot set and the round bulb to its left have exact formulas. Before iterating at all, test whether c is inside either one. If it is, the answer is "inside" and the loop never runs. In a view of the whole set, that removes most of the black area for the price of a few multiplications per pixel.

Periodicity checking. Inside the set, most orbits do not wander forever; they settle into a repeating cycle. If the orbit ever comes back to a value it has already had (within a tiny tolerance), it will repeat forever and can never escape. So save z every so often and compare. A common scheme saves it at iterations 1, 2, 4, 8, 16 and so on, which catches a cycle of any length without storing the whole orbit (it is the same idea as Brent's cycle-finding method). Deep inside a bulb the loop stops after a few dozen iterations instead of thousands.

Two more tricks work at the level of the picture rather than the pixel. The Mandelbrot set is symmetric about the real axis, so a view that straddles it only needs half its rows computed. And because the set is connected, a rectangle whose whole border is inside the set is (in exact arithmetic) inside throughout; the Mariani–Silver method computes borders first and fills rectangles whose borders agree, subdividing those that do not. That last one is a gamble at finite resolution: a filament thinner than a pixel can slip through a border unseen, so careful renderers use it only with generous safety margins.

Computers store most real numbers in IEEE 754 floating point: a sign, an exponent and a fixed number of significant binary digits, the significand. It is scientific notation in base 2. The exponent lets the number be enormous or tiny; the significand decides how many digits of it are actually kept.

The two common sizes are float32 (single precision), with 24 significant bits, about 7 decimal digits, and float64 (double precision), with 53 significant bits, about 16 decimal digits. Graphics cards are built around float32. Ordinary programming languages, JavaScript included, use float64 by default.

A fixed number of digits means a fixed relative spacing between neighbouring numbers. Near 1, float32 can represent 1 and the next number up, 1.00000012, and nothing in between; that gap is 2⁻²³, about 1.2 × 10⁻⁷. For float64 the gap near 1 is 2⁻⁵², about 2.2 × 10⁻¹⁶. The gap is called the ulp, a unit in the last place.

Now zoom. A view 1,000 pixels wide and w units across puts adjacent pixels w/1,000 apart. While that is much bigger than the ulp, every pixel gets its own coordinate. Once it falls below the ulp, several adjacent pixels round to the same coordinate, compute the same orbit, and get the same colour. The image turns into blocks, and nothing about the iteration count can fix it.

Plug in numbers. For float32 near |c| ≈ 1 and a 1,000-pixel image, the view runs out of digits at a width of about 10⁻⁴, a magnification of a few times 10⁴ compared with the whole set. Where the coordinates are a bit smaller than 1, like the point in the instrument below, it holds out to around 10⁵ to 10⁶. For float64 the same arithmetic gives a width of about 10⁻¹³: magnifications around 10¹³ to 10¹⁵. That is the wall every simple deep-zoom program hits, and it is why famous zoom videos that go to 10¹⁰⁰ or 10¹⁰⁰⁰ need something else.

The straightforward way past 10¹⁴ is to use more digits. Arbitrary-precision libraries such as GMP and MPFR store a number as a long array of machine words and do the arithmetic in software, as many digits as you ask for. A zoom to 10¹⁰⁰ needs a little over 100 decimal digits, roughly 340 bits, just to tell neighbouring pixels apart.

It works, and it is slow. Each operation now loops over many words, and the cost of a multiplication grows faster than the number of digits (the schoolbook method is quadratic; cleverer methods help only for very long numbers). Deep zooms also need more iterations, so the two costs multiply.

Then K. I. Martin published a method, with a free program called SuperFractalThing, that changed the field. The insight is that all the pixels in a deep view are extraordinarily close together, so their orbits are almost the same. Compute one orbit, the reference orbit, at full arbitrary precision for a point C in the view. For every other pixel c = C + δ, do not compute its orbit; compute only how far it is from the reference, ε, which is tiny, and which can be tracked in ordinary float64.

Each pixel now costs about the same as a float64 pixel, plus one expensive orbit shared by the whole image. Martin's paper added a second trick, series approximation: for the first many iterations, ε is a smooth function of δ and can be written as a short polynomial, εₙ ≈ Aₙδ + Bₙδ² + Cₙδ³, whose coefficients are computed once alongside the reference. Every pixel can jump straight to iteration, say, 40,000 by evaluating the polynomial, and only iterate from there. For deep zooms that skips the vast majority of the work.

The price is glitches. The approximation assumes the pixel's orbit behaves like the reference's. When it does not (typically when the pixel's full value zₙ = Zₙ + εₙ passes much closer to zero than the reference does), the float64 ε loses its precision, and a whole patch of pixels comes out as a flat blob or a wrong copy of something nearby. Renderers detect suspicious pixels with a size test and redo them with a new reference point chosen inside the glitch. Newer renderers avoid most glitches by "rebasing": whenever a pixel's orbit gets closer to zero than its difference from the reference, restart that pixel against the reference orbit from its beginning.

One more wall remains. float64's exponent stops at about 10⁻³⁰⁸, so past a zoom of 10³⁰⁸ even the tiny δ values underflow to zero. Deep-zoom programs handle that by storing δ and ε with a separate, wider exponent: ordinary float64 digits with an extra integer saying how many powers of two to shift.

Every pixel of an escape-time fractal is independent of every other: no pixel needs another's answer. Programmers call that embarrassingly parallel, and it is a perfect fit for a graphics card, which runs thousands of small threads at once. Splitting rows between processor cores helps too, and so do SIMD instructions that do four or eight pixels' arithmetic in one go.

GPUs have two catches for fractals. First, threads run in lockstep groups (32 threads in an NVIDIA "warp", 32 or 64 on AMD). A group finishes only when its slowest thread does, so one pixel that needs 10,000 iterations holds its 31 neighbours idle while they wait. Fractals are the worst case for this, since the boundary mixes fast and slow pixels at every scale.

Second, GPUs are float32 machines. Many consumer cards run float64 at a small fraction of their float32 speed, often 1/32 or 1/64. Apple's Metal shading language has no double-precision type at all, and WebGL shaders have only 32-bit floats. That is why GPU fractal explorers dissolve into blocks around 10⁵ unless they emulate more precision in software (two float32s as a "double-float", the GPU cousin of double-double) or, better, use perturbation, where the per-pixel work is small numbers that low precision can handle.

A real-time visualiser has a hard budget. At 60 frames a second each frame gets 16.7 milliseconds, on a 120 Hz display 8.3, and the fractal has to share that with everything else on screen. When a full picture costs more than the budget, the work has to be spread across frames. That is progressive rendering, and there are two classic ways to do it.

• Coarse to fine. Compute every eighth pixel and draw it as a big block, then every fourth, and so on, as the 3D viewer in Chapter 14 does. The PNG image format's Adam7 interlacing works on the same principle, in seven passes. You see the whole picture at once, blurry, then sharpening.
• Bands. Compute a band of rows per frame, top to bottom, and keep sweeping. Simple, cache-friendly, and the cost per frame is fixed by the band height.

For a still image either is fine. For an animated fractal, where the parameters change continuously, banding has a trap. The picture on screen is a patchwork: rows finished just now, rows finished one sweep ago, rows in between. If the display also fades between redraws (a common choice in visualisers, where trails look good on other patterns), then the freshest band is bright and everything older is dimmer, and the eye reads that as a bright stripe rolling down the screen, several times a second. It looks like flashing. It is not a bug in the fractal; it is a bug in how the picture is assembled.

The fix is to stop showing work in progress. Build each sweep off-screen, in a buffer the viewer never sees. When a sweep finishes, it becomes the newest complete picture. On screen, show a crossfade between the last two complete pictures, moving linearly from the older to the newer over the time the next sweep takes. Every frame changes by a small, even amount everywhere, so motion reads as flow instead of a moving seam. The cost is latency: what you see is about one sweep behind what is being computed.

Anti-aliasing belongs in the same budget. A fractal boundary has detail finer than any pixel, so one sample per pixel gives jagged, sparkling edges that crawl when animated. Supersampling computes several samples per pixel and averages them; n × n samples cost n² times as much, so 4 × 4 is sixteen times the work. Renderers save most of that by sampling adaptively, only where neighbouring pixels disagree, and by jittering sample positions so that the error looks like fine grain rather than stair-steps. In real time, a cheaper cousin is to jitter the sample point a little each frame and let the temporal crossfade average them.

JellyDazzle is J's audio-reactive kaleidoscope app for Mac, iPhone and iPad, on the App Store. It runs 24 engine modes written in ARM64 assembly and 638 patterns written in C, and about 20 of those patterns are fractals: a Mandelbrot deep zoom, a Julia morph, a Burning Ship dive, Newton basins, Phoenix, Tricorn, Lyapunov, a Buddhabrot, a flame IFS, an Apollonian gasket, a growing fern, a dragon-curve fold and Sierpinski rain. Everything in this chapter shows up in it, because a visualiser lives inside the real-time budget all the time.

Bands per frame

JellyDazzle's fractals render progressively: a band of rows per frame, sweeping down the screen. That is the second strategy above, chosen for the usual reason: a full fractal frame does not fit in a frame's time alongside everything else the app is doing.

The 3.6.0 flash

In version 3.6.0, that made the fractals flash. A bright band rolled down the screen 6 to 12 times a second, because the canvas decayed between redraws: the freshest rows were bright, the older ones had faded, and the boundary between them swept down the picture over and over. It is precisely the "bands on a fading canvas" row of the table, and the lab above reproduces it.

The 3.6.1 fix

Version 3.6.1 fixed it with the two ideas from this chapter: build each sweep off-screen, and crossfade the last two finished pictures linearly, so the shape flows instead of flashing. On the fractals that were measured, the frame-to-frame change dropped from about 9–16 to under 1. That is more than a tenfold drop in how much the image changes between frames, which is the difference between a picture that strobes and one that breathes.

Seeds

Each appearance of a fractal pattern is driven by a seed, a number that picks its parameters: which Julia constant, which part of the Burning Ship to dive toward, which flame variations. A seed is what makes the twentieth Julia morph of an evening different from the first. In 3.6.5 an audit found patterns that ignored their seed, so they looked identical every time they appeared, and fixed 61 of them. It is a reminder that "random" in a visualiser is a contract that has to be tested, not assumed.

Colour

The colours come from the live palette: 180 colour schemes, interpolated in OKLab, Björn Ottosson's perceptual colour space from 2020 (Chapter 16 explains why that matters), and the C key changes it. Palette cycling is part of the app's inheritance. JellyDazzle traces its heritage to DAZZLE.EXE, from 1990, and is an homage to James R. Shiflett's DOS program, from the era when changing a few palette registers could make a whole screen of colour move for free.