# C++ Mandelbulb-set generator

This mandelbulb-set generator creates a .xyz file which contains all the points that are in the set which lie on a 3D grid. It does what it is supposed to.

However, the number of calculations done is proportional to the cube of the resolution*, which means that the script takes very long if the resolution is in the order of 500 or 1000.

* The resolution is the amount of points that each axis is divided into. A resolution of 5 means that the x, y and z axis are each divided by 5. Therefore, there are $$\5^3 = 125\$$ points that are calculated.

Here is the entire code:

#include <iostream>
#include <cmath>
#include <fstream>

using namespace std;

/* Mandelbulb: https://en.wikipedia.org/wiki/Mandelbulb
*
* The sequence v(n+1) = v(n)^power + c is iterated several times for every point in 3d space. v(0) and c are equal to that point.
* If the length of  v(n) does not exceed 'maxlength' (i.e. it does not diverge) then the point is an element of the mandelbulb set
* and saved to the file.
*/

int main()
{

//settings:
int iter = 5;
int resolution = 50;
double power = 8.0;
double maxlength = 2.0;
double rangemin = -1.75;
double rangemax = 1.75;

//Points:
double xpos, ypos, zpos;    //positions to cycle through
double xpoint, ypoint, zpoint;  //point to be calculated at that position
double cx, cy, cz;  //C-point
double r, phi, theta;   //spherical coordinates

//calculated variables:
double div = (rangemax - rangemin) / double(resolution);

//other variables:
double length = 0;
double density = 0; //Named 'density' since the 3d-volume has density=1 if the point is in the set and density=0 if point is not in the set

double progress, progdiv = 100.0 / pow((double(resolution) + 1), 3.0);

//User interface:
cout << endl;
cout << "+--------------------+" << endl;
cout << "|Mandelbulb generator|" << endl;
cout << "+--------------------+" << endl << endl;

cout << "Standard settings:" << endl;
cout << "Iterations:    " << iter << endl;
cout << "Power:     " << power << endl;
cout << "Maxlength: " << maxlength << endl;
cout << "rangemin:  " << rangemin << endl;
cout << "rangemax:  " << rangemax << endl;
cout << "Resolution:    " << resolution << endl;
cout << endl;

//create output file:
ofstream file ("Mandelbulb.xyz");

file << div << endl;    //first line contains distance between points for other scripts to read

//x,y,z-loop:
for(xpos = rangemin; xpos <= rangemax; xpos += div){
for(ypos = rangemin; ypos <= rangemax; ypos += div){
for(zpos = rangemin; zpos <= rangemax; zpos += div){

//Display progress in console:
progress += progdiv;
cout << "               \r";
cout << progress << "%\r";
cout.flush();

//reset for next point:
xpoint = xpos;
ypoint = ypos;
zpoint = zpos;
cx = xpos;
cy = ypos;
cz = zpos;

//Sequence loop:
for(int i = 0; i <= iter; i++)
{
r = sqrt(xpoint*xpoint + ypoint*ypoint + zpoint*zpoint);
phi = atan(ypoint/xpoint);
theta = acos(zpoint/r);

xpoint = pow(r, power) * sin(power * theta) * cos(power * phi) + cx;
ypoint = pow(r, power) * sin(power * theta) * sin(power * phi) + cy;
zpoint = pow(r, power) * cos(power * theta) + cz;

length = r;

if(length >= maxlength)
{
density = 0.0;
break;
}
density = 1.0;

}

if(density == 1.0)
{
file << xpos << "   " << ypos << "  " << zpos << endl;
}

}//zpos loop end
}//ypos loop end
}//xpos loop end

cout << endl << "Done." << endl;
return 0;
}


I hope that the comments make clear what happens, otherwise please ask. Again, the goal is to make it execute as fast as possible. I am a beginner at C++ and might not see some obvious bottlenecks. I'd like to avoid using external libraries.

• Parallelism is the answer here. Generating a fractal is an embarrassingly parallel problem. Split the space into k threads and have each of them calculate the subspace autonomously. If you take this approach to the limit then you would use GPUs, and have a thread calculate a single point. I wrote an article some time ago on fractals on CUDA (Julia set in this case) davidespataro.it/cuda-julia-set-fractals. – Davide Spataro Dec 3 '19 at 21:22
• You should definitely avoid flush() and endl. They will slow your program down since they force it to write to the file/stdout right now. For this it has to wait on the system. Without flushing, it can buffer the writing for later (it will certainly do it at some point) and then can write more data in a single write, which minimizes waiting. See also here – n314159 Dec 3 '19 at 21:40

I probably can't say much about the algorithm but I see many things which can help you improve your use of C++ and programming in general.

• ### Don't use using namespace std

Why? You pollute the whole std namespace into your program which could lead to name clashes. More about it here.
• ### Use const and constexpr

Consider these declarations:

//settings:
int iter = 5;
int resolution = 50;
double power = 8.0;
double maxlength = 2.0;
double rangemin = -1.75;
double rangemax = 1.75;


In your program you can now happily modify these parameters by accident. Since they should be constants say so with the keyword const:

const int iter = 5;
const int resolution = 50;
const double power = 8.0;
const double maxlength = 2.0;
const double rangemin = -1.75;
const double rangemax = 1.75;


If C++11 is available, even better: make them constexpr:

constexpr int iter = 5;
constexpr int resolution = 50;
constexpr double power = 8.0;
constexpr double maxlength = 2.0;
constexpr double rangemin = -1.75;
constexpr double rangemax = 1.75;


This makes these variables true constants. That means they eat up no memory. They are just aliases for the numbers.

• ### Use auto whenever possible

You don't have to mention types explicitly. Instead use auto. So we get:

constexpr auto iter = 5;
constexpr auto resolution = 50;
constexpr auto power = 8.0;
constexpr auto maxlength = 2.0;
constexpr auto rangemin = -1.75;
constexpr auto rangemax = 1.75;

• ### omit return 0 in main()

In C++, return 0 gets automatically inserted into the main function if you don't write it. So for this special case it is not necessary to write the return.

• ### Avoid declaring several variables on one line

double progress, progdiv = 100.0 / pow((double(resolution) + 1), 3.0);

//...

//Display progress in console:
progress += progdiv;


My compiler shows a warning in the last line:

variable 'progress' is uninitalized when used here

That means progress can have any value here and gets added here.

So better write two lines and initialize:

auto progress = 0.0;
auto progdiv = 100.0 / pow((double(resolution) + 1), 3.0);

• ### Avoid std::endl

This is really a very common beginners' trap. Even some books use std::endl all the time but it is wrong.

First of all, std::endl is nothing but this:

'\n' + std::flush


So it gives you a new line (What you want) and the performs a expensive flush of the output sequence which is in 99.9% of the cases not necessary. So stick to \n.

• ### Declare variables as late as possible:

This:

//reset for next point:
xpoint = xpos;
ypoint = ypos;
zpoint = zpos;
cx = xpos;
cy = ypos;
cz = zpos;


Should be:

//reset for next point:
auto xpoint = xpos;
auto ypoint = ypos;
auto zpoint = zpos;
auto cx = xpos;
auto cy = ypos;
auto cz = zpos;


Unlike in old C Standard, you can declare variables as late as possible. They should be declared when they are needed. This avoids mistakes.

Same goes for this:

auto r = sqrt(xpoint*xpoint + ypoint*ypoint + zpoint*zpoint);
auto phi = atan(ypoint/xpoint);
auto theta = acos(zpoint/r);

• ### Limit your line length

Consider this nice line:

//other variables:
double length = 0.0;
double density = 0.0; //Named 'density' since the 3d-volume has density=1 if the point is in the set and density=0 if point is not in the set


You should limit your line length to a value (80 or 100 are common) so you can read all the code without scrolling to the right. With this limit you can also open two code pages next to each other. More about it here.

Most IDEs support a visible line length marker. In mine it looks like this:

So better like this:

double length = 0.0;
//Named 'density' since the 3d-volume has density=1 if the point is in the
//set and density=0 if point is not in the set
double density = 0.0;

• ### Avoid long functions

Your main function is very long and it's hard to grasp from looking at it what is going on. Prefer shorter functions which only do one thing. The main does several things (user input, calculation etc). Split that into smaller functions. And the smaller functions into even smaller functions. This way you also have the option to run the algorithm without any user input, maybe for unit tests.

Refactored, the code looks like this:

#include <iostream>
#include <cmath>
#include <fstream>

/* Mandelbulb: https://en.wikipedia.org/wiki/Mandelbulb
*
* The sequence v(n+1) = v(n)^power + c is iterated several times for every
* point in 3d space. v(0) and c are equal to that point.
* If the length of  v(n) does not exceed 'maxlength' (i.e. it does not diverge)
* then the point is an element of the mandelbulb set
* and saved to the file.
*/

struct Parameters{
const int iter;
const int resolution;
const double power;
const double maxLength;
const double rangemin;
const double rangemax;
};

