The science
Iteration, escape and smooth colour
zn+1 = zn² + c, z0 = 0. If |zn| ever exceeds 2 the orbit escapes; c belongs to the set otherwise. The escape count n is an integer, so colouring by it gives bands. With a large escape radius the continuous count μ = n + 1 − log2(ln|zn|) removes them (it differs from the textbook normalised value only by a constant shift). Here the bailout is |z| = 256.
Distance estimate
Differentiating the iteration gives z'n+1 = 2 zn z'n + 1. The distance from c to the set boundary is about DE = |z| ln|z| / |z'| (Douady–Hubbard Green's function argument; a lower bound up to a factor 4 by Koebe). It drives the boundary glow and the thin filaments stay visible at any depth.
Analytic landmarks
The main cardioid is the set of c with an attracting fixed point: c = eit/2 − e2it/4 is its boundary, where the multiplier is eit. Where the multiplier is a root of unity e2πip/q a bulb of period q touches the cardioid; the overlay marks these root points for q ≤ 8. The period-2 bulb is the disc |c + 1| < 1/4. A Misiurewicz point such as c = i has a pre-periodic orbit (0, i, −1+i, −i, −1+i, ...) and the set around it is exactly self-similar and spiral shaped, which is why a deep zoom of it never runs out of structure.
Why double precision stops at about 10-15
A float64 has 53 bits. When the picture is narrower than about 10-15 neighbouring pixels have the same c and the image turns into blocks. float32 (what a phone GPU offers) fails at about 10-5.
Perturbation theory
Pick one reference point C and compute its orbit Zn once, with arbitrary precision (here BigInt fixed point, 100 bits more than the zoom needs). For a pixel c = C + δc write zn = Zn + δn. Subtracting the two iterations gives δn+1 = 2 Zn δn + δn² + δc. All quantities in this equation are small differences, so ordinary floating point is accurate to relative 10-7, however deep we are. The GPU stores δ divided by the view scale s (so values stay near 1 even at s = 10-34, far below the float32 range) and the iteration becomes δ' ← 2Zδ' + sδ'² + δc'.
Glitches and rebasing
Perturbation breaks when the pixel orbit passes much closer to 0 than the reference does, or when the reference orbit ends (it escaped): Pauldelbrot's criterion flags |Z+δ| ≪ |Z|. Instead of computing secondary references we use the rebasing trick of Zhuoran (2021): when |Zn + δn| < |δn| or the reference ends, set δ ← Z+δ and restart the reference index at 0. The rebase map view shows how many times each pixel did this.
Series approximation
Early on δn is a polynomial in δc: δn ≈ Anδc + Bnδc² + Cnδc³ with A' = 2ZA + 1, B' = 2ZB + A², C' = 2ZC + 2AB. The coefficients are computed along the reference orbit, and every pixel starts at the largest n where the cubic term stays below 10-8 of the linear term. Switch it off and the picture must not change.
Honest simplifications
- The reference orbit is computed in fixed point with 96 guard bits; a chaotic orbit could in principle lose more, so extreme depths are verified by re-rendering, not by proof.
- Pixel orbits are float32; counts can differ by ±1 at razor-thin boundaries from true values (see the Verify tab for the measured agreement).
- The depth is limited to 5×10-35 because the scaled values must stay inside the float32 range. Series approximation uses only 3 terms.
- The maximum iteration count is automatic (about 200 + 70 d1.3 for d decimal digits of depth); deeper structure may need the multiplier.
References
B. Mandelbrot (1980); A. Douady & J. H. Hubbard (1984–85) on connectedness and external rays; K. I. Martin (2013) "SuperFractalThing" introduced perturbation to this community; Pauldelbrot (2014) glitch criterion; Zhuoran (2021) rebasing; Ewing & Schober (1992) and Fisher & Hill (1993) on the area. Not affiliated with any software named here.