In late November 2024, I saw this video:

It looks at the Mandelbrot set with exponents other than 2 and constructs a fractal in which the coordinate is the exponent used in the Mandelbrot set. I was looking for a new cover image for my blog at the time and decided to use the complex-exponent Mandelbrot set, but with a twist: I'd render a cloudy histogram of all the trajectories of points that escape the Mandelbrot set, a technique known as the Buddhabrot.

The above image is from Wikipedia.¹ The result looks like this:

Refresher on the Mandelbrot set

I'll assume the reader knows what the Mandelbrot set is. The mathematical formulation of the basic Mandelbrot set is as follows:

1. Start with a complex number c equal to a normalized image coordinate and a complex number z equal to 0.
2. Iterate a few hundred times the equation z = z^p + c where p is the exponent.
3. If the magnitude of z exceeds 2, the point c is not in the Mandelbrot set.
4. Color the point on the image we took c from black if it is in the Mandelbrot set, and white if it is not.

This is fairly simple to implement in Rust with the png crate for writing images and rayon for parallelizing the main image creation loop:

use std::fs::File;
use std::io::BufWriter;
use std::path::Path;
use rayon::prelude::*;

fn main() {
    let (w, h) = (512, 512);
    let pixels = (0..w * h).into_par_iter().map(|i| {
        let (y, x) = (i / h, i % h);
        let (y, x) = (y as f32 / h as f32, x as f32 / w as f32);
        let (y, x) = ((y - 0.5) * 4., (x - 0.5) * 4.);

        let mut a = 0.0;
        let mut b = 0.0;
        for _ in 0..100 {
            // (a+bi)(a+bi)
            // a^2 - b^2 + 2abi
            (a, b) = (a*a - b*b + x, 2.*a*b + y);
            if (a.abs() + b.abs()) > 2. {
                return (255, 255, 255);
            }
        }

        (5, 5, 5)
    });
    let out_file = File::create(Path::new("assets/out.png")).unwrap();
    let out_writer = BufWriter::new(out_file);
    let mut encoder = png::Encoder::new(out_writer, w, h);
    encoder.set_color(png::ColorType::Rgb);
    encoder.set_depth(png::BitDepth::Eight);
    let mut writer = encoder.write_header().unwrap();
    let data = pixels.map(|x| [x.0, x.1, x.2]).flatten().collect::<Vec<_>>();
    writer.write_image_data(&data).unwrap();
}

As you can see, this code hardcodes the complex square function and uses 100 iterations to determine if a point is in the Mandelbrot set. The result is as follows:

Adding the complex power is a bit more tricky. There are no primitives for complex exponentiation, so here is a simple derivation of the formula for $z^{xi}$ for $x \in \mathbb{R}$:

$$z^{xi} = e^{xi \ln z}$$ $$= e^{xi (\ln |z| + i\tan^{-1}(\text{Re}(z), \text{Im}(z)))}$$ $$= e^{-x\tan^{-1}(\text{Re}(z), \text{Im}(z)) + xi \ln |z|}$$ $$= e^{-x\tan^{-1}(\text{Re}(z), \text{Im}(z))} e^{xi \ln |z|}$$ $$e^{xi \ln |z|} = \cos(x \ln |z|) + i \sin(x \ln |z|)$$

In plain words, we replace complex exponentiation with complex multiplication (rotation) of a complex logarithm. Code snippet that implements this modification for $x = 4$:

//...
    // z^i + c = e^(i ln z) + c

    let mut distance = (a*a + b*b).sqrt().ln();
    let mut angle = b.atan2(a);
    (distance, angle) = (distance * 4., angle * 4.);
    (distance, angle) = (-angle, distance);
    let distance = distance.exp();
    let (a_base, b_base) = (angle.cos(), angle.sin());
    (a, b) = (distance * a_base, distance * b_base);


    (a, b) = (a + x, b + y);
// ...

However, simply plugging this into the Mandelbrot renderer will not work. To see the issue, think about what a and b are at the start of the loop and what operations are done to them:

// ...
let mut a = 0.0;
let mut b = 0.0;
// ...
    let mut distance = (a*a + b*b).sqrt().ln();
        // = (0 + 0).sqrt().ln() = (0).ln() = very bad thing

Put simply, this code will create NaN values and will likely create a pure black or white image. I tried some things like clipping the square root to a minimum value, but it only recovered some of the pixels.

Of all places, I found the solution to this issue in TikTok. A video by Dylan Woodbrey shows an implementation of the complex exponent Mandelbrot set in Desmos along with a link to the formula. I looked through the math to see where I went wrong and realized the Desmos implementation used a slightly different definition for the Mandelbrot set:

1. Start with a complex number c equal to a normalized image coordinate and a complex number z equal to 0 that is the same as c.
2. Iterate a few hundred times the equation z = z^p + c where p is the exponent.
3. If the magnitude of z exceeds 2, the point c is not in the Mandelbrot set.
4. Color the point on the image we took c from black if it is in the Mandelbrot set, and white if it is not.

In other words, the starting coordinate sidesteps the numerical issue entirely by starting the iteration from the point we are testing. The change in the code is simple:

// let mut a = 0.0;
// let mut b = 0.0;
let mut a: f64 = x;
let mut b: f64 = y;
// ...

This is mathematically equivalent to assuming the first complex exponentiation produced a zero and we simply added the coordinate to it.

Following this change, the code produces the following image:

We can see the interesting grid-like structure from the video, but the image is black-and white, which makes it kind of boring. This means it's time to add the Mandelbuddha histogram counting.

1. Start with a complex number c equal to a normalized image coordinate and a complex number z that is the same as c
2. Iterate a few hundred times the equation z = z^p + c where p is the exponent.
3. If the magnitude of z exceeds 2, the point c is not in the Mandelbrot set.
4. Color the point on the image we took c from black if it is in the Mandelbrot set, and white if it is not. If the point escaped, record all points on it and make them slightly more bright on the image.

The current architecture assumes the color for all pixels is independent of each other and can be computed from just the coordinate. This is not the case for the Mandelbuddha, as a pixel can be colored from any point that has a trajectory passing through it. What we can do instead is sample N randomly positioned points, run the Mandelbrot iteration on them, and record the trajectory of each point that escapes in an atomic histogram:

let mut atomic_pixels: Vec<AtomicU32> = (0..w * h).map(|_| AtomicU32::new(0)).collect();
(0..w*h*10).into_par_iter().for_each(|i| {
    let mut rng = rand::thread_rng();
    let (y, x): (f64, f64) = (rng.gen(), rng.gen());
    let mut visited = Vec::<usize>::with_capacity(N_ITERATIONS);
    // ...
    if (a*a + b*b) > 4. {
        for v in visited {
            atomic_pixels[v].fetch_add(1, std::sync::atomic::Ordering::Relaxed);
        }
    }
});

This can create real Mandelbuddha images, albeit black-and-white:

I did a bit of tweaking to find an exponent I liked:

The image is grainy, like path tracing renders or a low-resolution photo. I like the look, but the same technique could produce smoother renders like the one on Wikipedia.

Adding color is simple: some channels will have lower cutoffs for the number of iterations and some will have higher, with the same array of recorded points used for all channels.

At this point, I added some features like better color normalization² and cropping of the resulting images. I put an image from mandelbub on the main page of the website and generated this wallpaper: