Newton's Fractal
Fractals
Newton's method applied to z³−1, running live on a Super Nintendo.
The complex plane is partitioned into three basins of attraction — one per
cube root of unity — coloured red, green, and blue and shaded by convergence speed
(bright = fast, dark = boundary). Each tile requires a full 32-bit complex division.
Written in C and compiled with the llvm-mos-based
65816 toolchain (+mos-a16, 16-bit accumulator). Fills progressively.
Self-running demo — tiles fill progressively, 4 per frame.
Click the screen to focus and watch. Tab away and it pauses.
What it is
Newton's method finds roots of a function by iterating z ← z − f(z)/f′(z). For f(z) = z³−1 there are three roots — the cube roots of unity:
z₁ = 1 z₂ = e^(2πi/3) ≈ −0.500 + 0.866i z₃ = e^(4πi/3) ≈ −0.500 − 0.866i
Starting from any complex initial point, repeated application of the Newton step converges to one of the three roots. The basin of attraction of each root is the set of all starting points that converge to it. The boundary between basins is a fractal — infinitely detailed, with all three basins meeting at every boundary point. This demo colours each 8×8-pixel tile by which root it converges to (red / green / blue) and shades by speed (bright = few iterations, dark = many, near the fractal boundary).
Fixed-point arithmetic on 16-bit hardware
The iteration runs in Q8.8 fixed-point (scale = 256; range ≈ [−128, 128)). Each Newton step requires:
z² = (zr²−zi², 2·zr·zi) >> 8 [two __mulsi3]
z³ = (z2r·zr − z2i·zi,
z2r·zi + z2i·zr) >> 8 [two more __mulsi3]
num = z³ − 1 (numerator)
den = 3·z² (denominator)
q = num / den in Q8.8 [__divsi3 — complex division]
z ← z − q To prevent 32-bit overflow during the complex division, both the numerator and denominator are pre-scaled by >>1 before the products are computed (4× shrinkage combined), keeping the intermediate values below INT32_MAX. The step is clamped to ±2.0 to avoid runaway near z=0 where the denominator degenerates.
Compiler stress-test
| Item | What it exercises |
|---|---|
| __divsi3 | One full 32-bit signed complex division per tile — z³−1 / 3z² |
| __mulsi3 | Four 32-bit multiplications per Newton step: z²r, z²i, z³r, z³i |
| rep / sep | Mode-switch brackets widen A to 16-bit for 16×16→16 operations |
| Convergence branch | Iteration terminates early when |z−rootₙ| < ε — variable-depth loop |
An 8×8 gate grid (64 sample points in [−1.2, 1.2]²) runs Newton's method with a
20-iteration cap and folds each result into a rotating-XOR CRC. The gate asserts
corpus_result (0x4D8B) is identical on host, default 8-bit SNES,
+mos-a16, and +mos-xy16 — five-way differential. Hit
Verify fidelity above to reproduce that assertion live in this tab via the
same cycle-accurate bsnes-jg core used in CI.
Written in C with the llvm-mos 65816 toolchain and verified against
two emulators (MAME + bsnes-jg). Hit Verify fidelity to reproduce the build
gate's WRAM assert (corpus_result = 0x4D8B) live in this tab. No far
pointers — same binary checked five ways: host == default == +mos-a16
== +mos-xy16 == bsnes-jg.