# Mandelbrot fractal

This is a code that outputs a .ppm file with the mandelbrot fractal.

How can I optimize this?

#include<bits/stdc++.h>
using namespace std;

int findMandelbrot(double cr, double ci, int max_iterations)
{
int i = 0;
double zr = 0.0, zi = 0.0;
while (i < max_iterations && zr * zr + zi * zi < 4.0)
{
double temp = zr * zr - zi * zi + cr;
zi = 2.0 * zr * zi + ci;
zr = temp;
++i;
}
return i;
}

double mapToReal(int x, int imageWidth, double minR, double maxR)
{
double range = maxR - minR;
return x * (range / imageWidth) + minR;
}

double mapToImaginary(int y, int imageHeight, double minI, double maxI)
{
double range = maxI - minI;
return y * (range / imageHeight) + minI;
}

int main()
{
ifstream f("input.txt");
int imageWidth, imageHeight, maxN;
double minR, maxR, minI, maxI;

if (!f)
{
cout << "Could not open file!" << endl;
return 1;
}

f >> imageWidth >> imageHeight >> maxN;
f >> minR >> maxR >> minI >> maxI;

ofstream g("output_image.ppm");
g << "P3" << endl;
g << imageWidth << " " << imageHeight << endl;
g << "255" << endl;

double start = clock();

for (int i = 0; i < imageHeight; i++)
{
for (int j = 0; j < imageWidth; j++)
{
double cr = mapToReal(j, imageWidth, minR, maxR);
double ci = mapToImaginary(i, imageHeight, minI, maxI);

int n = findMandelbrot(cr, ci, maxN);

int r = ((int)sqrt(n) % 256);
int gr = (2*n % 256);
int b = (n % 256);

g << r << " " << gr << " " << b << " ";
}
g << endl;

if(i == imageHeight / 2) break;
}

cout << "Finished!" << endl;

double stop = clock();

cout << (stop-start)/CLOCKS_PER_SEC;
return 0;
}


I go until imageHeight / 2 because in photoshop I can just copy the other half.

I was thinking about logarthmic power but tried something and only works with integers... At the end I am showing the time it takes to run some data. I tried some things from the internet but nothing really worked out or maybe I didn't how to put it right..

this is the output with the following input: 512 512 512 -1.5 0.7 -1.0 1.0

and it took 0.315s to produce this.

• On my machine it gives 0.13s on release build. Why would you want less? It gets down to 0.125s on -O3 and -march=native Jun 4 '17 at 15:16
• @Incomputable because I don't want to use it with data like: "512 512 512 -1.5 0.7 -1.0 1.0". I want to use it with data like: "40000 40000 40000 -1.5 0.7 -1.0 1.0" or even more and here the time part gets tricky. Jun 4 '17 at 15:43
• Well, it might help to move from ppm to bmp. Also, avoid flushing. Jun 4 '17 at 22:23
• @Deduplicator, avoiding flushing seems to have very small effect. It seems like at least on libc++ and clang that is not an issue. Jun 5 '17 at 6:05

Your code produces a "plain" PPM file where all sample values are represented as ASCII decimal numbers (separated by white space). As a consequence, a substantial amount of time is spent with formatting and writing the data.

An alternative is to write a "binary" PPM file. It has the magic "P6" instead of "P3", and all sample values are represented as a single byte. (See PPM Format Specification for the details.)

So you would create the file and write the header with (using a more descriptive variable name for the file stream):

ofstream ppmFile("output_image.ppm", ios::out | ios::binary);
ppmFile << "P6" << endl;
ppmFile << imageWidth << " " << imageHeight << endl;
ppmFile << "255" << endl;


(endl flushes the output file in addition to writing a newline character, which is not necessary here. On the other hand, it does no harm because this part is not performance critical.)

In the inner loop, an RGB triplet is appended with

int n = findMandelbrot(cr, ci, maxN);

