# Plotting the Mandelbrot set efficiently

This is my Mandelbrot set program - it includes Smooth Coloring, perodicity checks, and my approach to biomorphs.

I know of Edge Detection, but I couldn't find a simple explanation that I can understand. Any help appreciated!

How do I optimize the Mandelbrot set plotting algorithm implemented in escapesmooth() and display() further?

#include <GL/glut.h>
#include <iostream>
#include <cmath>
float lerp(float a,float b,float t){return a+t*(b-a);}
const int WIDTH =800;
const int HEIGHT=600;
float r1,r2,g1,g2,b1,b2;
bool biomorph=false;
double xPos=-3.2;
double yPos=-2.0;
double zoom=150;
double startX,startY,start_X,start_Y;
int maxIter =50;
int maxIter2=50;
double zoomFactor=1.0;
float period=0;
bool render=true;
float escapesmooth(double real, double imag){ // Optimized Escape Time Algorithm + smooth colors + perodicity checks
double r0=real;
double i0=imag;
float iter=0.0;
double tempreal,io,ro;
while(iter<maxIter&&real*real+imag*imag<=16){
tempreal=real*real-imag*imag+r0;
imag=2*real*imag+i0;
real=tempreal;
iter++;
if(real==ro&&imag==io){
iter=maxIter;
break;
}
//interesting effect
if(biomorph){
if(abs(real)>0.75&&abs(imag)>0.75){
iter=real+imag;
break;
}
}
period++;
if(period>20){
period=0;
ro=real;
io=imag;
}
}
if(iter<maxIter){
float log_zn=log((real*real)+(imag*imag))/2;
float nu=log(log_zn/log(2))/log(2);
iter=iter+1.0-nu;
}
return iter;
}

void display(){
if(render){
for(int x=0;x<WIDTH;x++){
glBegin(GL_POINTS);
for(int y=0;y<HEIGHT;y++){
float iterations=escapesmooth(x/zoom+xPos,y/zoom+yPos);
glVertex2i(x,y);
if(iterations==maxIter){
glColor3f(0,0,0);
}
else{
float r1=0.5+sin(floor(iterations)/10)/2;
float g1=floor(iterations)/maxIter;
float b1=1-floor(iterations)/maxIter;

float r2=0.5+sin((floor(iterations)+1)/10)/2;
float g2=  (floor(iterations)+1)/maxIter;
float b2=1-(floor(iterations)+1)/maxIter;
/*Sunburst:*//*
float r1=floor(iterations)/32;
float r2=floor(iterations+1)/32;
float g1=floor(iterations)/64;
float g2=floor(iterations+1)/64;
float b1=floor(iterations)/128;
float b2=floor(iterations+1)/128;
/*test*//*
float r1=floor(iterations)/32;
float r2=floor(iterations+1)/32;
float g1=floor(iterations)/64;
float g2=floor(iterations+1)/64;
float b1=floor(iterations)/128;
float b2=floor(iterations+1)/128;*/

float t=fmod(iterations,1);
glColor3f(lerp(r1,r2,t),lerp(g1,g2,t),lerp(b1,b2,t));
}
}
glEnd();
glFlush();}
}
//maxIter*=1.2;
}
void reshape(int w,int h){glutReshapeWindow(WIDTH,HEIGHT);}
bool isDragging = false;
int lastX,lastY;

void mouse(int button, int state, int x, int y){
double mouseX=static_cast<double>(x)/WIDTH;
double mouseY=static_cast<double>(y)/HEIGHT;
if(button==GLUT_LEFT_BUTTON&&state==GLUT_DOWN){
isDragging=true;
lastX=x;
lastY=y;
} else if(button==GLUT_LEFT_BUTTON&&state==GLUT_UP){
isDragging=false;
}
}
void motion(int x, int y) {
if (isDragging) {
double deltaX=(x-lastX)/zoom;
double deltaY=(y-lastY)/zoom;
xPos-=deltaX;
yPos+=deltaY;
lastX=x;
lastY=y;
maxIter=maxIter2;
glutPostRedisplay();
}
}
void key(unsigned char key, int x, int y) {
switch (key){
case 'q': {
zoom *= 1.1;
maxIter = maxIter2;
render = true;
glutPostRedisplay();
break;
}

case 'e': {
zoom *= 0.9;
maxIter = maxIter2;
render = true;
glutPostRedisplay();
break;
}
case '=': {maxIter2 += 1;maxIter = maxIter2;render = true;glutPostRedisplay();break;}
case '+': {maxIter2 += 1;maxIter = maxIter2;render = true;glutPostRedisplay();break;}
case '-': {maxIter2 -= 1;maxIter = maxIter2;render = true;glutPostRedisplay();break;}
case 'b': {if(biomorph){biomorph=false;}else{biomorph=true;}glutPostRedisplay();break;}
case 'r': {double xPos=-3.2;double yPos=-2.0;double zoom=150;maxIter = 50;render = true;glutPostRedisplay();break;}
}
}

