Mandelbrot

How to determine “in” or “out”

The Mandelbrot set is a visualization of a mathematical property of values on the complex plane. Every pixel location [x, y] (with e.g. x in [0, 1920) and y in [0, 1080)) is mapped to a complex value $c=\frac{x-1280}{540} + \frac{(y-540)}{540}i$.

Then let us iterate the sequence $z_{n+1} = z_n^2 + c$ with $z_0 = c$ until either a predefined number of iterations is reached, or the sequence tends to infinity (which is certain once $norm(z_n) = Re(z_n)^2+Im(z_n)^2 >= 2^2$). If the sequence tends to infinity, color the corresponding pixel with the number of iterations that were required to determine the condition. Otherwise color the pixel black.

In simpler terms, given cr and ci ($c = c_r + c_ii$) the iteration goes like this:

float zr = cr;
float zi = ci;
float norm = zr * zr + zi * zi;
unsigned iterations = 0;
while (norm < 4 and iterations < max_iterations) {
  // z = z² + c
  zr = zr * zr - zi * zi + cr;
  zi = 2 * zr * zi + ci;
  norm = zr * zr + zi * zi;
  ++iterations;
}
color = iterations == max_iterations ? BLACK : colorize(iterations);

Exercise

Where’s the data-parallelism we can use for SIMD?

TIP

You don’t have to start from scratch. There’s an existing (non-SIMD) implementation in a zip or tar.gz.

See below for an overview of the code.

Why doesn’t g++ -O3 -ffast-math -march=native vectorize the loop? Can you fix it? (I gave up… maybe std::execution::unseq can do it? )

TIP

GCC optimization dump options (try e.g. -fopt-info-vec-all)

Use stdx::native_simd<T> to explicitly vectorize the Mandelbrot computation.

TIP

vir::iota_v, vir::cvt, mask reductions, where expressions

Are you’re images not equal? Why? How can you make them equal? A discussion about floating-point and equality follows…

Compare speed-ups. Does this match your expectations? Try wider simd (using stdx::fixed_size_simd).

Render this: ./mandel part3.ppm 300000 -0.743643887037158704752191506114774 0.131825904205311970493132056385139 0.00000000002.

Why does it look bad? How can you fix it? Performance impact?

Hints

Starting point

1. Quickly scan over image.h and mycomplex.h.

   • Image pixels can be accessed via img[x, y] and are stored as std::uint32_t. The low 8 bits represent red, the bits 8:15 represent green, the bits 16:23 represent blue, and the high 8 bits are discarded.

   • mycomplex<T> is a very trimmed-down implementation of complex numbers, only providing the operations we need for this task. You can also use std::complex<float>, but for me at least that makes the program slower. Also std::complex<simd<float>> is not really a thing, whereas mycomplex<simd<float>> works just fine.

2. Ignore tsc.h. It provides the time_stamp_counter used in main().

3. Open up main.cpp and find the commented code at the bottom of main(). Uncomment the code and then the compiler will complain about mandelbrot.h.

4. mandelbrot.h is where the fun happens.

   • The file already includes vir-simd headers and makes the stdx convenience namespace work. It then defines two simd specializations for float and unsigned int.

   • The mandelbrot function constructs an Image and iterates over all pixels. For each pixel it determines the corresponding value for c and then iterates the sequence. The number of iterations is then turned into a color using three differently oscillating functions for the three colors.

How to write an image (write_image implements this already)

A very simple, straightforward solution is to output PBM (b/w), PGM (grayscale), or PPM (RGB)

E.g. for PGM, open a file (std::fstream) and write

constexpr int width = 1920;
constexpr int height = 1280;

file << "P5\n"                         // filetype
     << width << ' ' << height << '\n' // image size
     << "255\n";                       // max grayscale color value

for (y = 0; y < height; ++y) {
  for (x = 0; x < width; ++x) {
    unsigned char color = 256 - mandelbrot_iterations(x, y); // 0-255
    file.put(color);
  }
}

Bonus: Buddhabrot

The Buddhabrot is a derivative of the Mandelbrot iteration above. The Wikipedia section on “Rendering method” explains how it’s done. Your challenge is to find a way to apply what you learned about SIMD, ILP, and memory accesses to compute an image more efficiently.