uint8_t rgbTriplet[3];
rgbTriplet[0] = (int)sqrt(n);
rgbTriplet[1] = 2 * n;
rgbTriplet[2] = n;
ppmFile.write((char *)rgbTriplet, 3);


and no newlines are written after each row. Truncating of int to the unsigned type uint8_t is well-defined so that the explicit remainder operation % 256 is not needed.

In my test on a MacBook this reduces the running time for the input data "512 512 512 -1.5 0.7 -1.0 1.0" from ≈0.25 seconds to ≈0.14 seconds.

• Comments are not for extended discussion; this conversation has been moved to chat. Jun 5 '17 at 12:21

I agree with most of what @user140417 said (though I agree with Deduplicator on the call to pow). So sticking strictly to performance, I have the following recommendations.

# Precompute Everything You Can

You are calculating the real and imaginary parts on every iteration of your inner loop in main(). This is costly. Every row will have the same real values and every column will have the same imaginary values. So pre-calculate 1 row of reals and 1 row of imaginaries before starting the main loop.

It's the same situation with your coloring. Since you know maxN, you can have 3 arrays - one for the red channel, one for the green channel, and one for the blue channel. Just pre-calculate all maxN possible values. Then in the inner loop, it's just a table look-up rather than a calculation.

# Use SIMD Techniques

Much of your math is computing the same loop for different sets of floating point numbers. This is exactly the kind of thing that SIMD instruction sets like SSE, AVX, etc. were invented for. You can compute 2 or 4 values per instruction, allowing you to significantly improve the speed of your implementation. You can take it even further with…

If you break the input area up into either tiles or strips, you can have multiple threads, each working on a different tile or strip at the same time. On multi-core machines this will speed things up in a nearly linear fashion (in my experience) - up to a point. Combine that with SIMD instructions and you can really get things moving!

# Use the GPU

If you really want to make things fast, you can write a GPU implementation using OpenGL compute shaders, or OpenCL kernels, or whatever other library you like to use. GPUs generally have hundreds to thousands of cores of floating point processors making this a really good fit.

• About precalculation thing, clang seems to be able to outline and precompute the constants itself, as it made no difference for me. Jun 5 '17 at 6:02
• That's entirely possible. Compiler optimizations are getting better every day! Jun 5 '17 at 6:03

### Don't include bits/stdc++.h

This is a non-portable header file that pulls in way more than needed, increasing compilation time and binary size.

### Don't use using namespace std;

This combined with the above greatly pollutes the global namespace. If you ever decide to use a third party library (for precision math, image drawing, or what have you) there is a high probability you'll run into conflicts.

### Use more descriptive variable names

All of the zr, zi, et al. is not only hard to read but easy to make a fatal typo that's hard to debug. You may as well use std::complex or at the least a std::pair. If you later decide to add the ability to zoom or pan your code readability will rapidly diminish if you continue to use those types of variable names.

### Rename mapToReal and mapToImaginary

Maybe you are using the mathematically correct definition, but from what I can tell a better name for the function would be "normalize" or "scale", as "map" is too abstract/vague.

## Optimization

### Factor out the arithmetic inside findMandelbrot

You are doing a total of (I count) 6 multiplications which can be reduced to three if you refactor the code slightly:

std::complex z{0, 0};
std::complex square_temp{0, 0};
while (square_tmp.real + square_tmp.imag <= 4.0) {
z.imag = std::pow(z.real + z.imag, 2.0) - square_temp.real - square_temp.imag;
z.imag += c.imag;
z.real = square_temp.real - square_temp.imag + c.real;
square_temp.real = std::pow(z.real, 2);
square_temp.imag = std::pow(z.imag, 2);
}
// Note code untested


• std::pow is quite a lot more heavy than a simple multiplication. Also, the scaling-funtions are identical. Jun 4 '17 at 22:07
• @Deduplicator Pow is meant to make the example code easier to read, not in production code. For one thing, z and c are terrible variable names but I'm not writing their code for them. Jun 4 '17 at 22:15