int main(int argc, char** argv){
glutInit(&argc,argv);
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
glutInitWindowSize(WIDTH,HEIGHT);

glMatrixMode(GL_PROJECTION);
gluOrtho2D(0, WIDTH, 0, HEIGHT);
glMatrixMode(GL_MODELVIEW);
glClearColor(0.0,0.0,0.0,1.0);

glutDisplayFunc(display);
glutReshapeFunc(reshape);
glutMouseFunc(mouse);
glutMotionFunc(motion);
glutKeyboardFunc(key);

glutMainLoop();
return 0;
}



# Unnecessary use of OpenGL

Using OpenGL to draw individual pixels is the completely wrong thing to do. Apart from requiring OpenGL support which might not be present on all systems, this has a huge overhead.

Instead of using the rather low-level GLUT library, I recommend that you use a library like SDL to draw to a 2D image (in SDL, that is done using a Renderer) and to handle user input.

# float or double?

I see both float and double being used in the same expressions. Usually you want to pick one type and stick with it, otherwise unnecessary conversions between them will happen.

# Make better use of C++

As Toby Speight already mentioned, use std::complex<> types to store complex numbers. Most math functions will also directly work on those types, greatly simplifying your equations.

Since C++20 there is std::lerp(), so you don't have to write your own.

# Use a code formatting tool

Toby Speight already mentioned the high code density. It would be much more readable if you had consistent whitespace around operators, empty lines between functions, some more empty lines to separate blocks of code and to make the structure more clear.

Instead of doing that all manually, consider using a code formatting tool to automate this for you. Some IDEs and code editors might also have a formatting function built-in.

• There's some justification for not using std::complex when the goal is to make common-subexpression elimination between parts of the computations easier for the compiler to see. Specifically between the products in the magnitude-squared expression (r*r + i*i) and tempreal=real*real-imag*imag+r0;. That might be unnecessary worry since you can assume what will be inside a complex multiply function, but it's traditional and normal for Mandelbrot code to look like this. (Maybe that's a bad excuse!) Of course, either way you'd want to check the compiler-generated asm. Commented Sep 13, 2023 at 12:30
• When I said "normal for Mandelbrot code to look like that", I don't mean jamming the break condition into the while like while(iter<maxIter&&real*real+imag*imag<=16), or the poor choices of var names like i0 and io, just the fact that there are expressions like real*real flying around. There's certainly room for improvement. Commented Sep 13, 2023 at 12:35
• I want to say that the compiler should be able to see this due to inlining, but then again it might need -ffast-math for it to be allowed to do that. Commented Sep 13, 2023 at 12:35
• You don't need -ffast-math for two instances of real*real to CSE down to one computation, even if one was in a function being inlined, since of course you need to enable -O2 or higher optimization for it to perform non-terribly. Only the most paranoid strict-FP semantics would require it to be done twice, to raise two separate FP exceptions, and even GCC's bad default of -ftrapping-math doesn't do that. (It's broken anyway, but I don't think its goal was to preserve the exact number of FP exceptions, just which ones might have occurred or not, although it fails at that.) Commented Sep 13, 2023 at 12:39
• A code-formatting tool may be much better than nothing, but I have yet to see one that consistently makes good decisions, however customised. (A ‘good’ decision is of course one that matches how I personally would write it :-) IMHO, while it's a good first step if you have unformatted code, by all means then go through and improve it. (And don't run well-formatted code through it…) Commented Sep 14, 2023 at 15:00

First impression is that the code is very dense. Reading would be a bit easier if you give the code room to breathe!

Variable names in particular are short and uninformative (why do we have both i0 and io in the same scope, for example). Some of the longest names are real and imag, which we wouldn't even need if we included <complex> and used a std::complex type.

Though we include the C++ <cmath> header, we use the C names log, abs, etc. to call its functions. That's not portable - write std::log, std::abs, etc. in full.

Can you do anything to reduce the number of global variables? Shared state always makes it harder to reason about the behaviour of functions in a program.

Some code is unnecessarily convoluted - e.g. if(biomorph){biomorph=false;}else{biomorph=true;} where simple biomorph = !biomorph; would be clearer.

# General Observations

I would not want to maintain this code or delegate this code to someone else to maintain.

It is difficult to read and even more difficult to modify.

Before trying to optimize the code the code should be easy to read, write and maintain.

Commented out code such as the Sunburst code in the display indicates that the code is not ready for review, nor is it ready to be maintained.

# Avoid Global Variables

