Mandelbrot

One rule, z → z² + c, repeated from zero. Points whose sequence stays small are black; the colours show how fast the others escape. Drag to move, scroll or pinch to zoom, and tap a point to see its orbit. How it works.

Una regla, z → z² + c, repetida desde cero. Los puntos cuya secuencia se mantiene pequeña son negros; los colores muestran lo rápido que escapan los demás. Arrastra para moverte, usa la rueda o pellizca para hacer zoom, y toca un punto para ver su órbita. Cómo funciona.

The idea in plain words

The whole picture comes from one rule: take a number, square it, add c. Start at zero and repeat. Each point c of the plane gives a different sequence. For some points the numbers stay small forever. For others they grow without limit. Paint the first kind black and you have the Mandelbrot set. That is the entire definition.

The numbers are complex numbers, which are simply points on a plane. Squaring one has a picture: it doubles the point's angle and squares its distance from the centre. So points near the centre get pulled in, points far away get flung out, and adding c nudges them every step. The border between "pulled in" and "flung out" is where all the detail lives.

Tap any point of the picture to see its journey. The white path is the sequence z₀ = 0, z₁, z₂… Inside the set it settles into a loop or a spiral. Outside it crosses the dashed circle and runs away. Near the border it wanders for a long time before deciding.

The colours outside the set are not part of the set. They show how fast each point escapes, which makes the invisible fine structure near the border visible. Zoom into the edge and the detail never ends; tiny copies of the whole set keep appearing.

Why it matters

It is a clean example of emergence: a rule you can write on a napkin produces structure no one designed. The same rule restricted to real numbers is the logistic map, a simple model of population growth. The famous cascade of period doublings that leads it into chaos appears here as the chain of bulbs along the horizontal axis.

Robert Brooks and Peter Matelski drew the first rough picture of it in 1978. Benoit Mandelbrot, working at IBM, made it famous from 1980 on with computer graphics. Adrien Douady and John Hubbard studied it in depth and gave it Mandelbrot's name.

Is there a theorem?

The circle of no return. If the sequence ever gets farther than 2 from the centre, it escapes for sure. This is easy to prove (see step 2 below), and it is what makes the picture computable at all.
It is connected (Douady and Hubbard, 1982). The little islands you see when you zoom in are all joined to the main body by filaments too thin to draw.
Its border is as wiggly as a border can be (Shishikura, 1998). It has Hausdorff dimension 2: a curve so crinkled that, in this sense, it fills area.
Still open: whether it is "locally connected" (the MLC conjecture). It sounds technical, but it is one of the central open problems in complex dynamics.

The equations

The rule: z_n+1 = z_n² + c, with z₀ = 0
The set: M = { c : |z_n| ≤ 2 for every n }
Squaring, with z = x + yi: z² = (x² − y²) + 2xy·i
Smooth escape count: ν = n + 1 − log₂(log |z_n|)

What the picture does and does not say

Black means "did not escape within the iteration limit". Some black points near the border are really outside; raise Iterations and they light up.
The colours are a choice. Only black and not-black belong to the mathematics.
The GPU uses 32-bit floats, so a direct zoom stops at about ×100,000. Perturbation (step 9) takes it to about ×10¹³.

How it works

This part builds the explorer step by step. It covers:

Complex numbers as plain pairs of numbers.
The escape-time algorithm, and why the number 2.
Why a pixel-by-pixel loop is too slow, and how the GPU fixes it.
Three ways to colour: inside/outside, bands and smooth.
Zooming around the cursor, panning and pinching.
Drawing the orbit of a point.
Measuring and fixing performance on phones.
Zooming ten billion times deeper with perturbation.

Plain JavaScript and raw WebGL, no libraries. The scene is about 285 lines; the plumbing it shares with the Julia scene (sizing, navigation, adaptive resolution, GPU reset) is a module of about 300.

1. Complex numbers are points

A complex number x + yi is two ordinary numbers. We never need a special type: we keep x and y as two variables. The only rule to remember is i² = −1, which gives the formula for squaring.

// A complex number z = x + yi is just a point (x, y).
// z² = (x + yi)² = x² − y² + 2xy·i, because i² = −1.
function square(x, y) {
  return [x * x - y * y, 2 * x * y];
}

2. Escape time

Start at z = 0, apply z = z² + c again and again, and count the steps until |z| passes 2. We compare x² + y² with 4 to avoid a square root. If it never escapes within maxIter steps, we call the point inside.

function escapeTime(cx, cy, maxIter) {
  let x = 0, y = 0;                       // z₀ = 0
  for (let n = 0; n < maxIter; n++) {
    [x, y] = [x * x - y * y + cx, 2 * x * y + cy];   // z = z² + c
    if (x * x + y * y > 4) return n;      // |z| > 2: it escapes
  }
  return maxIter;                         // never escaped: inside (as far as we know)
}

Switch the scene to Inside / outside: this function, and nothing else, draws that picture.

Why exactly 2?

Suppose |z| > 2 and |z| ≥ |c|. Then |z² + c| ≥ |z|² − |c| ≥ |z|² − |z| = |z|(|z| − 1). Since |z| − 1 > 1, the next value is farther out than this one, by a factor that keeps growing. Points of the set all have |c| ≤ 2, so once an orbit passes 2 it is gone for good.

3. Painting the plane, slowly

To draw the picture, run escapeTime once per pixel. Each pixel maps to a point c of the plane: the centre of the canvas is the centre of the view, and the canvas height covers height units.

for (let py = 0; py < H; py++) {
  for (let px = 0; px < W; px++) {
    const cx = centerX + ((px - W / 2) / H) * height;
    const cy = centerY - ((py - H / 2) / H) * height;
    const n = escapeTime(cx, cy, maxIter);
    setPixel(px, py, n === maxIter ? INSIDE : colorFor(n));
  }
}

This works, but a 1000 × 600 canvas at 300 iterations is up to 180 million steps per frame. Panning becomes a slideshow.

4. Moving to the GPU

Every pixel is independent of the others, which is exactly the kind of work a graphics card is built for. A fragment shader is a small program the GPU runs once for every pixel, thousands at a time. We move the same loop into it, written in GLSL.

precision highp float;
uniform vec2 uRes;       // canvas size in pixels
uniform vec2 uCenter;    // centre of the view in the plane
uniform float uHeight;   // height of the view in the plane
uniform float uMaxIter;

