Fractal Generator

What It Is

Every day at midnight EST, the site generates one fractal image and caches it as a PNG. The same date always produces the same image — there are no parameters, no user inputs, and no randomness beyond what is seeded by the date itself. The only decision the system makes is which of 28 fractal algorithms to apply, which pre-catalogued region of the complex plane to examine, and which of 20 colour palettes to map the result through.

This document explains the mathematical theory behind each of the four main fractal families in the system, how they are implemented efficiently with NumPy, how colour is applied, and the engineering decisions behind the seeding and caching infrastructure.

The Complex Plane

All fractal algorithms in this system operate on the complex plane ℂ. A complex number c = a + bi maps directly to a 2D point: the real part a gives the horizontal position and the imaginary part b gives the vertical. The entire 800×800 pixel image corresponds to a rectangular window in ℂ, and each pixel is tested independently.

In code, NumPy represents this naturally: a 2D array of complex64 values, one per pixel, is constructed from two linspace arrays and the j operator:

x = np.linspace(xmin, xmax, 800)
y = np.linspace(ymin, ymax, 800)
X, Y = np.meshgrid(x, y)
C = X + 1j * Y      # shape (800, 800) of dtype complex128

This single array represents 640,000 complex numbers simultaneously. Every subsequent operation is applied to the entire array at once via NumPy's vectorised arithmetic — no Python-level loop over pixels exists in any of the generators.

Escape-Time Algorithms

The majority of the 28 fractal types are built on the escape-time principle: for each pixel, apply a recurrence repeatedly and record how many iterations pass before the value exceeds some threshold. That count becomes the pixel colour.

The Mandelbrot Set

The canonical escape-time fractal. For each point c ∈ ℂ, define the orbit of zero under the quadratic map:

The Mandelbrot set M is the set of all c for which this orbit remains bounded forever. Visualising how quickly orbits escape — rather than just whether they do — produces the familiar banded structure at the boundary.

The escape condition is |z| > 2. This threshold has a clean mathematical justification: if |z| > 2 and |z| > |c|, then |z²| = |z|² > 2|z| > |z| + |c| ≥ |z² + c| − |c|, which can be used to show the orbit diverges to infinity. An escape radius larger than 2 would still be correct but would slow convergence for no visual benefit at standard zoom levels.

The iteration limit is 256. Pixels that have not escaped after 256 iterations are treated as belonging to the set (coloured black in most palettes). The iteration count M[r][c] for each pixel is then used to index into a 256-entry colour lookup table.

Julia Sets

Julia sets are closely related to the Mandelbrot set but are parameterised differently. Instead of varying c per pixel and starting from z = 0, the Julia set for a fixed c ∈ ℂ is computed by varying the starting point z across the plane and using the same constant c for every pixel:

The Mandelbrot set and Julia sets are dual to each other in a precise sense: the Julia set for a given c is connected if and only if c belongs to the Mandelbrot set. Points deep inside M produce solid-looking Julia sets; points near the boundary produce intricate connected structures; points outside produce "dust" (totally disconnected sets).

The system catalogues 12 named Julia set constants chosen for visual richness — the Douady rabbit (c = -0.8 + 0.156i), the Siegel disk (c = 0.285 + 0.01i), the dendrite (c = -0.4 + 0.6i), and others. The seed selects one and applies a small perturbation so that the same c-value doesn't produce an identical image on every occurrence.

The Vectorised Iteration Loop

The core iteration is the same for all escape-time types. The key implementation detail is the boolean mask that prevents already-escaped pixels from being updated:

Z = np.zeros_like(C)       # all orbits start at 0 (Mandelbrot) or C itself (Julia)
M = np.zeros(C.shape, dtype=int)

for i in range(256):
    mask = np.abs(Z) <= 2          # pixels still inside escape radius
    Z[mask] = Z[mask]**2 + C[mask]  # apply recurrence only to unescaped pixels
    M[mask] = i                      # overwrite iteration count until pixel escapes