It is very difficult to read, write, debug and maintain programs that use global variables. Global variables can be modified by any function within the program and therefore require each function to be examined before making changes in the code. In C and C++ global variables impact the namespace and they can cause linking errors if they are defined in multiple files. The answers in this stackoverflow question provide a fuller explanation.

It is important to note that not all of the global variables are being used. All of the variables on the following lines are not used in the code.

float r1, r2, g1, g2, b1, b2;
double startX, startY, start_X, start_Y;


It is possible that some of the global variables are of the wrong type, for instance

double xPos = -3.2;


should probably be declared as

GLdouble xPos = -3.2;


for portability.

# Declare the Variables as Needed

In the original version of C back in the 1970s and 1980s variables had to be declared at the top of the function. That is no longer the case, and a recommended programming practice to declare the variable as needed. In C the language doesn't provide a default initialization of the variable so variables should be initialized as part of the declaration. For readability and maintainability each variable should be declared and initialized on its own line.

# Vertical Spacing

While there are some blank lines in the code, the first section of code is almost unreadable due to the lack of vertical spacing.

Rather than cram multiple statements onto a single line the code would be more maintainable if each statement was on a separate line.

    case '=': {maxIter2 += 1; maxIter = maxIter2; render = true; glutPostRedisplay(); break; }


would be much easier to read and maintain if it looked like this:

    case '=':
{
maxIter2 += 1;
maxIter = maxIter2;
render = true;
glutPostRedisplay();
break;
}


# DRY Code

There is a programming principle called the Don't Repeat Yourself Principle sometimes referred to as DRY code. If you find yourself repeating the same code multiple times it is better to encapsulate it in a function. If it is possible to loop through the code that can reduce repetition as well.

In the switch statement in the key() function there is duplication of code:

    case '=': {maxIter2 += 1; maxIter = maxIter2; render = true; glutPostRedisplay(); break; }
case '+': {maxIter2 += 1; maxIter = maxIter2; render = true; glutPostRedisplay(); break; }


In this case the code duplication is not necessary, because you can put both of the cases together:

    case '=':
case '+':
{
maxIter2 += 1;
maxIter = maxIter2;
render = true;
glutPostRedisplay();
break;
}


The could be more reduction of code duplication if the maxIter2 assignment was removed from the block.

    case '=':
case '+':
maxIter2 += 1;
case '-':
maxIter2 -= 1;
{
maxIter = maxIter2;
render = true;
glutPostRedisplay();
break;
}


It is possible that there could be even more reduction of duplicate code in this function.

• If you're going to have the + and = cases fall into the - case, they need maxIter2 += 2 to offset the maxIter2-- in that block, for a net increment of += 1. Also, maxIter = --maxIter2; is worth considering, although changing one variable and copying to another seem like two logically separate things that shouldn't be combined in one statement. Commented Sep 13, 2023 at 12:45

Since you want to do this efficiently - measure how long each OpenGL call takes. I wouldn't be surprised at all if that takes most of the time. Find out how to store various pixel values int a buffer and draw them with one call.

Measure how much time you actually save by checking that the starting value repeats. I wouldn't be surprised if the amount is negative.

Figure out where the main black blob is and find a formula to check whether a value is inside it. You can also check for the next black circle. If you think that takes too long do the check only if the previous pointer hit maxIter. So if you check points horizontally, you do this check after the first black point until you leave the black area again.

Don't check for escape on every iteration. Do say 8 iterations, then check. If you have escaped, redo the last 8 iterations to find exactly when you escaped. That saves lots of checks and lots of looping.

Use fma() whenever possible. It is faster and more precise (but not 100% identical, so the compiler doesn't do it automatically). fma(a, b, c) = a * b + c with a single rounding, and usually takes the same time as the multiplication on its own.

Check which is faster - fabs(x) > 0.75 or x * x > 0.75 * 0.75.

Use vector registers to calculate multiple values at the same time. And then try doing two or for sets of points at the same time.

float vs double: If you use vector registers then float is faster BUT float isn't every precise. In areas with 200 or 500 iterations, your images will suffer badly.

• The compiler may in fact use fused multiply-accumulate instructions automatically in some cases. Using -ffast-math with GCC and Clang will allow the compiler to do more agressive math optimizations that might affect precision. Commented Sep 13, 2023 at 16:37
• @G.Sliepen: GCC's default (without -ffast-math) is -ffp-contract=fast - contract a*b+c into FMA even if it's split across separate statements. (ISO C and C++ only allow contraction within one expression). Clang's default is IIRC -ffp-contract=no; unlike GCC, it supports on (within expressions but not across statements) as well as fast. -ffast-math sets fp-contract=fast, so yes probably good since I don't expect it to hurt precision in this case. Commented Sep 13, 2023 at 23:52
• IIRC, -std=c++?? might change the FP-contract default, as opposed to the default -std=gnu++11` or whatever. Commented Sep 14, 2023 at 13:30