void main() {
  vec2 p = (gl_FragCoord.xy - 0.5 * uRes) / uRes.y;
  vec2 c = uCenter + p * uHeight;
  vec2 z = vec2(0.0);
  float n = 0.0;
  bool escaped = false;
  for (int i = 0; i < 5000; i++) {
    if (float(i) >= uMaxIter) break;
    z = vec2(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c;
    if (dot(z, z) > 4.0) { escaped = true; break; }
    n += 1.0;
  }
  …
}

To make the shader run on every pixel, we draw two triangles that cover the whole canvas. The shader receives each pixel's position in gl_FragCoord and turns it into c, exactly as the CPU loop did.

// Two triangles that cover the screen; the fragment shader does the rest.
gl.bufferData(gl.ARRAY_BUFFER,
  new Float32Array([-1, -1, 1, -1, -1, 1,  -1, 1, 1, -1, 1, 1]),
  gl.STATIC_DRAW);
…
gl.drawArrays(gl.TRIANGLES, 0, 6);

Why the strange loop with 5000 and a break?

In WebGL 1, a loop must have a fixed limit known when the shader is compiled. So we loop up to a constant and break out when we reach the iteration count chosen on the slider.

5. Colour

Painting each escape count with its own colour gives bands: you can see each step as a stripe. The palette is a cosine wave per colour channel, a trick by Inigo Quilez that gives smooth, cyclic gradients with a single line.

For the smooth mode we want a fractional step count. Far from the set, each step roughly squares |z|, so log(log|z|) grows by log 2 per step. Measuring how far along that last step we got gives the fraction. It only works when |z| is already large, so in smooth mode we let points run until |z| > 256 before stopping. That changes the colours, not the set.

vec3 palette(float t) {
  // Cosine palette (Inigo Quilez): a + b·cos(2π(c·t + d))
  return vec3(0.5) + vec3(0.5) * cos(6.28318 * (vec3(1.0) * t + vec3(0.62, 0.70, 0.80)));
}
…
float bail = uMode == 2 ? 256.0 : 2.0;
…
if (uMode == 2) n = n + 1.0 - log2(log(length(z)));   // smooth
float t = uMode == 1 ? floor(n) : n;                   // bands keep whole steps
gl_FragColor = vec4(palette(0.02 * t), 1.0);

6. Zoom and pan

Zooming should keep the point under your cursor or fingers fixed. The trick is to remember which point of the plane was under the cursor, change the scale, and then move the view so that point is back under the cursor.

// Zoom by `factor` keeping the plane point under (px, py) fixed on screen.
view.zoomAt = (px, py, factor) => {
  const before = view.toPlane(px, py);
  view.h = Math.min(maxH, Math.max(minH(), view.h * factor));
  const after = view.toPlane(px, py);
  view.cx += before.x - after.x;
  view.cy += before.y - after.y;
  view.moved();
};

The mouse wheel and the +/− buttons call zoomAt (the buttons zoom around the centre). Pinching uses the same idea with two fingers: the point that was between them stays under their midpoint, so you can zoom and move at once. Whenever a finger is added or lifted, the gesture starts again from where the view is, or the view would jump. A drag moves the centre; a tap that does not move picks the point for the orbit.

Why does zoom stop at about ×100,000?

GPU floats have about 7 significant digits. When the view is 0.00002 units tall and sits near −0.75, neighbouring pixels differ in the 8th digit, so they get the same c and the picture turns into blocks. Deep-zoom explorers use double precision, or "perturbation": one precise orbit on the CPU and small differences on the GPU. Step 9 does exactly that.

7. Orbits

The orbit is the same iteration, run on the CPU for one point, keeping every step. We draw it on a transparent 2D canvas stacked on top of the WebGL one, so it can change without redrawing the fractal.

function orbitOf(cx, cy, maxIter) {
  const pts = [[0, 0]];
  let x = 0, y = 0;
  for (let i = 0; i < maxIter; i++) {
    [x, y] = [x * x - y * y + cx, 2 * x * y + cy];   // z = z² + c
    pts.push([x, y]);
    if (x * x + y * y > 4) return { pts, escaped: true };
  }
  return { pts, escaped: false };
}

Try points in the main cardioid (they spiral to one fixed point), in the big left bulb (they alternate between two), and just outside the border (they wander, then leave).

8. Fast enough for a phone

The first version froze on phones after some zooming. Three things add up. Near the border almost every pixel needs many iterations. Phone screens have two or three device pixels per CSS pixel, so there are four to nine times more pixels to compute. And if a single frame takes seconds, the browser decides the GPU has hung and resets it: the WebGL "context" is lost and the picture stops updating for good.

Before optimising, measure. Tick Show diagnostics: it shows how long the last frame took on the GPU, at what resolution, and whether this GPU supports high-precision floats.

const t0 = performance.now();
…
gl.drawArrays(gl.TRIANGLES, 0, 6);
// Reading one pixel waits for the GPU to finish, so the timing includes the real work.
gl.readPixels(0, 0, 1, 1, gl.RGBA, gl.UNSIGNED_BYTE, new Uint8Array(4));
view.frameMs = performance.now() - t0;
…
if (!view.moving && view.frameMs > SLOW_MS && idleScale > 1) {
  idleScale = Math.max(1, idleScale * 0.7);
  setScale(idleScale);
  view.requestDraw();
}

Skip the easy black. The big black shapes are the slowest pixels, because they run every iteration and never escape. But two of them have exact formulas: the main cardioid and the circle of the left bulb. A pixel inside either one is black, with no loop at all.

// Points in the main cardioid or the big left bulb never escape: skip the loop.
vec2 d = c - vec2(0.25, 0.0);
float q = dot(d, d);
vec2 b = c + vec2(1.0, 0.0);
if (q * (q + d.x) <= 0.25 * c.y * c.y || dot(b, b) <= 0.0625) {
  gl_FragColor = vec4(inside, 1.0);
  return;
}

Fewer pixels while moving. While you drag or pinch, the eye cannot see fine detail anyway. So we start at half a pixel per CSS pixel (four times fewer pixels than a desktop at full size, and up to sixteen times fewer than a phone) and let the browser stretch the image. Then the moving resolution adapts to the GPU: sharper if moving frames are cheap, coarser if not (step 8 of the Julia scene). When the view has been still for 180 ms, we render once more at full resolution. If even that resting frame takes longer than 250 ms, the resting resolution goes down a notch.

const MOVE_SCALE = 0.5;   // starting render pixels per CSS pixel while moving

function setScale(s) {
  view.scale = s;
  canvas.width = Math.max(1, Math.round(view.W * s));
  canvas.height = Math.max(1, Math.round(view.H * s));
}

// Called on every pan, zoom or parameter drag: drop to low resolution, and go back to full
// detail once nothing has moved for a moment.
view.moved = () => {
  if (!view.moving) { view.moving = true; lastDraw = 0; }
  if (view.scale !== moveScale) setScale(moveScale);
  clearTimeout(settleTimer);
  settleTimer = setTimeout(() => { view.moving = false; setScale(idleScale); view.requestDraw(); }, SETTLE_MS);
  view.requestDraw();
};

Survive a reset. If the browser still resets the GPU, we ask for it back, rebuild the shader and continue at a lower resolution instead of freezing.

canvas.addEventListener("webglcontextlost", (e) => {
  e.preventDefault();   // "I want it back"
  view.lost = true;
  paintOverlay();
});
canvas.addEventListener("webglcontextrestored", () => {
  view.lost = false;
  idleScale = 1;        // come back at a lower resolution
  initGL();
  setScale(idleScale);
  view.requestDraw();
});

Why not just lower the iteration count?

It helps, but it changes the picture: with too few iterations the border fills with false black. Lowering resolution only changes sharpness, and only while you move.

9. Deep zoom with perturbation

The precision wall from step 6 is about the numbers, not the GPU's speed. A float32 has about 7 significant digits. At ×10¹⁰, two neighbouring pixels differ in the 12th digit, so every pixel gets the same c. Untick Deep zoom at the "Dive ×10¹⁰" place and the picture turns into a single flat colour.

The fix, found by K. I. Martin in 2013, is to stop asking every pixel for its full position. Think of a satellite photo. Measuring the exact latitude of every house to ten decimals is hard. Instead, measure one landmark very precisely, and for each house measure only "12 metres east of the landmark". The small offsets are easy, and the precise part is done once.

Here the landmark is the reference orbit: the centre of the view, iterated on the CPU in float64 (16 digits), which is cheap because it is a single point. Every pixel then tracks only its difference from that orbit. Those differences are tiny, and tiny numbers are exactly what floats store well: a float32 can hold 10⁻³⁰ with its 7 digits intact.

// Reference orbit: one point (the centre of the view) iterated in float64 on the CPU.
function referenceOrbit(cx, cy, maxIter) {
  const zs = [0, 0];
  let x = 0, y = 0;
  for (let i = 0; i < maxIter; i++) {
    [x, y] = [x * x - y * y + cx, 2 * x * y + cy];
    if (x * x + y * y > 65536) break;
    zs.push(x, y);
  }
  return zs;   // [X0, Y0, X1, Y1, …], uploaded to the GPU as a float texture
}

The algebra is short. Write the pixel's orbit as the reference plus a difference, z = Z + δ, and its point as c = C + Δc. Put that into z² + c and subtract the reference step Z² + C. What is left is the rule for the difference: δ' = 2·Z·δ + δ² + Δc. It never needs the big numbers, only Z, which the texture provides.

// z = Z + dz, c = C + dc. From z² + c: dz' = 2·Z·dz + dz² + dc. Only small numbers here.
vec2 dc = p * uHeight;
vec2 dz = vec2(0.0);
float m = 0.0;
for (int i = 0; i < 5000; i++) {
  if (float(i) >= uMaxIter) break;
  dz = 2.0 * cmul(refAt(m), dz) + cmul(dz, dz) + dc;
  m += 1.0;
  z = refAt(m) + dz;
  if (dot(z, z) > bail * bail) { escaped = true; break; }
  n += 1.0;
  // Rebasing: if the pixel is nearer to 0 than to the reference, or the reference has run
  // out, restart the reference from Z0 = 0 and carry the whole value as the difference.
  if (dot(z, z) < dot(dz, dz) || m >= uRefLen - 1.0) { dz = z; m = 0.0; }
}

What are "glitches", and what does rebasing do?

The trick assumes each pixel stays close to the reference. Sometimes a pixel's orbit passes near 0 while the reference is far away, and the difference becomes larger than the value it describes. Precision is lost and early perturbation renderers showed flat blobs called glitches. Since 2021, renderers use rebasing: when the pixel is nearer to 0 than to the reference, or when the reference orbit has escaped and ended, restart the reference from its beginning (Z₀ = 0) and carry the whole value as the new difference. It costs one comparison.

Why does it stop at about ×10¹³?

Now the limit is the reference and the position of the view, both stored in float64 (about 16 digits). Going deeper needs arbitrary-precision numbers for that one orbit, for example JavaScript's BigInt as fixed-point. The GPU side would not change.

Deep places need more iterations, because points near the border take longer to decide. The "Dive ×10¹⁰" preset raises the count to 2000. That is the expensive part on a phone, and the adaptive resolution from step 8 is what keeps it usable.

Going further

Fix c and vary the starting point instead: that draws a Julia set, one for every c. It is the next scene: Julia.
Compute the reference orbit with arbitrary precision to zoom past ×10¹³.
Colour by the orbit's shape (orbit traps) instead of its escape time.

La idea en palabras simples

Toda la imagen sale de una regla: toma un número, elévalo al cuadrado y súmale c. Empieza en cero y repite. Cada punto c del plano da una secuencia distinta. En algunos puntos los números se quedan pequeños para siempre. En otros crecen sin límite. Pinta de negro los primeros y tienes el conjunto de Mandelbrot. Esa es toda la definición.

Los números son números complejos, que no son más que puntos en un plano. Elevar uno al cuadrado tiene un dibujo: duplica el ángulo del punto y eleva al cuadrado su distancia al centro. Así, los puntos cercanos al centro se ven atraídos, los lejanos salen despedidos, y sumar c los empuja un poco en cada paso. La frontera entre "atraído" y "despedido" es donde vive todo el detalle.

Toca cualquier punto de la imagen para ver su recorrido. El camino blanco es la secuencia z₀ = 0, z₁, z₂… Dentro del conjunto acaba en un ciclo o en una espiral. Fuera, cruza el círculo discontinuo y se escapa. Cerca de la frontera vaga mucho tiempo antes de decidirse.

Los colores de fuera no son parte del conjunto. Muestran lo rápido que escapa cada punto, y así hacen visible la estructura fina cerca del borde. Acércate al borde y el detalle no se acaba nunca; siguen apareciendo copias diminutas del conjunto entero.

Por qué importa

Es un ejemplo limpio de emergencia: una regla que cabe en una servilleta produce una estructura que nadie diseñó. La misma regla restringida a los números reales es el mapa logístico, un modelo sencillo de crecimiento de poblaciones. La famosa cascada de duplicaciones de periodo que lo lleva al caos aparece aquí como la cadena de bulbos sobre el eje horizontal.

Robert Brooks y Peter Matelski dibujaron la primera imagen tosca en 1978. Benoit Mandelbrot, que trabajaba en IBM, lo hizo famoso desde 1980 con gráficos por ordenador. Adrien Douady y John Hubbard lo estudiaron a fondo y le dieron el nombre de Mandelbrot.

¿Hay un teorema?

El círculo sin retorno. Si la secuencia se aleja alguna vez más de 2 del centro, se escapa con seguridad. Es fácil de demostrar (ver el paso 2 más abajo) y es lo que permite calcular la imagen.
Es conexo (Douady y Hubbard, 1982). Las islitas que ves al acercarte están todas unidas al cuerpo principal por filamentos demasiado finos para dibujarlos.
Su borde es tan arrugado como puede serlo un borde (Shishikura, 1998). Tiene dimensión de Hausdorff 2: una curva tan plegada que, en ese sentido, llena área.
Aún abierto: si es "localmente conexo" (la conjetura MLC). Suena técnico, pero es uno de los problemas abiertos centrales de la dinámica compleja.

Las ecuaciones

La regla: z_n+1 = z_n² + c, con z₀ = 0
El conjunto: M = { c : |z_n| ≤ 2 para todo n }
Elevar al cuadrado, con z = x + yi: z² = (x² − y²) + 2xy·i
Conteo de escape suave: ν = n + 1 − log₂(log |z_n|)

Qué dice la imagen y qué no

Negro significa "no escapó dentro del límite de iteraciones". Algunos puntos negros cerca del borde en realidad están fuera; sube Iteraciones y se encienden.
Los colores son una elección. A las matemáticas solo les pertenecen el negro y el no negro.
La GPU usa números de 32 bits, así que el zoom directo se detiene en torno a ×100 000. La perturbación (paso 9) lo lleva a unos ×10¹³.

Cómo funciona

Esta parte construye el explorador paso a paso. Trata:

Los números complejos como simples pares de números.
El algoritmo de tiempo de escape, y por qué el número 2.
Por qué un bucle píxel a píxel es demasiado lento, y cómo lo resuelve la GPU.
Tres maneras de colorear: dentro/fuera, bandas y suave.
Hacer zoom alrededor del cursor, desplazar y pellizcar.
Dibujar la órbita de un punto.
Medir y arreglar el rendimiento en móviles.
Acercarse diez mil millones de veces más con perturbación.

JavaScript puro y WebGL a pelo, sin librerías. La escena tiene unas 285 líneas; la fontanería que comparte con la escena del Julia (tamaño, navegación, resolución adaptativa, reinicio de la GPU) es un módulo de unas 300. Los comentarios del código están en inglés, como en el repositorio.

1. Los números complejos son puntos

Un número complejo x + yi son dos números corrientes. No hace falta un tipo especial: guardamos x e y en dos variables. La única regla que hay que recordar es i² = −1, y de ella sale la fórmula del cuadrado.

// A complex number z = x + yi is just a point (x, y).
// z² = (x + yi)² = x² − y² + 2xy·i, because i² = −1.
function square(x, y) {
  return [x * x - y * y, 2 * x * y];
}

2. Tiempo de escape

Empieza en z = 0, aplica z = z² + c una y otra vez y cuenta los pasos hasta que |z| pase de 2. Comparamos x² + y² con 4 para ahorrarnos la raíz cuadrada. Si no escapa en maxIter pasos, decimos que el punto está dentro.

function escapeTime(cx, cy, maxIter) {
  let x = 0, y = 0;                       // z₀ = 0
  for (let n = 0; n < maxIter; n++) {
    [x, y] = [x * x - y * y + cx, 2 * x * y + cy];   // z = z² + c
    if (x * x + y * y > 4) return n;      // |z| > 2: it escapes
  }
  return maxIter;                         // never escaped: inside (as far as we know)
}

Cambia la escena a Dentro / fuera: esta función, y nada más, dibuja esa imagen.

¿Por qué justo 2?

Supón que |z| > 2 y |z| ≥ |c|. Entonces |z² + c| ≥ |z|² − |c| ≥ |z|² − |z| = |z|(|z| − 1). Como |z| − 1 > 1, el siguiente valor queda más lejos que este, por un factor que no deja de crecer. Todos los puntos del conjunto cumplen |c| ≤ 2, así que una órbita que pasa de 2 ya no vuelve.

3. Pintar el plano, despacio

Para dibujar la imagen, se ejecuta escapeTime una vez por píxel. Cada píxel corresponde a un punto c del plano: el centro del canvas es el centro de la vista, y la altura del canvas abarca height unidades.

for (let py = 0; py < H; py++) {
  for (let px = 0; px < W; px++) {
    const cx = centerX + ((px - W / 2) / H) * height;
    const cy = centerY - ((py - H / 2) / H) * height;
    const n = escapeTime(cx, cy, maxIter);
    setPixel(px, py, n === maxIter ? INSIDE : colorFor(n));
  }
}

Funciona, pero un canvas de 1000 × 600 a 300 iteraciones son hasta 180 millones de pasos por fotograma. Desplazarse se vuelve una presentación de diapositivas.

4. Pasar a la GPU

Cada píxel es independiente de los demás, y ese es justo el trabajo para el que está hecha una tarjeta gráfica. Un fragment shader es un programa pequeño que la GPU ejecuta una vez por píxel, miles a la vez. Llevamos a él el mismo bucle, escrito en GLSL.

precision highp float;
uniform vec2 uRes;       // canvas size in pixels
uniform vec2 uCenter;    // centre of the view in the plane
uniform float uHeight;   // height of the view in the plane
uniform float uMaxIter;

void main() {
  vec2 p = (gl_FragCoord.xy - 0.5 * uRes) / uRes.y;
  vec2 c = uCenter + p * uHeight;
  vec2 z = vec2(0.0);
  float n = 0.0;
  bool escaped = false;
  for (int i = 0; i < 5000; i++) {
    if (float(i) >= uMaxIter) break;
    z = vec2(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c;
    if (dot(z, z) > 4.0) { escaped = true; break; }
    n += 1.0;
  }
  …
}

Para que el shader se ejecute en cada píxel, dibujamos dos triángulos que cubren todo el canvas. El shader recibe la posición de cada píxel en gl_FragCoord y la convierte en c, igual que hacía el bucle en la CPU.

// Two triangles that cover the screen; the fragment shader does the rest.
gl.bufferData(gl.ARRAY_BUFFER,
  new Float32Array([-1, -1, 1, -1, -1, 1,  -1, 1, 1, -1, 1, 1]),
  gl.STATIC_DRAW);
…
gl.drawArrays(gl.TRIANGLES, 0, 6);

¿Por qué ese bucle raro con 5000 y un break?

En WebGL 1, un bucle tiene que tener un límite fijo conocido al compilar el shader. Así que el bucle llega hasta una constante y se sale cuando alcanza el número de iteraciones elegido en el slider.

5. Color

Pintar cada número de pasos con su propio color da bandas: cada paso se ve como una franja. La paleta es una onda coseno por canal de color, un truco de Inigo Quilez que da degradados suaves y cíclicos con una sola línea.

Para el modo suave queremos un número de pasos fraccionario. Lejos del conjunto, cada paso eleva |z| más o menos al cuadrado, así que log(log|z|) crece log 2 por paso. Medir cuánto avanzamos en ese último paso da la fracción. Solo funciona cuando |z| ya es grande, así que en el modo suave dejamos correr los puntos hasta |z| > 256 antes de parar. Eso cambia los colores, no el conjunto.

vec3 palette(float t) {
  // Cosine palette (Inigo Quilez): a + b·cos(2π(c·t + d))
  return vec3(0.5) + vec3(0.5) * cos(6.28318 * (vec3(1.0) * t + vec3(0.62, 0.70, 0.80)));
}
…
float bail = uMode == 2 ? 256.0 : 2.0;
…
if (uMode == 2) n = n + 1.0 - log2(log(length(z)));   // smooth
float t = uMode == 1 ? floor(n) : n;                   // bands keep whole steps
gl_FragColor = vec4(palette(0.02 * t), 1.0);

6. Zoom y desplazamiento

El zoom debe dejar quieto el punto que está bajo el cursor o los dedos. El truco es recordar qué punto del plano estaba bajo el cursor, cambiar la escala y mover la vista para que ese punto vuelva a quedar bajo el cursor.

// Zoom by `factor` keeping the plane point under (px, py) fixed on screen.
view.zoomAt = (px, py, factor) => {
  const before = view.toPlane(px, py);
  view.h = Math.min(maxH, Math.max(minH(), view.h * factor));
  const after = view.toPlane(px, py);
  view.cx += before.x - after.x;
  view.cy += before.y - after.y;
  view.moved();
};

La rueda del ratón y los botones +/− llaman a zoomAt (los botones acercan alrededor del centro). El pellizco usa la misma idea con dos dedos: el punto que estaba entre ellos se queda bajo su punto medio, así que puedes acercarte y moverte a la vez. Cada vez que se apoya o se levanta un dedo, el gesto vuelve a empezar desde donde está la vista; si no, la vista saltaría. Arrastrar mueve el centro; un toque sin movimiento elige el punto de la órbita.

¿Por qué el zoom se detiene en torno a ×100 000?

Los números de la GPU tienen unas 7 cifras significativas. Cuando la vista mide 0,00002 unidades de alto y está cerca de −0,75, los píxeles vecinos difieren en la octava cifra, así que reciben el mismo c y la imagen se llena de bloques. Los exploradores de zoom profundo usan doble precisión o "perturbación": una órbita precisa en la CPU y pequeñas diferencias en la GPU. El paso 9 hace justo eso.

7. Órbitas

La órbita es la misma iteración, ejecutada en la CPU para un solo punto y guardando cada paso. La dibujamos en un canvas 2D transparente encima del de WebGL, así puede cambiar sin volver a dibujar el fractal.

function orbitOf(cx, cy, maxIter) {
  const pts = [[0, 0]];
  let x = 0, y = 0;
  for (let i = 0; i < maxIter; i++) {
    [x, y] = [x * x - y * y + cx, 2 * x * y + cy];   // z = z² + c
    pts.push([x, y]);
    if (x * x + y * y > 4) return { pts, escaped: true };
  }
  return { pts, escaped: false };
}

Prueba puntos en el cardioide principal (giran en espiral hacia un punto fijo), en el bulbo grande de la izquierda (alternan entre dos) y justo fuera del borde (vagan y luego se van).

8. Lo bastante rápido para un móvil

La primera versión se congelaba en el móvil después de algo de zoom. Se suman tres cosas. Cerca del borde casi todos los píxeles necesitan muchas iteraciones. Las pantallas de móvil tienen dos o tres píxeles físicos por píxel CSS, así que hay entre cuatro y nueve veces más píxeles que calcular. Y si un solo fotograma tarda segundos, el navegador cree que la GPU se colgó y la reinicia: se pierde el "contexto" de WebGL y la imagen deja de actualizarse para siempre.

Antes de optimizar, hay que medir. Marca Mostrar diagnóstico: muestra cuánto tardó el último fotograma en la GPU, a qué resolución, y si esta GPU admite números de alta precisión.

const t0 = performance.now();
…
gl.drawArrays(gl.TRIANGLES, 0, 6);
// Reading one pixel waits for the GPU to finish, so the timing includes the real work.
gl.readPixels(0, 0, 1, 1, gl.RGBA, gl.UNSIGNED_BYTE, new Uint8Array(4));
view.frameMs = performance.now() - t0;
…
if (!view.moving && view.frameMs > SLOW_MS && idleScale > 1) {
  idleScale = Math.max(1, idleScale * 0.7);
  setScale(idleScale);
  view.requestDraw();
}

Saltarse el negro fácil. Las grandes formas negras son los píxeles más lentos, porque hacen todas las iteraciones y nunca escapan. Pero dos de ellas tienen fórmula exacta: el cardioide principal y el círculo del bulbo izquierdo. Un píxel dentro de cualquiera de los dos es negro, sin bucle.

// Points in the main cardioid or the big left bulb never escape: skip the loop.
vec2 d = c - vec2(0.25, 0.0);
float q = dot(d, d);
vec2 b = c + vec2(1.0, 0.0);
if (q * (q + d.x) <= 0.25 * c.y * c.y || dot(b, b) <= 0.0625) {
  gl_FragColor = vec4(inside, 1.0);
  return;
}

Menos píxeles mientras te mueves. Mientras arrastras o pellizcas, el ojo no ve el detalle fino. Así que empezamos a medio píxel por píxel CSS (cuatro veces menos píxeles que un ordenador a tamaño completo, y hasta dieciséis veces menos que un móvil) y dejamos que el navegador estire la imagen. Después, la resolución en movimiento se adapta a la GPU: más nítida si los fotogramas en movimiento son baratos, más basta si no (paso 8 de la escena del Julia). Cuando la vista lleva 180 ms quieta, se pinta una vez más a resolución completa. Si incluso ese fotograma en reposo tarda más de 250 ms, la resolución en reposo baja un escalón.

const MOVE_SCALE = 0.5;   // starting render pixels per CSS pixel while moving

function setScale(s) {
  view.scale = s;
  canvas.width = Math.max(1, Math.round(view.W * s));
  canvas.height = Math.max(1, Math.round(view.H * s));
}

// Called on every pan, zoom or parameter drag: drop to low resolution, and go back to full
// detail once nothing has moved for a moment.
view.moved = () => {
  if (!view.moving) { view.moving = true; lastDraw = 0; }
  if (view.scale !== moveScale) setScale(moveScale);
  clearTimeout(settleTimer);
  settleTimer = setTimeout(() => { view.moving = false; setScale(idleScale); view.requestDraw(); }, SETTLE_MS);
  view.requestDraw();
};

Sobrevivir a un reinicio. Si aun así el navegador reinicia la GPU, la pedimos de vuelta, reconstruimos el shader y seguimos a menor resolución en lugar de congelarnos.

canvas.addEventListener("webglcontextlost", (e) => {
  e.preventDefault();   // "I want it back"
  view.lost = true;
  paintOverlay();
});
canvas.addEventListener("webglcontextrestored", () => {
  view.lost = false;
  idleScale = 1;        // come back at a lower resolution
  initGL();
  setScale(idleScale);
  view.requestDraw();
});

¿Por qué no bajar simplemente las iteraciones?

Ayuda, pero cambia la imagen: con pocas iteraciones el borde se llena de negro falso. Bajar la resolución solo cambia la nitidez, y solo mientras te mueves.

9. Zoom profundo con perturbación

El muro de precisión del paso 6 tiene que ver con los números, no con la velocidad de la GPU. Un float32 tiene unas 7 cifras significativas. A ×10¹⁰, dos píxeles vecinos difieren en la cifra 12, así que todos reciben la misma c. Desmarca Zoom profundo en "Inmersión ×10¹⁰" y la imagen se vuelve un solo color plano.

La solución, que encontró K. I. Martin en 2013, es dejar de pedirle a cada píxel su posición completa. Piensa en una foto de satélite. Medir la latitud exacta de cada casa con diez decimales es difícil. En cambio, mides un punto de referencia con mucha precisión y, para cada casa, solo "12 metros al este de la referencia". Las distancias pequeñas son fáciles, y la parte precisa se hace una sola vez.

Aquí la referencia es la órbita de referencia: el centro de la vista, iterado en la CPU con float64 (16 cifras), algo barato porque es un solo punto. Después cada píxel solo sigue su diferencia con esa órbita. Esas diferencias son diminutas, y los números diminutos son justo lo que los float guardan bien: un float32 puede guardar 10⁻³⁰ con sus 7 cifras intactas.

// Reference orbit: one point (the centre of the view) iterated in float64 on the CPU.
function referenceOrbit(cx, cy, maxIter) {
  const zs = [0, 0];
  let x = 0, y = 0;
  for (let i = 0; i < maxIter; i++) {
    [x, y] = [x * x - y * y + cx, 2 * x * y + cy];
    if (x * x + y * y > 65536) break;
    zs.push(x, y);
  }
  return zs;   // [X0, Y0, X1, Y1, …], uploaded to the GPU as a float texture
}

El álgebra es corta. Escribe la órbita del píxel como la referencia más una diferencia, z = Z + δ, y su punto como c = C + Δc. Mete eso en z² + c y resta el paso de la referencia, Z² + C. Lo que queda es la regla de la diferencia: δ' = 2·Z·δ + δ² + Δc. Nunca necesita los números grandes, solo Z, que viene de la textura.

// z = Z + dz, c = C + dc. From z² + c: dz' = 2·Z·dz + dz² + dc. Only small numbers here.
vec2 dc = p * uHeight;
vec2 dz = vec2(0.0);
float m = 0.0;
for (int i = 0; i < 5000; i++) {
  if (float(i) >= uMaxIter) break;
  dz = 2.0 * cmul(refAt(m), dz) + cmul(dz, dz) + dc;
  m += 1.0;
  z = refAt(m) + dz;
  if (dot(z, z) > bail * bail) { escaped = true; break; }
  n += 1.0;
  // Rebasing: if the pixel is nearer to 0 than to the reference, or the reference has run
  // out, restart the reference from Z0 = 0 and carry the whole value as the difference.
  if (dot(z, z) < dot(dz, dz) || m >= uRefLen - 1.0) { dz = z; m = 0.0; }
}

¿Qué son los "glitches" y qué hace el rebasing?

El truco supone que cada píxel se mantiene cerca de la referencia. A veces la órbita de un píxel pasa cerca de 0 mientras la referencia está lejos, y la diferencia se vuelve más grande que el valor que describe. Se pierde precisión, y los primeros programas con perturbación mostraban manchas planas llamadas glitches. Desde 2021 se usa el rebasing: cuando el píxel está más cerca de 0 que de la referencia, o cuando la órbita de referencia ya escapó y se acabó, la referencia vuelve a empezar (Z₀ = 0) y todo el valor pasa a ser la nueva diferencia. Cuesta una comparación.

¿Por qué se detiene en torno a ×10¹³?

Ahora el límite es la referencia y la posición de la vista, guardadas en float64 (unas 16 cifras). Para ir más hondo hacen falta números de precisión arbitraria para esa única órbita, por ejemplo BigInt de JavaScript como coma fija. La parte de la GPU no cambiaría.

Los lugares profundos necesitan más iteraciones, porque cerca del borde los puntos tardan más en decidirse. "Inmersión ×10¹⁰" sube el conteo a 2000. Eso es lo caro en un móvil, y la resolución adaptativa del paso 8 es lo que lo mantiene usable.

Para seguir

Fijar c y variar el punto de partida: eso dibuja un conjunto de Julia, uno por cada c. Es la escena siguiente: Julia.
Calcular la órbita de referencia con precisión arbitraria para pasar de ×10¹³.
Colorear según la forma de la órbita (orbit traps) en vez del tiempo de escape.

SourcesFuentes

B. B. Mandelbrot, "Fractal aspects of the iteration of z → λz(1−z) for complex λ and z", Annals of the New York Academy of Sciences 357 (1980) 249–259.
R. Brooks, J. P. Matelski, "The dynamics of 2-generator subgroups of PSL(2,C)", in Riemann Surfaces and Related Topics, Annals of Mathematics Studies 97 (1981), from a 1978 conference.
A. Douady, J. H. Hubbard, "Itération des polynômes quadratiques complexes", C. R. Acad. Sci. Paris 294 (1982) 123–126.
M. Shishikura, "The Hausdorff dimension of the boundary of the Mandelbrot set and Julia sets", Annals of Mathematics 147 (1998) 225–267. arxiv.org/abs/math/9201282
I. Quilez, "Palettes". iquilezles.org/articles/palettes/
Ultra Fractal documentation, "Perturbation calculations" (on K. I. Martin's method, first implemented in SuperFractalThing). ultrafractal.com/help/formulas/perturbationcalculations.html
R. Munafo, "Perturbation Methods", Mu-Ency. mrob.com/pub/muency/perturbationmethods.html
The Art of Code, "Making a Mandelbrot fractal explorer in Unity" (video). youtube.com/watch?v=kY7liQVPQSc
D!NG, "The Mandelbrot Set" (video). youtube.com/watch?v=MwjsO6aniig

Tutorial format inspired by Catlike Coding.Formato de tutorial inspirado en Catlike Coding.