Mandelbrot experiments

Disclosure: The compute used to produce the renders in this post was sponsored in full by Daytona all compute and equipment used to produce renders belongs to Daytona; go check them out!

A 30s Mandelbrot zoom at 1080p and 60 frames per second contains 1,800 frames. Rendering every frame with four samples per pixel requires 14,929,920,000 sample calculations to produce, even before any calibration or retries. Each sample can require many iterations of the same equation:

The rule fits on one line. Rendering it at extreme magnifications takes rather more work.

At first, an ordinary complex number is enough to locate each sample. As the view shrinks, adjacent coordinates become indistinguishable in floating-point arithmetic. Near the boundary, some orbits take much longer to escape than their neighbours. A movie adds another requirement: the renderer has to preserve its progress when a worker disappears halfway through a run.

This project develops the same calculation through an interactive explorer, a progressive deep renderer, and an adaptive version that distributes frames across GPUs. Following those stages explains both the image and the machinery needed to produce it.

Four rendering runs

Each clip is 30 seconds long, rendered at 1920 × 1080 and 60 frames per second zooming in on 0 Real, 1 Imaginary on a Re, Im coordinate system (the normal mandelbrot set axis) representing a Misiurewicz point.

NVIDIA T4 (~100 USD)

NVIDIA A100 (~5,000 USD)

NVIDIA H200 (~50,000 USD)

20-GPU NVIDIA H100 cluster (~850,000 USD)

The Mandelbrot equation

A complex number has two components:

The real component gives horizontal position; the imaginary component gives vertical position. Each image sample chooses a different value of c, starts its orbit at zero, and repeatedly applies the update.

Writing $z=x+yi$ turns complex squaring into two real equations:

That is the arithmetic underneath the GPU kernel. The same behaviour is easy to inspect without drawing anything:

def orbit(c, updates=8):
    z = 0j
    values = [z]

    for _ in range(updates):
        z = z * z + c
        values.append(z)

    return values


print(orbit(0))
print(orbit(-1))
print(orbit(1))
print(orbit(1j))

#include <cuda_runtime.h>
#include <cstdio>
#include <cstdlib>

void checkCuda(cudaError_t error) {
    if (error != cudaSuccess) {
        std::fprintf(stderr, "%s\n", cudaGetErrorString(error));
        std::exit(EXIT_FAILURE);
    }
}

// Each GPU thread computes one complex orbit.
__global__ void orbit(
    const double2* inputs,
    double2* values,
    int count,
    int updates
) {
    const int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= count) return;

    const double cr = inputs[index].x;
    const double ci = inputs[index].y;

    double zr = 0.0;
    double zi = 0.0;

    const int offset = index * (updates + 1);
    values[offset] = make_double2(zr, zi);

    for (int step = 0; step < updates; ++step) {
        const double next_r = zr * zr - zi * zi + cr;
        const double next_i = 2.0 * zr * zi + ci;

        zr = next_r;
        zi = next_i;

        values[offset + step + 1] = make_double2(zr, zi);
    }
}

int main() {
    constexpr int count = 4;
    constexpr int updates = 8;
    constexpr int stride = updates + 1;

    // Unified memory is accessible from both CPU and GPU.
    double2* inputs = nullptr;
    double2* values = nullptr;

    checkCuda(cudaMallocManaged(
        &inputs, count * sizeof(double2)
    ));
    checkCuda(cudaMallocManaged(
        &values, count * stride * sizeof(double2)
    ));

    inputs[0] = make_double2( 0.0, 0.0);
    inputs[1] = make_double2(-1.0, 0.0);
    inputs[2] = make_double2( 1.0, 0.0);
    inputs[3] = make_double2( 0.0, 1.0);

    // Launch one block; four threads compute the four orbits.
    orbit<<<1, 32>>>(inputs, values, count, updates);
    checkCuda(cudaGetLastError());
    checkCuda(cudaDeviceSynchronize());

    // Print the results on the CPU after the GPU finishes.
    for (int index = 0; index < count; ++index) {
        std::printf("[");

        for (int step = 0; step <= updates; ++step) {
            const double2 z = values[index * stride + step];

            std::printf(
                "%s(%.17g%+.17gj)",
                step == 0 ? "" : ", ",
                z.x,
                z.y
            );
        }

        std::printf("]\n");
    }

    checkCuda(cudaFree(inputs));
    checkCuda(cudaFree(values));

    return 0;
}