When a pixel first exceeds the escape radius, mask becomes False for that index and it is never written again — its final value in M is the iteration count at which it escaped. The loop continues for all 256 iterations even if most pixels have already escaped, because NumPy's vectorised operations have no early-exit mechanism. The mask operations themselves are cheaper than branching, so this is still significantly faster than any pixel-by-pixel Python loop.

Escape-Time Variants — One-Line Changes to the Recurrence

Most of the 28 fractal types are structurally identical to Mandelbrot: the same vectorised loop, the same escape condition, the same output format. The only difference is a single line — the recurrence applied on each step. This table shows how small mathematical modifications produce completely different geometries:

The Burning Ship is particularly interesting: taking absolute values of the real and imaginary parts before squaring is equivalent to folding the complex plane onto its positive quadrant before each iteration. This introduces reflective symmetry and produces the characteristic "ship" outline with sharp angular features that the smooth Mandelbrot never exhibits.

The Celtic variant is more subtle. Expanding z² = (a+bi)² = (a²-b²) + 2abi, the Celtic map takes the absolute value of only the real part of z²: the recurrence is (|a²-b²| + 2abi) + c. This preserves the imaginary part of the square unchanged while flipping the real part across zero, creating a fractal "knot" structure that is related to but visually distinct from both Mandelbrot and Burning Ship.

The Tricorn uses the complex conjugate z̄ = a - bi instead of z before squaring. This breaks the Cauchy-Riemann equations — the map is no longer complex-differentiable (holomorphic) — which is why the Tricorn has fundamentally different topology from the Mandelbrot set, including the characteristic three-pronged shape that gives it its name.

Memory Fractals: Phoenix and Man o' War

These two types extend the standard orbit by carrying additional state between iterations.

The Phoenix fractal adds a feedback term using the previous iteration's z value, introducing a form of short-term memory into the orbit:

where p is a real constant (typically small). This second-order recurrence creates spiral structures and "phoenix-wing" patterns that don't appear in first-order maps.

The Man o' War fractal goes further: both z and c are updated each iteration. Starting with C_0 = C (the pixel coordinates), the recurrence is:

C is no longer constant — it accumulates the orbit history. This means the effective "parameter" shifts as the iteration progresses, producing a richer but less geometrically regular structure.

Transcendental Fractals

Several types replace the polynomial z² with transcendental functions. The recurrences are:

The larger escape radii are needed because these functions grow much faster than polynomials — e^z can reach enormous magnitudes quickly along the real axis, so the threshold of 2 would incorrectly classify most of the plane as "escaped" before any structure forms. The sine fractal uses a different escape criterion entirely: |\text{Im}(z)| > 50, because \sin(z) grows without bound in the imaginary direction but not along the real axis.

The Collatz fractal is a smooth analytic extension of the Collatz conjecture map to the complex plane. The standard Collatz map on integers is f(n) = n/2 if n is even, 3n+1 if odd — the smooth version interpolates between these using a cosine:

When z is an even integer this reduces to z/2; when odd, to (3z+1)/2. The resulting fractal has a completely different character from polynomial types — irregular, coral-like structures without the smooth boundary gradients of the Mandelbrot set.

Zubieta Fractal — Rational Maps with Poles

The Zubieta fractal uses a rational map (ratio of polynomials) rather than a polynomial: z ← z² + c/z. This introduces a pole at z = 0: points close to the origin are mapped to very large values. The implementation must guard against division by zero:

safe = mask & (np.abs(Z) > 1e-10)
Z[safe] = Z[safe]**2 + C[safe] / Z[safe]   # only update away from the pole

The pole creates "tearing" structures in the fractal that polynomial maps cannot produce — circular holes and radial filaments emanating from the origin.

Newton Fractals — Convergence-Based Colouring

Newton fractals use a completely different criterion: instead of checking whether an orbit escapes to infinity, they apply Newton's root-finding method and check which root the orbit converges to.

Newton's Method in the Complex Plane

Given a polynomial f(z), Newton's method iterates:

For f(z) = z³ − 1, the three roots are the cube roots of unity: 1, e^{2πi/3}, and e^{4πi/3}. The derivative is f′(z) = 3z², so the step simplifies to:

Each pixel in the plane is a starting point for this iteration. After enough steps, it will converge to one of the three roots (or occasionally oscillate near a critical point). The pixel colour encodes two pieces of information: which root it converged to (the colour region — three broadly distinct hues), and how many steps it took (the shading within that region).

Implementation: Tracking Two Arrays

Unlike escape-time fractals, which only need a single integer array for iteration counts, Newton fractals require tracking which root each pixel converges to. The implementation uses two output arrays and updates them as soon as convergence is detected:

roots    = [1, np.exp(2j*np.pi/3), np.exp(4j*np.pi/3)]
root_map = np.full(Z.shape, -1, dtype=int)    # which root (-1 = not yet converged)
iter_map = np.zeros(Z.shape, dtype=int)         # how many steps to convergence

for i in range(256):
    denom = 3 * Z**2
    safe  = np.abs(denom) > 1e-10            # guard the critical point z=0
    Z[safe] = Z[safe] - (Z[safe]**3 - 1) / denom[safe]

    for r_idx, root in enumerate(roots):
        converged = (np.abs(Z - root) < 1e-6) & (root_map == -1)
        root_map[converged] = r_idx
        iter_map[converged] = i

# encode root region + convergence speed into a palette-indexable value
M = root_map * (MAX_ITER // 3) + (iter_map % (MAX_ITER // 3))
M = np.clip(M, 0, MAX_ITER - 1)

The encoding maps root 0 → palette range [0, 85], root 1 → [86, 170], root 2 → [171, 255]. Within each range, the convergence speed provides fine shading. Pixels that never converge (e.g. starting on the boundary itself) are assigned root_map = 0 by the final clipping operation.

Nova Fractal — Newton's Method with a Perturbation Term

The Nova fractal extends the basic Newton iteration by adding a per-pixel constant c to each step:

This is equivalent to asking: "for which starting points does Newton's method, when perturbed by this constant amount each step, still converge?" The perturbation breaks the perfect three-fold symmetry of the standard Newton fractal and creates spiral structures winding around the roots. The system also implements a Newton Z⁴ variant using f(z) = z⁴ − 1, which has four roots and produces 4-fold symmetric basins.

IFS Fractals — The Chaos Game

The Barnsley Fern, Sierpiński Triangle, and Dragon Curve use a fundamentally different algorithm called the chaos game or iterated function system (IFS). Rather than testing every pixel in a grid, the algorithm iterates a single point through a random sequence of affine transformations and plots each landing position. The set of all landing positions, as the number of iterations approaches infinity, converges to the attractor of the IFS.

Affine Transformations

An affine transformation in 2D is a function of the form:

where A is a 2×2 matrix and b is a translation vector. Each transformation is a combination of rotation, scaling, and shearing — a contraction mapping that shrinks distances by some factor less than 1. The Barnsley Fern uses four such maps:

These four transformations were chosen by Michael Barnsley through a process of collage theorem: find contractions whose union of images covers the target shape. The stem transform is a pure contraction toward the base; the main frond transform is a near-identity with a slight rotation and scaling, reproducing the self-similar nature of a real fern leaflet.

Why the Chaos Game Converges

The mathematical guarantee is the Contractive Mapping Theorem (Banach Fixed Point Theorem) applied to sets. Each transformation T_i is contractive, meaning there exists s < 1 such that d(T_i(x), T_i(y)) ≤ s·d(x, y) for all points x, y. The union of a finite collection of such contractions has a unique fixed-point attractor — a set A such that A = T₁(A) ∪ T₂(A) ∪ T₃(A) ∪ T₄(A). Starting from any point and applying randomly chosen transforms (with the specified probabilities), the orbit converges to this attractor with probability 1.

In practice, 500,000 iterations is sufficient to densely approximate the attractor. The seed introduces a slight tilt variation by adjusting the off-diagonal entries of T₂'s matrix by a small amount proportional to seed % 7 * 0.003, producing ferns with slightly different leaflet angles on different days.

Rendering: Log-Scaled Density Grid

After generating 500,000 points, each point is mapped to a pixel coordinate and accumulated into a hit-count grid. The distribution is highly non-uniform: the stem receives only 1% of points while the main frond body receives 85%, creating a density ratio of roughly 85:1. Rendering the raw counts directly would make the stem nearly invisible while oversaturating the frond.

The solution is log-scaling: replacing each count n with log(1 + n) before normalising to the palette range. This compresses the high-density regions and amplifies the low-density ones, making all structural features visible simultaneously. The +1 prevents taking the log of zero for empty pixels.

# map the point cloud to a pixel density grid
M = np.zeros((HEIGHT, WIDTH), dtype=float)
np.add.at(M, (yi, xi), 1)        # accumulate hits at each pixel

M_log = np.log1p(M)               # log(1+n) — preserves zero cells as zero
M     = (M_log / M_log.max() * (MAX_ITER - 1)).astype(int)   # normalise to [0, 255]

The Sierpiński Triangle and Dragon Curve use the same framework with different sets of affine transforms. The Sierpiński Triangle is three transforms that each map the triangle onto one of its three half-size sub-triangles. The Dragon Curve (Heighway dragon) uses two transforms corresponding to the two half-size copies that tile the full dragon.

Rational Maps — Magnet Fractals

The Magnet I fractal uses a rational map derived from models of phase transitions in ferromagnetism (hence the name). The recurrence is:

This comes from the renormalisation group equations for a 2D Ising model near its critical temperature, with c representing the temperature-like parameter. The squaring of a rational function makes the denominator a source of complexity — pixels near its zeros are mapped to large values, creating additional structural features beyond simple polynomial escape-time fractals. The escape radius is 4 rather than 2 because the squaring step can produce values larger than 2 even for bounded orbits near the real axis.

Colour Mapping

The output of each fractal generator is a 2D integer array M with values in [0, 255]. Applying a colour palette is a single NumPy fancy-index operation: image = palette[M], which maps each of the 640,000 integers to an RGB triple simultaneously.

Palette Construction

Twenty named palettes are defined as lists of 5–6 RGB anchor colours. At startup (and on each generation call), the selected palette is expanded from its anchors to a full 256-entry lookup table by linear interpolation between consecutive anchor pairs:

# expand anchor colours to 256-entry LUT
palette = np.zeros((256, 3), dtype=np.uint8)
for i in range(256):
    idx    = (i / 256) * (len(anchors) - 1)
    lo, hi = int(idx), min(int(idx) + 1, len(anchors) - 1)
    t      = idx - lo                               # interpolation parameter in [0,1)
    palette[i] = anchors[lo] * (1 - t) + anchors[hi] * t

The interpolation distributes the anchor colours evenly across the 256 steps, so a palette with 5 anchors divides 256 into 4 bands of 64 entries each. The result is a smooth gradient from the first anchor to the last — or a cyclic one if the first and last anchor are the same colour.

Why Iteration Count Produces Structure

Near the boundary of any escape-time fractal, small changes in the starting position produce large changes in how long the orbit takes to escape. Points just inside the boundary may take 200 iterations; points immediately adjacent may escape in 5. This extreme sensitivity is what creates the intricate banding and spiral features — the boundary itself is the locus of maximum iteration count variation, and mapping that variation through a colour gradient makes the structure visible.

Far from the boundary (well inside or well outside the set), iteration counts are either uniformly high (inside) or uniformly low (outside). These regions appear as solid blocks of a single colour regardless of palette choice. All visual interest is concentrated in a narrow band around the boundary.

Seeding and Caching

Deterministic Date Seeds

The date string (e.g. "2026-03-30") is passed through MD5 and the first 8 hex digits are interpreted as a 32-bit integer. MD5 is not used here for its cryptographic properties (collision resistance, pre-image resistance) but purely as a hash function that maps strings to integers with good distribution — similar inputs produce dissimilar outputs, so adjacent dates produce unrelated images rather than slowly evolving sequences.

This seed then drives two separate randomness sources. First, the fractal type is selected deterministically via seed % len(fractal_types). Second, np.random.seed(seed) is called inside each generator, which seeds NumPy's Mersenne Twister PRNG. All subsequent calls to np.random.rand() within that generator are deterministic — the zoom variation, centre offset, and Julia set perturbation are all fixed by the date.

Cache Design

Each generated image is saved as cache/YYYY-MM-DD.png. A sidecar JSON file records metadata for each date: the fractal type name, palette name, seed value, and filename. On any request for a given date, the cache file is checked first — if it exists, the image is served directly from disk. Generation only occurs once per date, regardless of traffic.

The cache is persistent across restarts. With 28 types and generation taking roughly 1–3 seconds on the server, backfilling the full archive for a year of past dates would take around 20 minutes — manageable as a one-time operation if needed. The current deployment has been running long enough that the cache contains several hundred images.

There is no cache invalidation or expiry. Once a date's image is generated, it never changes. This is intentional: the whole point of the system is that any given date always shows the same fractal.

Smooth Colouring

The integer iteration count produces visible "banding" — discrete colour steps at the boundaries between different escape counts. The standard fix is smooth colouring (also called continuous iteration count or the normalised escape count method), which produces a fractional escape value rather than an integer one.

The idea comes from the observation that once a point escapes, the final value of |z| at that point carries additional information about how close it was to escaping earlier. Specifically, using the formula:

where i is the integer escape count, |z| is the magnitude of z at escape, and R is the escape radius (2 for Mandelbrot), gives a real-valued escape count that varies smoothly across the escape boundary. The double logarithm accounts for the quadratic growth of the Mandelbrot map — for degree-2 polynomials, the potential function near infinity behaves like log log |z|.

In NumPy, this can be applied after the main iteration loop by using the final magnitudes of escaped pixels:

# after the main escape-time loop, smooth the integer count for escaped pixels
escaped = M > 0   # pixels that actually escaped (not stuck at max iteration)
Z_mag   = np.abs(Z)
Z_mag   = np.where(Z_mag < 1.0, 1.0, Z_mag)   # clamp to avoid log(0)

smooth  = M.astype(float)
smooth[escaped] = (M[escaped]
    - np.log2(np.log2(Z_mag[escaped]))
    + np.log2(np.log2(2.0)))

The result can then be mapped through the palette using linear interpolation between the two palette entries at floor(smooth) and ceil(smooth), weighting by the fractional part. This completely eliminates visible banding at the cost of slightly more expensive palette application.

The current system uses integer iteration counts rather than smooth values — the banding is a deliberate aesthetic choice that gives the images a more graphic, almost data-visualisation feel. The smooth colouring approach is available in the codebase as an optional mode.

Hybrid and Variant Types

Several fractal types in the system combine ideas from two families or introduce structural variations that don't fit neatly into the main categories.

Burning Julia

The Burning Julia is a Julia set using the Burning Ship recurrence rather than the standard quadratic map. Like a standard Julia set, a fixed complex constant c is used for all pixels, and the starting value z varies across the plane. Unlike standard Julia sets, the Burning Ship absolute-value fold is applied at each step:

This combines the connectivity structure of Julia sets (which depends on whether c is inside or outside the Burning Ship set) with the sharp angular geometry of the Burning Ship recurrence. The visual result is Julia-like connected filaments with the characteristic "hard edges" of the Burning Ship rather than the smooth curves of standard Julia sets.

Glynn Fractal

The Glynn fractal is a Julia set with a non-integer exponent. The recurrence is:

where the exponent 1.5 gives the Glynn fractal its characteristic 3-lobe structure (a p+1 = 2.5 ≈ 2.5-fold structure, rounding to 3 lobes). Non-integer exponents are implemented using z^p = e^{p \log z} in the complex domain. NumPy's complex power handles this correctly for the principal branch of the logarithm, though the branch cut along the negative real axis introduces a visible seam in some orientations:

# non-integer power via complex logarithm: z^p = exp(p * log(z))
# np.power handles this correctly for complex z
Z[mask] = np.power(Z[mask], 1.5) + c   # c is a fixed complex constant for Glynn

Lambda Fractal

The Lambda fractal uses a different form of parameterisation. Rather than the standard additive form z^n + c, it uses a multiplicative one:

This is the complex quadratic family in a different coordinate. By expanding: λz(1-z) = λz - λz², and making the substitution z = ½ + w/λ, it can be put into the standard form w² + c with c = λ/2 - λ²/4. The Lambda form has a fixed point at z = 1 - 1/λ and another at z = 0. The fractal structure arises at the boundary between initial values that converge to each fixed point and those that diverge, coloured by the iteration count before exceeding the escape radius of 4.

Mandelbrot Power Variants

The system includes fractional-power Mandelbrot variants using powers 1.5, 2.5, and 3.5. Like the Glynn fractal, these use np.power(z, p) for non-integer p, which computes via the principal complex logarithm. The intermediate value p = 2.5 produces a fractal with a "melted" structure between the 2-fold Mandelbrot and the 3-fold Multibrot, with irregular asymmetry caused by the branch cut of the non-integer power.

Scheduling and Serving

APScheduler Integration

The fractal is generated each day by an APScheduler job running inside the Quart application process. The scheduler fires at 00:05 EST (05:05 UTC) — slightly after midnight so any timezone offset or clock drift doesn't accidentally generate on the wrong date. The job generates the fractal for today's date if it hasn't been cached yet, then sleeps until the next day.

from apscheduler.schedulers.asyncio import AsyncIOScheduler

scheduler = AsyncIOScheduler(timezone='America/New_York')

@scheduler.scheduled_job('cron', hour=0, minute=5)
async def scheduled_generate():
    today = datetime.now('America/New_York').strftime('%Y-%m-%d')
    path = os.path.join(CACHE_DIR, f"{today}.png")
    if not os.path.exists(path):
        # generate runs in a thread pool to avoid blocking the event loop
        await asyncio.to_thread(generate_fractal, today)

scheduler.start()

Since generate_fractal is CPU-bound (NumPy operations on a large array), it's dispatched via asyncio.to_thread so it doesn't block the Quart event loop during the seconds it takes to compute. HTTP requests continue to be served normally while generation is in progress.

API Routes

The blueprint exposes three routes. The metadata endpoint returns the fractal type, palette name, and seed for a given date. The image endpoint serves the cached PNG directly. The history endpoint lists available dates in the cache, which the front end uses to build the calendar view.

@fractal_bp.route('/fractal/api/today')
async def api_today():
    today = datetime.now().strftime('%Y-%m-%d')
    _, meta = get_fractal(today)
    return jsonify({'date': today, 'image_url': f'/fractal/api/image/{today}', **meta})

@fractal_bp.route('/fractal/api/image/<date_str>')
async def api_image(date_str):
    if not re.match(r'^\d{4}-\d{2}-\d{2}$', date_str):
        return "Invalid date", 400
    filepath, _ = get_fractal(date_str)
    return await send_file(filepath, mimetype='image/png',
                           max_age=86400)   # cached images never change

@fractal_bp.route('/fractal/api/history')
async def api_history():
    cached = sorted(
        f.removesuffix('.png') for f in os.listdir(CACHE_DIR) if f.endswith('.png')
    )
    return jsonify(cached)

Performance Considerations

The system renders at 800×800 pixels with 256 maximum iterations. This means the escape-time loop runs at most 256 × 640,000 = 163,840,000 complex operations. NumPy's vectorised operations execute in C, making this feasible in a few seconds on a single server CPU despite the scale. A Python-level loop over pixels would be approximately 100× slower and would take minutes.

The IFS fractals (Barnsley Fern, Sierpiński, Dragon Curve) use 500,000 random iterations of a sequential algorithm and cannot be vectorised in the same way — each iteration depends on the previous result. They run in a Python loop and are consequently slower than the escape-time types, taking 5–10 seconds on the server. Since generation is cached after the first run, this cost is paid at most once per date.

The resolution and iteration count were chosen empirically. At 800×800, the image is large enough for detailed structure to be visible without being impractically slow. Increasing to 1024×1024 would roughly double the computation time for comparable visual improvement. The 256 iteration limit is a standard choice for escape-time fractals at moderate zoom levels; very deep zooms would require higher limits to avoid the set boundary "collapsing" into a single dark stripe.