Smerity (n):

One of the most upvoted posts on the Rust subreddit has been "A fractal I rendered with rust without any external libraries". Whilst the output is beautiful there was a minor complaint - zero dependencies still included Tokio, image, and rand!

In recent days I had a desire to recreate my experience of an early foray into numerical computing. Rust seemed like an obvious candidate!

I decided to use the opportunity to replicate two aspects with zero dependencies: parallel computation and image export. The approaches used will be anything but optimal but may be of use to others writing fractal visualizers, scanline renderers, raytracers, and other graphical projects that benefit from being "from the ground up".

To follow along or run the code yourself, see the GitHub repository.

In brief: chaos and the logistic map

Rather than fractals we'll be playing with chaos by producing a plot - a bifurcation diagram, specifically the logistic map. This diagram represents the behaviour of a simple equation, \(x_{t+1} = rx_t(1 - x_t)\), repeated over and over.

For certain values of r (such as between 2.4 and 3) the resulting series of values for x converges to a single number. Beyond 3 however the result gets quite interesting. The series of values for x dance around a small but growing set of points.

With the help of the Python terminal:

>>> r = 3.2
>>> 0.5
0.5
>>> r * _ * (1 - _)
0.8
>>> r * _ * (1 - _)
0.512
>>> r * _ * (1 - _)
0.7995392
>>> r * _ * (1 - _)
0.512884056522752
>>> r * _ * (1 - _)
0.7994688034800593
>>> r * _ * (1 - _)
0.5130189943751092
>>> r * _ * (1 - _)
0.7994576185134749
>>> r * _ * (1 - _)
0.5130404310855622

Here we see that for r = 3.2 the series of values for x bifurcates (forks or splits) between 0.799 and 0.513.

What might this look like if we sweep across the range of values for r and keep track of a heat map of where the resulting points end up?

For more on the "chaotic" aspect of the logistic map, many others in the Rust community suggested Veritasium's video.

Multi-threading for speed

In graphic generation tasks such as fractal rendering or ray tracers, the task is easily sped up by adding more computational power (see: embarrassingly parallel).

Threads are the simplest way of taking advantage of multiple cores. Each thread can be given a small chunk of work and then the underlying operating system shuffles it across the physical CPUs available to you.

// Each step of the loop will farm work out to upwards of NTHR threads
for chunk in (0..RESX).collect::<Vec<usize>>().chunks(NTHR) {
    let mut threads = vec!();

    // Each thread creates a slice of the image, fills it using the `bifurcate` function, then returns the slice
    for x in chunk {
        let handle = thread::spawn(move || {
            let mut slice: Vec<u32> = vec!(0; RESY);
            bifurcate(&mut slice);
            slice
        });
        threads.push(handle);
    }

    // Once all the threads are spawned we then retrieve the `slice` result from them in order
    // `into_iter` is necessary as it may yield (T, &T, &mut T) vs iter (&T) and JoinHandles are not Copy-able
    let slices = threads.into_iter().map(|x| x.join().unwrap());
    img.extend(slices);
}

Is this better than tokio, crossbeam, or rayon? Absolutely not. Does it work well enough for our task? Absolutely.

The main disadvantage here is that we're assuming each of the blocks of work are similarly sized. For our task we know this will be true. If the tasks have a high degree of variability however we're likely to "starve" our CPUs - easy threads finish early and leave only one or two "hard" threads worth of work.

When running with N threads on an Intel i7-1065G7 processor with 4 CPUs and 8 hyperthreads:

• 1 thread: 38.9 seconds
• 2 threads: 21.4 seconds
• 4 threads: 17.1 seconds
• 8 threads: 9.8 seconds

Exporting images: the Portable GrayMap (PGM)

Exporting images as either grayscale or RGB is a basic requirement for any graphical application.

The PGM format is the simplest way of producing an image, representing the pixels in a standard space delimited text format. If you're interested in colour you can look to the Portable PixMap (PPM) as well.

Rather than performing anti-aliasing in the application itself we'll be exporting at a higher resolution and then downsampling.

Converting from PGM to PNG

To convert a PGM image to a more standard format such as PNG we can use ImageMagick. The file names are what set the format and we can specify downsampling but setting a lower resolution than what we produce.

convert -resize 2160x1440 img.pgm img.png

Depending on your image viewer this might not be necessary for local use but will be definitely needed if you intend on sharing the images you produce with others. Eye of GNOME, the default on Ubuntu, supports PPM and PGM for example.

What about rand?

In this we've provided (minimal) replacements for the tokio and imagecrates. What about rand however?

Assuming you don't need a strong random number generated, and especially if you don't have any intent of letting it within a billion bytes of cryptographic usage, you have a number of options for psuedo random number generation (PRNG).

For a simple project a PRNG such as the linear congruential generator is likely your best bet.

#[derive(Debug)]
struct LCG {
    state: u64
}

impl LCG {
    fn randint(&mut self, m: u32) -> u32 {
        // Random from [0, m] (i.e. including the end point)
        // LCG from Knuth's MMIX
        // See: https://en.wikipedia.org/wiki/Linear_congruential_generator
        let a: u64 = 6364136223846793005;
        let c: u64 = 1442695040888963407;
        self.state = a.wrapping_mul(self.state).wrapping_add(c);
        // Low bits don't have much entropy
        let low: u32 = self.state.reverse_bits() as u32;
        let high: u32 = (self.state >> 32).reverse_bits() as u32;
        (low ^ high) % (m + 1)
    }
}

fn main() {
    // Note: the PRNG is deterministic unless you modify this state on each run
    let mut lcg = LCG { state: 42 };
    println!("Random dice rolls: {}, {}", lcg.randint(6 - 1) + 1, lcg.randint(6 - 1) + 1);
}

Conclusion

In this code we've approximated the functionality gained from importing dependencies.

It's reassuring to know we can recreate such features ourselves even if that's rarely necessary. Beyond understanding how something works from the ground up (i.e. educational purpose) this may be necessary when deploying your code in a particularly limited environment / requirement such as no_std.

Luckily Rust provides an easy ecosystem to pull in dependencies through Crates - and I'd strongly recommend that instead!