import numpy as np
import ray

CUDA_SOURCE = r"""
extern "C" __global__
void orbit(
    const double* inputs,
    double* values,
    int count,
    int updates
) {
    const int index = blockIdx.x * blockDim.x + threadIdx.x;
    if (index >= count) return;

    // Complex numbers are stored as adjacent real/imaginary pairs.
    const double cr = inputs[2 * index];
    const double ci = inputs[2 * index + 1];

    double zr = 0.0;
    double zi = 0.0;

    const int offset = index * (updates + 1) * 2;
    values[offset] = zr;
    values[offset + 1] = zi;

    for (int step = 1; step <= updates; ++step) {
        const double next_r = zr * zr - zi * zi + cr;
        const double next_i = 2.0 * zr * zi + ci;

        zr = next_r;
        zi = next_i;

        values[offset + 2 * step] = zr;
        values[offset + 2 * step + 1] = zi;
    }
}
"""


# Each persistent worker reserves one CPU core and one GPU.
@ray.remote(num_cpus=1, num_gpus=1)
class OrbitWorker:
    def __init__(self):
        import cupy as cp

        # Ray exposes the assigned GPU as device 0 in this process.
        cp.cuda.Device(0).use()
        self.kernel = cp.RawKernel(CUDA_SOURCE, "orbit")

    def compute(self, points, updates=8):
        import cupy as cp

        count = len(points)
        if count == 0:
            return np.empty((0, updates + 1), dtype=np.complex128)

        # Transfer this worker's batch from CPU memory to its GPU.
        inputs = cp.asarray(
            np.column_stack((points.real, points.imag)),
            dtype=cp.float64,
        )
        output = cp.empty(
            (count, updates + 1, 2),
            dtype=cp.float64,
        )

        threads = 256
        blocks = (count + threads - 1) // threads

        self.kernel(
            (blocks,),
            (threads,),
            (
                inputs,
                output,
                np.int32(count),
                np.int32(updates),
            ),
        )

        # This transfer waits for the GPU results to be ready.
        result = cp.asnumpy(output)
        return result[..., 0] + 1j * result[..., 1]


if __name__ == "__main__":
    # Connect to an existing cluster.
    ray.init(address="auto")

    worker_count = 20
    if ray.cluster_resources().get("GPU", 0) < worker_count:
        ray.shutdown()
        raise RuntimeError("This example requires a 20-GPU Ray cluster.")

    # Include the original four examples, followed by a grid of points.
    # The larger grid gives all 20 workers independent work.
    examples = np.array([0, -1, 1, 1j], dtype=np.complex128)
    real = np.linspace(-2.0, 1.0, 640)
    imag = np.linspace(-1.5, 1.5, 360)
    grid = real[None, :] + 1j * imag[:, None]
    points = np.concatenate((examples, grid.ravel()))

    workers = [OrbitWorker.remote() for _ in range(worker_count)]

    try:
        # Assign a separate batch of complete orbits to each worker.
        batches = np.array_split(points, worker_count)

        pending = [
            worker.compute.remote(batch, updates=8)
            for worker, batch in zip(workers, batches)
        ]

        # Workers execute concurrently. Gathering in submission order
        # restores the original ordering of the input points.
        results = np.concatenate(ray.get(pending), axis=0)

        for values in results[:4]:
            print(values.tolist())

    finally:
        # Release the workers and their GPU reservations.
        for worker in workers:
            ray.kill(worker)
        ray.shutdown()

At zero, the sequence stays at zero. At minus one, it alternates between zero and minus one. At one, it grows through 1, 2, 5, and 26, then rapidly beyond them (this can be shown through a code dry run above).

The Mandelbrot set contains the parameters whose orbits remain bounded:

For these orbits starting at zero, reaching a magnitude greater than two establishes escape. The renderer uses a larger bailout radius of 100,000, continuing further because those later values help with the color smoothing and distance estimation. It tests the squared magnitude to avoid taking a square root at every step.

There is an unavoidable difference between this definition and a finished image. Boundedness concerns an infinite sequence; a program stops after a finite number of iterations or a condition is met (in the case of our production code it was 3600 GPU seconds, and the cluster case it was 3600 CPU seconds from the allocator / scheduler). A black sample may belong to the set, or it may be an exterior point that has not escaped yet. Raising the iteration limit helps distinguish those cases.

The deep kernels can also recognise the main cardioid and period-two bulb algebraically in the wide view. Those tests save work in known bounded regions. At deep scales, the code disables the direct floating point tests so rounded coordinates do not decide membership.

The familiar boundary comes from nearby parameters following similar orbits for a long time before their behaviour separates. Its shape belongs to the recurrence. The palette and apparent surface texture are additional choices.

When neighbouring pixels become the same number

The deep renderer measures depth in binary doublings. At depth d, it uses

where V is the horizontal field of view. At 1,000 doublings, the magnification is approximately $10^{301.030}$. The program can display that number through its exponent:

The counter records how much the viewport has shrunk. It says nothing by itself about the accuracy of the rendered pixels.

A float64 has 53 bits of significand precision, roughly 15–16 decimal digits although even with this it can't represent the entire discrete real space of those ~16 decimal digits, this is due to a fundemental limition of the summations of 2s with negative exponentials. Around one, adjacent representable values are separated by about $2^{-52}$. In ordinary Python, 1.0 + 1e-30 == 1.0 evaluates to True at very high precision levels the ability to not actually render floating point numbers exactly becomes a problem. Now there are a variety of ways to fix this issue, we do not have to use the traditional float64 representation for a decimal, we can use different representation forms that use vastly more memory but allow for exact representation (not approximation), furthermore if my memory serves me correctly we can use a dynamically allocated consistency buffer which gives a fixed memory size increase per item rather than the variable of exact replication; to be noted these are both not features of the standard python library to my knowledge and exist as third party work arounds and in all production cases I would not recommend python for this 😂.

That becomes a rendering problem when two pixel coordinates differ by less than rounding can preserve and we are using approximation representation for decimal numbers. The intended points are distinct, but the computer feeds the same parameter into both orbit calculations. More iterations cannot recover a coordinate difference that was discarded before iteration began. The worst case of this would be a infinite cycle where we calculate the exact same image for infinity. Precision and range impose separate limits. Precision determines how finely nearby values can be distinguished. Range determines how large or small a value can be represented at all. Goldberg's discussion of floating-point arithmetic explains the underlying rounding behaviour.

The deep renderer represents a complex quantity using float64 components and a separate integer exponent:

For example, (0.5, 0.0, -4999) represents $2^{-5000}$ without constructing that value as an ordinary float. Multiplication adds exponents; normalization keeps the components in a manageable range. This extends the range while retaining finite significand precision. The other necessary step is to preserve a tiny coordinate offset separately from the much larger centre. Adding them too early would recreate the original rounding problem.

Following a nearby orbit

Perturbation provides that separation.

Choose a reference parameter C with orbit Z, and write a nearby parameter and its orbit as

The reference follows the same rule:

Substituting the displaced orbit into the original equation and subtracting the reference update gives

No term has been discarded. This recurrence is an exact algebraic rearrangement. Numerical error enters when its terms are evaluated with finite arithmetic, and additional approximation enters when the renderer enables iteration skips. Instead of repeatedly forming the full parameter, each sample carries its small offset from the reference. The approach is described in mathr's perturbation notes.

The default destination makes the reference especially cheap. At $C=i$, its orbit is

Every component is a small integer, so the renderer can reproduce the reference values exactly with a short lookup. After the initial steps, only two values repeat. That cycle is repelling. The derivative of a complete trip around it is

Its magnitude describes local amplification, and its argument gives a rotation of $\pi/4$. These properties help explain the nearby structure. The full parameter-dependent calculation still includes nonlinear terms, so successive frames need not be exact rotated copies. The renderer also rebases when the actual orbit value becomes smaller than its displacement from the current reference value. It stores the actual value as the new displacement and restarts the reference index at zero. The sample's orbit and iteration count continue; only their representation changes. This reduces avoidable cancellation. Exact references are a useful special case. The custom-destination path computes a reference using mpmath, then supplies float64 arrays to the GPU. It imposes additional depth limits, and those GPU reference values retain finite precision. The capabilities of the exact periodic destination should not be assumed for every arbitrary coordinate.

Skipping predictable iterations

Keeping tiny offsets representable solves one problem. Calculating all their intermediate steps can still be expensive. When the displacement is sufficiently small, its square contributes very little:

Several linearised steps can be combined into a block:

At the default reference, starting from $Z=-1+i$, two such steps give

Repeating that two-step map eight times produces

The renderer can therefore advance sixteen iterations with

Those integer coefficients are exact for the linearised map. The omission of the quadratic terms remains an approximation.

The adaptive implementation can use larger blocks, up to a configured maximum of 1,024 iterations. It checks the reference phase, remaining iteration count, and displacement exponents before taking a shortcut. The exponent guard keeps the predicted displacement small throughout the selected block. When the guard fails, explicit perturbation steps resume. Position is only part of the state. A skip also has to advance the parameter derivative and iteration count, and account for the stripe contributions that explicit steps would have accumulated. Otherwise an optimisation could preserve escape classification while changing the image's lighting or texture. This is why comparisons between skipped and unskipped calculations need to inspect more than the final black-or-coloured decision.

Making colour and apparent depth

An integer escape count produces abrupt bands. To smooth them, the renderer estimates a fractional count from the orbit's growth beyond the bailout radius.

If escape is detected after N updates, it uses

The motivation is that squaring dominates far outside the set, making the logarithm of the magnitude approximately double at each update. The extra logarithm estimates a position between integer iterations.

A cyclic sinusoidal palette turns the smoothed count into colour. Angular stripe signals add texture by averaging information from the orbit's history. Changing those settings changes the appearance of the same underlying calculation.

The apparent lighting comes from another quantity carried alongside the orbit: its derivative with respect to the pixel parameter.

This derivative measures how sensitively the orbit responds to a small parameter change. The implementation uses it in an exterior distance estimate:

The factor of one half is the convention used here. This is a local estimate for shading, not an exact shortest-distance calculation.

The direction of $z/p$ supplies a direction for the lighting model. Combined with ambient, diffuse, and specular contributions, it produces an embossed appearance without constructing a three-dimensional mesh.

At deep scales, the derivative can become enormous while the estimated distance becomes tiny. The renderer keeps the derivative in extended-exponent form and computes normalized distances through logarithms. That lets the shader use relative scale without materialising values that would overflow or underflow.

A pixel also covers an area, so one sample can miss a narrow branch. The default deep-rendering configuration takes four subpixel samples and approximates averaging in linear light:

The exponent 2.2 is a practical gamma approximation. Increasing the sample count reduces some aliasing and flicker, at the cost of more orbit calculations.

Where the GPU helps

For each sample pass, a CUDA thread handles a pixel. It constructs the sample coordinate, follows the orbit, shades the result, and updates that pixel's accumulated colour. Across GPUs, workers receive complete frames to render independently.

The work inside one orbit remains sequential: its next value depends on its current value. Parallelism comes from the number of samples and frames available.

Some samples escape quickly while others remain active near the iteration limit. When neighbouring threads take different paths or execute different numbers of iterations, part of a GPU warp can sit idle while the remaining threads continue. The image's mathematical structure therefore affects execution efficiency.

The implementation uses ordinary FP64 orbit arithmetic. It does not use Tensor Core matrix operations, and adding GPU memory alone does not resolve the arithmetic and control-flow costs.

GpuRenderer reuses its device buffers across frames. It accumulates colours in float32, resolves them to RGB bytes on the GPU, and downloads that smaller image. At 1080p, the float32 colour array occupies about 24.9 MB; the RGB byte array occupies about 6.2 MB.

A thirty-second movie makes those repeated decisions consequential. Every extra sample, transfer, or explicit orbit step is multiplied across the output frames.

Spending compute on a camera path

The progressive renderer advances by a fixed zoom increment until its measured GPU-time budget is used. The completed source frames determine the deepest view available, and the exporter fits that journey into the requested playback duration.

The adaptive version begins with the final frame count and adjusts future camera depths to fit the remaining compute. It measures calibration probes at the final resolution and sample settings, then estimates rendering time at depth d with

This is an empirical scheduling model. Escape behaviour and iteration skips can change the cost, so the measurements continue to inform future assignments.

The scheduler subtracts completed work and the predicted cost of jobs already assigned before planning additional frames. This reservation prevents it from allocating the same remaining budget to several workers.

Once a job receives its camera positions, those positions stay fixed. A retry renders the same planned frames even if later measurements change the estimated affordable depth. Adaptation affects unassigned frames.

Visual continuity places another limit on the camera. If the magnification ratio between adjacent frames is capped at $r_z$, the depth step must satisfy

With 1,800 frames, there are only 1,799 intervals over which to spread the zoom. A deeper endpoint eventually requires larger jumps or more frames. Additional compute cannot remove that relationship.

The budget's unit also matters. CUDA events measure the sample and resolve kernels. Image downloads, encoding, setup, and uploads fall outside that interval. Allocated machines can continue costing money while doing those things, so measured GPU time, wall time, and billed resource time need separate accounting.

Keeping completed work after a failure

Rendering frames independently makes recovery possible, provided the system remembers what each worker was supposed to produce.

The adaptive queue saves an immutable plan before assigning it. A worker receives a temporary lease, renews that lease through heartbeats, renders a short clip, and publishes the clip and its receipt to shared storage. The receipt records the camera metadata, timing, and artifact information needed to accept the result.

If the worker disappears, its lease can expire and another worker can receive the same plan. The coordinator rejects stale completions and handles repeated commits without counting an accepted job twice. Physical rendering may happen more than once; accepted accounting is deduplicated.

That distinction matters when a worker finishes its GPU work but loses the connection before receiving confirmation. Durable plans and receipts let the coordinator reconstruct progress after a restart without depending entirely on its local database.

Encoding runs in the background while the GPU proceeds to another batch. A bounded queue limits uncompressed frames waiting for the encoder. When that queue fills, rendering waits instead of consuming memory indefinitely.

The intermediate clips use H.264, so they preserve usable movie frames with compression loss. During final assembly, the pipeline checks clip hashes and validates frame order and timing. It can concatenate encoded clips directly, or rebuild timestamps through re-encoding when needed. If an incomplete render requires a repeated final frame, the output records that padding.

A completed MP4 is therefore the result of both numerical work and bookkeeping. A correct camera plan is of little use if its frames are missing or arrive in the wrong order.

Checking what the image claims

A plausible fractal image can conceal a numerical error. A finite colour value can still be wrong, and two implementations can agree because they share the same mistake.

The project's numerical audits use a direct high-precision calculation of z = z*z + c as an independent comparison. That calculation avoids perturbation, rebasing, and iteration skips. Selected samples are compared on escape classification, smoothed count, distance estimates, and lighting direction. Skipped and unskipped paths are also compared for stripe accumulation.

The project documentation records direct comparisons through 16,384 binary doublings. Its adaptive CPU audit reports 162 sampled cases with matching escape classifications and a maximum smooth-count error of approximately $1.88\times10^{-9}$. Separate skipped-versus-unskipped checks reach 100,000 doublings.

Those are recorded results for selected cases, with different checks establishing different things. The 100,000-doubling comparison is not a direct high-precision certification of every pixel at that depth. A scheduler simulation does not establish physical multi-GPU speedup, and a matching file checksum establishes byte integrity rather than mathematical accuracy.

The zoom counter remains useful: it reports the scale of the requested view. Confidence in the picture comes from comparing calculations and checking convergence as iteration limits or sampling settings increase.

At an extreme magnification, the equation has not changed. What has changed is the care needed to preserve the distinction between neighbouring parameters. The renderer has to carry that distinction through its arithmetic and into the finished frames. Otherwise it can produce a smooth, convincing movie of an error.