![Mandelbrot](assets/mandelbrot-banner.jpg)
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… maybestd::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
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
(usingstdx::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
-
Quickly scan over
image.h
andmycomplex.h
.-
Image pixels can be accessed via
img[x, y]
and are stored asstd::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 usestd::complex<float>
, but for me at least that makes the program slower. Alsostd::complex<simd<float>>
is not really a thing, whereasmycomplex<simd<float>>
works just fine.
-
-
Ignore
tsc.h
. It provides thetime_stamp_counter
used inmain()
. -
Open up
main.cpp
and find the commented code at the bottom ofmain()
. Uncomment the code and then the compiler will complain aboutmandelbrot.h
. -
mandelbrot.h
is where the fun happens.-
The file already includes vir-simd headers and makes the
stdx
convenience namespace work. It then defines twosimd
specializations forfloat
andunsigned int
. -
The
mandelbrot
function constructs anImage
and iterates over all pixels. For each pixel it determines the corresponding value forc
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.