void printParameters(std::ostream &os, const Parameters &parameters);

void calculation(
std::ostream &calcOutput,
std::ostream &progressOutput,
bool progressOutputOn,
const Parameters &parameters);

bool sequenceLoop(double xpos, double ypos, double zpos, int iter,
double power, double maxLength);

bool sequence(double &xpoint, double &ypoint, double &zpoint, double power,
double maxLength, double cx, double cy, double cz);

void printParameters(std::ostream &os, const Parameters &parameters)
{
os << '\n'
<< "+--------------------+" << '\n'
<< "|Mandelbulb generator|" << '\n'
<< "+--------------------+" << '\n'
<<'\n'
<< "Standard settings:" << '\n'
<< "Iterations:    " << parameters.iter << '\n'
<< "Power:     " << parameters.power << '\n'
<< "Maxlength: " << parameters.maxLength << '\n'
<< "rangemin:  " << parameters.rangemin << '\n'
<< "rangemax:  " << parameters.rangemax << '\n'
<< "Resolution:    " << parameters.resolution << '\n'
<<'\n';
}

void calculation(
std::ostream &calcOutput,
std::ostream &progressOutput,
bool progressOutputOn,
const Parameters &parameters)
{
auto div = (parameters.rangemax - parameters.rangemin)
/ double(parameters.resolution);
//first line contains distance between points for other scripts to read
calcOutput << div << '\n';

auto progress = 0.0;
auto progdiv = 100.0 / pow((double(parameters.resolution) + 1), 3.0);

//x,y,z-loop:
for(auto xpos = parameters.rangemin;
xpos <= parameters.rangemax; xpos += div) {

for(auto ypos = parameters.rangemin;
ypos <= parameters.rangemax; ypos += div) {

for(auto zpos = parameters.rangemin;
zpos <= parameters.rangemax; zpos += div) {

if(progressOutputOn) {
//Display progress in console:
progress += progdiv;
progressOutput << "               \r"
<< progress << "%\r";
progressOutput.flush();
}

if(sequenceLoop(xpos, ypos, zpos, parameters.iter,
parameters.power, parameters.maxLength))
{
calcOutput << xpos << "   " << ypos << "  " << zpos << '\n';
}
}
}
}

if(progressOutputOn) {
progressOutput << '\n'
<< "Done." << '\n';
}
}

bool sequenceLoop(
double xpos,
double ypos,
double zpos,
double power,
int iter,
double maxLength)
{
auto xpoint = xpos;
auto ypoint = ypos;
auto zpoint = zpos;

for(auto i = 0; i <= iter; i++)
{
if(!sequence(xpoint, ypoint, zpoint, power, maxLength,
xpos, ypos, zpos)) {
return false;
}
}
return true;
}

bool sequence(
double &xpoint,
double &ypoint,
double &zpoint,
double power,
double maxLength,
double cx,
double cy,
double cz)
{
auto r = sqrt(xpoint*xpoint + ypoint*ypoint + zpoint*zpoint);
auto phi = atan(ypoint/xpoint);
auto theta = acos(zpoint/r);

xpoint = pow(r, power) * sin(power * theta) * cos(power * phi) + cx;
ypoint = pow(r, power) * sin(power * theta) * sin(power * phi) + cy;
zpoint = pow(r, power) * cos(power * theta) + cz;

return r < maxLength;
}

int main()
{
constexpr Parameters parameters{ 5, 50, 8.0, 2.0, -1.75, 1.75 };

printParameters(std::cout, parameters);

constexpr auto outputFileName = "Mandelbulb.xyz";
std::ofstream ofs (outputFileName);

calculation(ofs, std::cout, true, parameters);
}


I hope this gives some idea how to make the code more readable. I think it can be still improved further.

For a bit of speedup, I suggest to turn off the progress echoing in the console with progressOutputOn = false for the method calculation.

As mentioned in the comments, maybe you can split up the calculation into several threads. For example:

xpoint = pow(r, power) * sin(power * theta) * cos(power * phi) + cx;
ypoint = pow(r, power) * sin(power * theta) * sin(power * phi) + cy;
zpoint = pow(r, power) * cos(power * theta) + cz;


These 3 points get calculated one after each other. But they are independent of each other, so probably they can be excecuted in parallel. I sould try that out with std::async.

That's my first thought but maybe even the whole algorithm needs to be set up differently to have it ready for parallel computation.

Now that you have the code cleaned up, you can measure the time between the calculations to see how the speed is (also with other parameters).