Optimization Instance: Mandelbrot Set, half 1

Introduction

This put up paperwork the method of optimizing a small drawback – producing pictures of the Mandelbrot Set. Though I’ll present a speedup of ~8x over the naive implementation (there’s a half 2 that may carry this to ~100x in some circumstances), that is largely meant for example of how I am going about utilizing my data of software program design and laptop structure to optimize software program. This put up is fairly lengthy, however lots of the sections are unbiased of one another, so be happy to make use of the Contents to go to the sections of curiosity. Particularly I spend a while establishing the issue, and the true optimizing begins here.

Contents

Common notes

• This isn’t a cookbook – there isn’t a cookbook for this. Optimization is an exploration, not a recipe, it’s typically extra of an artwork than a science.
• Whereas many strategies would possibly be capable to be utilized to completely different issues, every issues can have its personal, typically distinctive, set of strategies that may be usefully utilized.
• Measure, measure, measure to know whether or not a change is an enchancment, and profile to know the place to enhance. Simply because one thing appears to be like prefer it would possibly enhance efficiency doesn’t imply that it does – I’ve had a few years of apply and I nonetheless typically go down blind alleys.
• Optimization is rarely completed (besides possibly for tiny issues, or the place an answer is completely reminiscence or I/O certain). Usually what occurs is that diminishing returns kick in, in that it begins taking an increasing number of effort to get smaller and smaller features. Generally attempting a complete completely different path would possibly yield higher outcomes.
• This isn’t an instructional paper – there was no literature evaluate, no peer evaluate, there isn’t a bibliography, and I’m not claiming that any of those concepts are novel.
• Though I exploit C++ for the examples, the strategies I exploit would work in any compiled language or assembler – C++, C, assembler, Fortran, Rust, Go, Java (if utilizing a JIT), though vectorizing requires entry to the intrinsics or a very good vectorizing compiler. They’d not provide something helpful in interpreted languages akin to Python or Javascript.
• There’s a repository with the code used for this here. The code proven right here shall be completely different because it has been simplified considerably for this put up.

Organising the issue

With a view to display the thought processes that go into optimizing code, it is very important firstly to have a properly outlined drawback to unravel. The instance chosen here’s a Mandelbrot Set generator. This has the benefits of being well-known and understood, appropriate for optimization, sufficiently small to discover the issue in an inexpensive period of time, and generates some fairly pictures. It has the drawback that there aren’t quite a bit algorithmic enhancements that may be made, and is simply too small an issue to display any attention-grabbing knowledge buildings.

Definition of the Mandelbrot Set

Given some extent (c) within the complicated airplane, it belongs within the Mandelbrot set if the recurrence relation (z_{n+1}=z_n+c, z_0=c) doesn’t diverge. It may be demonstrated that not diverging is equal to testing (|z|leq 4) for all (n geq 0).

Unpacking this a bit of,

As a result of we don’t have a manner of iterating an infinite variety of occasions, let’s decide some most variety of iterations (n_{max}).

For some extent ((c_x,c_y)), begin with (x=0, y=0) (that is equal to beginning with (x=c_x, y=c_y), it would simply require one additional iteration) and apply

(x_{n+1}={x_n}^2-{y_n}^2+c_x)
(y_{n+1}=2{x_n}{y_n}+c_y)

till (x_n^2+y_n^2>4) or (n>=n_{max}).

If (n_{max}) is reached, we are saying {that a} level is within the Mandelbrot Set, in any other case it isn’t. The upper the restrict (n_{max}), the much less the false positives.

It’s customary when displaying this in a picture to set the factors within the Mandelbrot set to be black, and use a gradient of colours to signify what number of iterations (n) it took to achieve (x_n^2+y_n^2>4). It finally ends up trying like this for (-2 leq x,y < 2):

Parameters of drawback

There are nonetheless some issues I have to nail down. Firstly, what kind of {hardware} am I going to run this on? Particularly, am I going to make use of a CPU, or a GPU with one thing like OpenCL or CUDA? In apply, this drawback is embarrassingly parallel (sure, that may be a technical time period), which makes it extremely suited to a GPU, however this put up is about CPU optimization. Particularly, I’m going to make use of x86-64, with the idea AVX2 is out there. Some notes shall be included for a way this might switch to different processors, akin to ARM, POWER or RISC-V. The actual CPU I used is a Intel(R) Core(TM) i5-4340M CPU @ 2.90GHz (Haswell).

Numeric illustration

There are two kinds of numbers that that should be processed and/or saved: the purpose coordinates (x),(y), and the iteration rely.

Illustration of (x),(y)

What kind of numerical illustration ought to be used for (x) and (y)? Choices embody 32 bit float, 64 bit double, 32, 64, 128 (or bigger) fastened level numbers saved in integers. The upper the precision, the deeper it’s attainable to zoom into the set earlier than rounding errors distort the photographs and make them grainy.

If fastened level integers are used, word that every one numbers are within the vary (-16 < x < 16), so 5 bits are wanted for the integer half, the remainder can be utilized for the fractional half.

Kind	Precision (bits)	Storage (bits)	Notes
float	~24	32	Precision will depend on distance from 0
double	~53	64	Precision will depend on distance from 0
32 bit fastened	27	32
64 bit fastened	59	64	Restricted AVX2 help
128 bit fastened	123	128	No direct CPU help

Be aware that the storage dimension isn’t going to matter till it turns into time to vectorize this – at which level the smaller sizes imply that extra might be packed right into a vector, which can dramatically assist with velocity.

I begin by ruling out float and 32 fastened level – they don’t permit significantly excessive ranges of zoom, and low zoom ranges don’t require many ranges of iteration, so a naive answer will already be very quick.

I additionally dominated out 128 bit integers or larger, as there isn’t a direct CPU help for this, and the issue turns into that of coping with excessive velocity multiprecision arithmetic, which is an attention-grabbing drawback, however takes me away from contemplating a wide range of optimizations.

64 bit integers would look to be a more sensible choice than doubles due to barely larger precision, and is definitely my choice (I favor to make use of fastened level to floating level the place possible as then I’ve a lot better management over error propagation), however there’s a wrinkle: I’m going to wish to vectorize this, and x86-64 AVX2 doesn’t straight help any manner of multiplying vectors of 64 bit integers and getting the excessive half, which might make the issue costly to vectorize (the identical is true of ARM NEON), which might carry me again to coping with multiprecision arithmetic to do multiplications.

For the remainder of this dialogue, I’ll use double. In apply this appears to restrict the workable zoom stage to about (2^{33} occasions image_width) which for 1000×1000 pictures means a most zoom of 8589934592000.

Illustration of the variety of iterations

One other factor to find out the scale of the weather of the array that iteration counts are saved in – the 2 choices are 16 bit integers and 32 bit integers. Smaller means much less reminiscence, however a smaller most variety of iterations. If the issue was reminiscence certain on entry to this array this selection is likely to be essential for efficiency, however there’s a lot numerical processing happening that it’s not going to matter right here. I punted on this determination by defining a sort iterations_t that may be modified later, and set it to be a 32 bit unsigned integer for now.

typedef uint32_t iterations_t;

First do the easy model

Step one is to put in writing a easy model. Beginning with this has a number of benefits:

• It’s a verify that the issue is properly understood,
• It may be used to supply a dataset for testing extra sophisticated variations,
• I think about it good apply to depart a commented copy in our closing code (or inside documentation or some place else accessible) in order that the algorithm is evident when doing upkeep,
• It is likely to be quick sufficient, wherein case optimization is just not even required.

Here’s a naive model of a Mandelbrot set generator written in C++:

iterations_t mandelbrot_point(double cx, double cy, iterations_t m) {
	iterations_t rely = 0;
	double x = 0, y = 0;
	whereas(rely < m && x * x + y * y <= 4.0) {
		double nx = x * x - y * y + cx;
		y = 2.0 * x * y + cy;
		x = nx;
		rely++;
	}
	return rely;
}

void mandelbrot_render_simple(iterations_t* p, double cx, double cy, 
		double zoom, int width, int top, iterations_t iterations) {
	double xs = cx - 0.5 / zoom;
	double ys = cy + 0.5 * top / (zoom * width);
	double inc = 1.0 / (zoom * width);
	for(int j = 0; j < top ; j++) {
		iterations_t* py = p + j * width;
		double y = ys - inc * j;
		for(int i = 0; i < width; i++) {
			py[i] = mandelbrot_point(xs + inc * i, y, iterations);
		}
	}
}

A word about mandelbrot_render_simple: It might be quicker to do the next for the internal loop:

		double x = xs;
		for(int i = 0; i < width; i++) {
			py[i] = mandelbrot_point(x, y, iterations);
			x += inc;
		}

This replaces xs + inc * i (an integer to double conversion, a multiply and an addition) by x += inc (addition), but when inc is far smaller than x, rounding errors will accumulate quickly. That is an instance of how slippery it may be to take care of floating level arithmetic – What Every Computer Scientist Should Know About Floating-Point Arithmetic by David Goldberg (ACM Computing Surveys 23 situation 1, 1991) ought to be required studying for anybody coping with floating level in a non-trivial manner. Additionally, the time spent in mandelbrot_render_simple is tiny in comparison with the time spent in mandelbrot_point, so there’s little or no profit to this in any case.

Arrange some exams

You will need to have some concept that our preliminary code works, and to arrange some exams to measure efficiency with.

I picked 4 areas of the Mandelbrot set, zoom set on the efficient most of 8589934592000 and max iterations of 50000, to make use of in testing, the primary three (A,B,C) to signify the form of workloads which may sometimes be discovered zooming in, and the fourth (D) to signify the slowest attainable instance.

check	(c_x)	(c_y)	black	element
A	-0.57245092932760	0.563219321276942	heaps	heaps
B	-0.57245092932763	0.563219321276842	little	heaps
C	-0.57245092932663	0.563219321276852	little	little
D	0	0	all	none

Check A

Check B

Check C

Check D

For these exams, I used clang++-13, with compiler flags

-Werror -Wall -Wextra -std=c++20 -mavx2 -mfma -mtune=haswell

Exams are run on an Intel(R) Core(TM) i5-4340M CPU @ 2.90GHz with 16GB of RAM, on Debian Linux with KDE and most functions shut down. I ran every case 20 occasions, which is sufficient to give a tough velocity comparability.

Preliminary outcomes:

	Check A	Check B	Check C	Check D
time (seconds)	973.4	68.7	41.2	3550.7
whole iterations	273760647440	19289413880	11543441620	1000000000000

Whether or not to optimize

I’ve executed plenty of tramping in New Zealand (tramping is the phrase we use for climbing and backpacking), and lots of tracks have rivers to cross with out the advantage of bridges. A part of the coaching for a way to do that safely is to first ask yourselves (plural, as a result of you aren’t doing this by yourself, proper?) the query “Will we cross?”, the purpose being that you just by no means have to cross a river. Solely after that do the opposite questions get thought-about: “The place can we cross?”, “How can we cross?”, and so on.

Optimization is analogous, that the primary query ought to at all times be “Do I optimize?”. Clearly the benefits of optimization are that code might race quicker, however there are causes to not:

• The code would possibly already be quick sufficient,
• The code won’t be the bottleneck, so optimizing it would have little impact till the bottlenecks are handled,
• It takes plenty of engineering time to do properly,
• Optimized code is usually a lot more durable to keep up.

For the needs of this dialogue, I determine to go forward with optimization.

Are there any compiler flags that may assist?

It’s good to get the compiler to do as a lot of the work as attainable, so with that in thoughts, are there any optimization flags which may assist that aren’t already lined by -O3?

One possibilty is -ffast-math for the mandelbrot_point operate solely. This can be a flag to actually be prevented except you’re positive of what you’re doing, however this specific operate is properly behaved, there aren’t any NaNs, no divisions or sq. roots, I don’t care about order of operations or associativity, and that loop by its nature is propagating rounding errors anyway, so -ffast-math might lead to it simply propagating completely different ones. Within the subsequent part, I’m going to have a look at the assembler produced by the compiler, which is an additional verify that ffast-math is protected to make use of.

Outcomes with -ffast-math

	Check A	Check B	Check C	Check D
time (seconds)	932.5	65.9	39.5	3404.8
whole iterations	273729279660	19298214380	11543330380	1000000000000

This can be a small enchancment. Be aware that the whole variety of iterations has modified barely, reflecting the completely different rounding within the calculation, however these adjustments are tiny, and may very well be prevented by specifying the order of operations in the issue, or by utilizing fastened level.

Have a look at the assembler

A long time in the past, when CPU architectures have been a lot easier and compilers weren’t pretty much as good, it could typically be price utilizing assembler for tight internal loops. On a number of the extra fashionable CPUs, significantly the x86 household with “out of order execution”, “μops” and “ports”, the machine code directions are transformed into one thing that bears little resemblance to the unique earlier than being run. The compilers have been written to output code that provides superb efficiency, and are very arduous to beat. I’m not saying it’s unimaginable to do higher – an individual can have a greater understanding of the issue being solved and could possibly take shortcuts that the compiler can’t, and an individual can do not less than in addition to the compiler by utilizing the compiler output as a place to begin, however this might solely be definitely worth the effort in excessive circumstances, and would require a lot of trial and error and testing for probably small features, and even then there’s the danger that one thing that works nice on one processor performs badly on one other processor in the identical household. It is usually a lot more durable to keep up. I selected to not go that route for this drawback, and let a compiler do this stage of lifting.

However, it is very important bear in mind at some stage of what’s going on, each on the meeting and structure stage. I can’t suggest Agner Fog’s Software Optimization Resources extremely sufficient for the x86 household of microprocessors. If there are equal sources for ARM (I do know that is unlikely, as there are much more selection in ARM chips), I might be delighted to search out out the place they’re.

For a good loop it is extremely a lot price trying on the assembler produced by the compiler (significantly after utilizing -ffast-math to verify that hasn’t had surprising penalties). This may be executed both with the -S flag with the compiler (gcc or clang, different compilers can have equivalents) or use the Compiler Explorer. The compiler explorer additionally has the benefits of displaying how completely different compilers deal with the identical code.

Doing this with mandelbrot_point with -ffast-math, the result’s (simply displaying the internal loop, and I’ve line numbered and annotated it):

0 .LBB0_3
1	vmulsd  %xmm4, %xmm4, %xmm6      x_squared = x * x
2	vmulsd  %xmm2, %xmm2, %xmm5      y_squared = y * y
3	vaddsd  %xmm6, %xmm5, %xmm7      mag_squared=x_squared + y_squared
4	vucomisd        %xmm3, %xmm7     department if mag_squared > 4.0
5	ja     .LBB0_6
6	vaddsd  %xmm0, %xmm6, %xmm6      x_squared_plus_cx = x_squared + cx
7	vaddsd  %xmm4, %xmm4, %xmm4      x_times_2 = x + x
8	vsubsd  %xmm5, %xmm6, %xmm5      x_next = x_squared_plus_cx - y_squared
9	vfmadd213sd %xmm1, %xmm4, %xmm2  y = x_times_2 * y + cy
10	incl    %eax                     rely++
11	vmovapd %xmm5, %xmm4             x =x _next
12	cmpl    %eax, %edi               department if rely != m
13	jne     .LBB0_3

This appears to be like actually environment friendly, and there a doesn’t appear to be something unhealthy launched by -ffast-math – largely it simply allowed a multiply and an add to be bundled into one vfmadd213sd instruction, and suppleness concerning the associativity of 2*x*y.

Sheep race optimization

Anybody that has spent any time round sheep will know that they often attempt to cluster collectively and transfer away from individuals. Say that you’ve got a race of sheep that you just wish to transfer forwards.

The apparent factor to attempt (possibility a within the diagram) is to face on the again within the hope that they are going to all transfer away from you. That doesn’t work, because the again sheep are blocked from transferring by the sheep in entrance of them (they might bunch up a bit of), and the entrance sheep don’t wish to transfer as you’re too distant to be a menace, and so they don’t wish to separate from the remainder of the sheep behind them.

What does work (possibility b within the diagram) is to begin in entrance of all of the sheep and stroll backwards alongside the race. As you stroll previous the entrance sheep, they are going to transfer forwards to get away from you, adopted by the subsequent, then the subsequent, till you’re on the again and they’re all transferring forwards.

This idea of transferring backwards to unstick a pipeline motivates the subsequent optimization – reversing the order of directions to stretch the gap between dependent directions, in order that the pipeline stays full.

a portion of the Instruction tables: Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs for Haswell (different generations are related sufficient that the identical reasoning nonetheless applies):

Instruction	Ports	Latency	Reciprocal Throughput
vaddsd/vsubsd	p1	3	1
vucomisd	p1	3	1
vmulsd	p01	5	0.5
vfmadd213sd	p01	5	0.5

and searching on the meeting above, there are plenty of stalls happening, for example %xmm5 is generated on line 1, however gained’t be prepared for line 2 till 5 clock cycles later, and %xmm4 is generated on line 7 and used on line 9. The processor internally will be capable to mitigate this a bit of by out of order execution.

The loop construction (utilizing <- to indicate “will depend on”) is one thing like:

whereas(D) {
	A <- B
	B <- A
	C <- A,B
	D <- C
}

the place A and B signify the calculation of the subsequent x and y, C and D signify the calculation and comparability of mag_squared to 4.

If calculations of C and D are reordered in order that the hole between them and what they rely upon is so long as attainable, it adjustments to

whereas(D) {
	D <- C
	C <- A,B
	A <- B
	B <- A
}

D now has almost a full loop between being calculated and being wanted (though shall be one loop behind), as does C.

I name this optimization a ‘sheep race’ – there’s nearly definitely a extra widespread identify for that, and I’d welcome being instructed what it’s.

Turning this into code (by the way, rely has been modified to lower fairly than enhance – this ends in no change in efficiency, however will assist later):

iterations_t mandelbrot_sheeprace(double cx, double cy, iterations_t m) {
	double x = 0, y = 0, x_squared = 0, y_squared = 0, mag_squared = 0;
	iterations_t rely = m + 2;
	whereas(rely != 0 && mag_squared <= 4.0) {
		mag_squared = x_squared + y_squared;
		y_squared = y * y;
		x_squared = x * x;
		double newx = x_squared - y_squared + cx;
		y = 2 * y * x + cy;
		x = newx;
		count--;
	}
	return m - rely;
}

This now has the comparability of the mag_squared to 4 happen earlier than the calculation of the subsequent worth of mag_squared, and the calculation of of mag_squared to happen earlier than the calculation of the subsequent values of x_squared and y_squared. Consequently, the worth used into comparability is now two loops behind, and rely must be adjusted accordingly (therefore the m + 2). The additional two x,y values generated are innocent, and now the dependencies are additional aside.

	Check A	Check B	Check C	Check D
time (seconds)	674	47.7	28.7	2462.3
whole iterations	273743913640	19294157420	11543469820	1000000000000

This can be a substantial enchancment, and the one value is that the code is much less apparent.

Interleaving

One other factor to do to fill within the pipeline bubbles is to do the calculations for 2 factors on the identical time. I’ll name every calculation a stream, so there are two streams, (s_1) and (s_2). One loop handles the case of each streams (s_1) and (s_2) working, and the opposite loop handles the case the place one of many streams has completed and solely (s_1) is working. If (s_1) finishes earlier than (s_2), the (s_1) consequence will get returned, and the contents of (s_2) get transferred to (s_1). The pipeline terminates when all streams are completed.

void mandelbrot2_sheeprace(double s1_cx, double s1_cy, double s2_cx, double s2_cy, uint32_t m, 
			iterations_t* s1_r, iterations_t* s2_r) {
	double s1_x = 0, s1_y = 0, s1_x_squared = 0, s1_y_squared = 0, s1_mag_squared = 0;
	double s2_x = 0, s2_y = 0, s2_x_squared = 0, s2_y_squared = 0, s2_mag_squared = 0;
	iterations_t rely = m + 2;
	whereas(rely != 0) {
		if(s2_mag_squared > 4.0) {
			// write out stream 1
			*s2_r = m - rely;
			// now solely stream 0 left
			goto single_point;
		}
		if(s1_mag_squared > 4.0) {
			// write out stream 0
			*s1_r = m - rely;
			// switch stream 1 to stream 0
			s1_r = s2_r;
			s1_x = s2_x;
			s1_y = s2_y;
			s1_cx = s2_cx;
			s1_cy = s2_cy;
			s1_x_squared = s2_x_squared;
			s1_y_squared = s2_y_squared;
			s1_mag_squared = s2_mag_squared;
			// now solely stream 0 left
			goto single_point;
		}
		s1_mag_squared = s1_x_squared + s1_y_squared;
		s2_mag_squared = s2_x_squared + s2_y_squared;
		count--;
		s1_y_squared = s1_y * s1_y;
		s2_y_squared = s2_y * s2_y;
		s1_x_squared = s1_x * s1_x;
		s2_x_squared = s2_x * s2_x;
		s1_y = 2 * s1_y * s1_x + s1_cy;
		s2_y = 2 * s2_y * s2_x + s2_cy;
		s1_x = s1_x_squared - s1_y_squared + s1_cx;
		s2_x = s2_x_squared - s2_y_squared + s2_cx;
	}
	*s1_r = m;
	*s2_r = m;
	return;
single_point:
	whereas(rely != 0 && s1_mag_squared <= 4.0) {
		s1_mag_squared = s1_x_squared + s1_y_squared;
		count--;
		s1_y_squared = s1_y * s1_y;
		s1_x_squared = s1_x * s1_x;
		s1_y = 2 * s1_y * s1_x + s1_cy;
		s1_x = s1_x_squared - s1_y_squared + s1_cx;
	}
	*s1_r = m - rely;
}

The corresponding outer loop is an easy extension of mandelbrot_render_simple – factors are simply processed two at a time:

template<auto F> void mandelbrot_render_simple2(iterations_t* p, double cx, double cy,
			double zoom, uint32_t width, uint32_t top, uint32_t iterations) {
	double xs = cx-0.5 / zoom;
	double ys = cy + 0.5 * top / (zoom * width);
	double inc = 1.0 / (zoom * width);
	for(uint32_t j = 0; j < top; j++) {
		iterations_t* py=p + j * width;
		double y=ys - inc * j;
		for(uint32_t i=0; i < width; i += 2) {
			F(xs + inc * i, y, xs + inc * (i + 1), y, iterations,
					&py[i], &py[i + 1]);
		}
	}
}

This may be expanded to work with (n) streams (s_1 cdots s_n). Examine the github repository for examples with 3 or 4 streams.

Here’s a diagram, displaying a pipeline working with 4 streams. Be aware that stream (s_4) transfers to stream (s_2) then (s_1) as they end.

Outcomes with pipeline of width 2

	#streams	Check A	Check B	Check C	Check D
time (seconds)	2	433.7	30.6	18.2	1562.7
whole iterations	2	273729279660	19298214380	11543330380	1000000000000
time (seconds)	3*	393.7	28.0	16.7	1424.2
whole iterations	3*	273773335620	19336345300	11565894240	1020000000000
time (seconds)	4	393.2	27.9	16.7	1424.4
whole iterations	4	273729279660	19298214380	11543330380	1000000000000

* – There are edge results the place the variety of streams doesn’t evenly divide the row width, leading to barely overflowing the rows on this case (I made the array p a bit of bigger so this might not lead to undefined habits). This is able to be comparatively simple to take care of, however the method taken within the subsequent part gained’t have this situation.

This can be a substantial velocity enchancment – the velocity enhance is now ~2.5x for a pipeline with 3 or 4 streams, however there are diminishing returns for growing the variety of streams. Potential causes for this:

• As extra streams are added, the calculations are interleaved, the time between dependent calculations will get to be larger than or equal to the latency, and which level including streams contributes nothing.

• Extra streams requires extra lively variables. Ideally the internal loop variables ought to all happen in CPU registers – if there are extra variables wanted at any given second than there are registers, the compiler will spill these values to reminiscence, which is slower. The x86_64 household has 16 AVX registers, that are those used for floating level (AVX-512 will increase that quantity to 32). PowerPC, ARM and RISC-V have 32 of every register sort, so might probably deal with a bigger variety of streams.

• The way in which that this code is ready up requires an additional loop for every attainable variety of lively streams (besides zero, wherein case the operate returns), so if there are (n) streams, there are (n) loops. This ends in the code dimension being (O(n^2)). Elevated code dimension would possibly put strain on the execution and department caches on the CPU – though the code must get fairly massive for this to occur on fashionable x86-64 CPUs.

• The operate doesn’t return till all of the streams are executed, so if one stream takes for much longer than the others, the opposite streams gained’t be processing (the black space within the determine above), which signifies that sources are sitting idle.

Interleaving, take two

The following step is to make use of SIMD to arrange extra streams, however the above method gained’t scale to this – the circumstances within the code would get large.

The issue might be restructured be inverting the order of the loops and placing what was the internal loop (mandelbrot_point) on the surface, and turning what was the outer loop (mandelbrot_render_simple) right into a state engine to get the subsequent level to course of:

First the state engine, which merely iterates by the factors (there’s a wrinkle with the tail which shall be described beneath, however might be assumed to be zero and ignored for now):

class MandelbrotStateEngineSimple {
public:
	MandelbrotStateEngineSimple(iterations_t* pp, int pwidth, int pheight,
				double pxs, double pys, double pinc, int ptail) {
		p = pp;
		width = pwidth;
		top = pheight;
		inc = pinc;
		xs = pxs;
		ys = pys;
		w = 0;
		h = 0;
		tail = ptail;
	};
	// get_next level units px, py, x, y and returns true if there
	// is a subsequent level, returns false in any other case
	bool get_next_point(double& px, double& py, int&x, int& y) {
		if(h == top) {
			if(tail == 0) {
				return false;
			}
			tail--;
			px = py = 0;
			// dummy level
			x = 0;
			y = top;
			return true;
		}
		x = w;
		y = h;
		px = xs + w * inc;
		py = ys - h * inc;
		w++;
		if(w >= width) {
			w = 0;
			h++;
		};
		return true;
	}
personal:
	iterations_t* p;
	int width;
	int top;
	int w;
	int h;
	double xs;
	double ys;
	double inc;
	int tail;
};

Right here is an easy render loop that makes use of this state engine – it first calls get_next_point to get the primary level, then each time some extent is completed, it writes it and will get the subsequent level.

void mandelbrot_render1(iterations_t* p, double cen_x, double cen_y, 
		double zoom, int width, int top, iterations_t iterations) {
	double xs = cen_x - 0.5 / zoom;
	double ys = cen_y + 0.5 * top / (zoom * width);
	double inc = 1.0 / (zoom * width);
	MandelbrotStateEngineSimple mse(p, width, top, xs, ys, inc, 0);
	double x = 0, y = 0, x_squared = 0, y_squared = 0;
	double mag_squared = 0;
	iterations_t rely = iterations;
	int dx = 0, dy = 0;
	double cx = 0, cy = 0;
	mse.get_next_point(cx, cy, dx, dy);
	whereas(true) {
		if(rely ==  0 || mag_squared >  4) {
			p[dx + dy * width] = iterations - rely;
			if(!mse.get_next_point(cx, cy, dx, dy)) {
				return;
			}
			// reset the iterators
			x = y = 0;
			rely = iterations;
		}
		count--;
		y = 2 * y * x + cy;
		x = x_squared - y_squared + cx;
		x_squared = x * x;
		y_squared = y * y;
		mag_squared = x_squared + y_squared;
	}
}

The extension of this to a number of streams is simple – every stream is primed with a name to get_next_point, and every time some extent in a stream is completed, it’s written out and changed by the subsequent level, which different streams preserve going.

void mandelbrot_render2_sheeprace(iterations_t* p, double cen_x, double cen_y, 
		double zoom, int width, int top, iterations_t iterations) {
	double xs = cen_x - 0.5 / zoom;
	double ys = cen_y + 0.5 * top / (zoom * width);
	double inc = 1.0 / (zoom * width);
	MandelbrotStateEngineSimple mse(p, width, top, xs, ys, inc, 1);
	double s1_x = 0, s1_y = 0, s1_x_squared = 0, s1_y_squared = 0, s1_mag_squared = 0;
	double s2_x = 0, s2_y = 0, s2_x_squared = 0, s2_y_squared = 0, s2_mag_squared = 0;
	iterations_t s1_count = iterations + 2;
	iterations_t s2_count = iterations + 2;
	int s1_dx = 0, s1_dy = 0;
	int s2_dx = 0, s2_dy = 0;
	double s1_cx = 0, s1_cy = 0;
	double s2_cx = 0, s2_cy = 0;
	mse.get_next_point(s1_cx, s1_cy, s1_dx, s1_dy);
	mse.get_next_point(s2_cx, s2_cy, s2_dx, s2_dy);
	whereas(true) {
		if(s1_count ==  0 || s1_mag_squared >  4) {
			p[s1_dx + s1_dy * width] = iterations - s1_count;
			if(!mse.get_next_point(s1_cx, s1_cy, s1_dx, s1_dy)) {
				return;
			}
			s1_x = s1_y = s1_x_squared = s1_y_squared = s1_mag_squared = 0;
			s1_count = iterations + 2;
		}
		if(s2_count ==  0 || s2_mag_squared >  4) {
			p[s2_dx + s2_dy * width] = iterations - s2_count;
			if(!mse.get_next_point(s2_cx, s2_cy, s2_dx, s2_dy)) {
				return;
			}
			s2_x = s2_y = s2_x_squared = s2_y_squared = s2_mag_squared = 0;
			s2_count = iterations + 2;
		}
		s1_mag_squared = s1_x_squared + s1_y_squared;
		s2_mag_squared = s2_x_squared + s2_y_squared;
		s1_count--;
		s2_count--;
		s1_y_squared = s1_y * s1_y;
		s2_y_squared = s2_y * s2_y;
		s1_x_squared = s1_x * s1_x;
		s2_x_squared = s2_x * s2_x;
		s1_y = 2 * s1_y * s1_x + s1_cy;
		s2_y = 2 * s2_y * s2_x + s2_cy;
		s1_x = s1_x_squared - s1_y_squared + s1_cx;
		s2_x = s2_x_squared - s2_y_squared + s2_cx;
	}
}

That is what the pipeline appears to be like like:

Be aware that the order of entry of factors to the pipeline will not be the identical because the order of exit. For example, a enters the pipeline earlier than b, however b finishes earlier than a.

There must be some care taken when there aren’t any extra factors to course of – if the operate exits instantly, there could also be different factors left within the pipeline that don’t get processed. The answer to that is to flush the pipeline on the finish by processing some additional factors. First the array p is prolonged by one to be of dimension width*top+1, in order that the final index width*top similar to (0,top) can be utilized to dump the flush values. If if there are (n) streams, then (n+1) factors that take the utmost variety of iterations ((c_x=0), (c_y=0) works for this ) is enough to flush the pipeline. The tail parameter is used to insert these pipeline flushing factors.

	Check A	Check B	Check C	Check D
render1	934.5	66.1	39.6	3406.9
render1_sheeprace	855.5	60.6	36.4	3124.9
render2	467.9	33.3	20.1	1706.1
render2_sheeprace	428.5	30.4	18.3	1567.9
render3	494.6	35.1	21.0	1807.4
render3_sheeprace	520.2	39.9	22.1	1899.6
render4	527.4	37.4	22.4	1923.0
render4_sheeprace	560.8	40.1	24.1	2044.3

One shocking factor from these outcomes is that the candy spot appears to be 2 streams, whereas the earlier pipeline had a candy spot of three streams. Wanting on the assembler, get_next_point is getting inlined, which provides to the register strain, which results in extra register spills with 3 streams. I speculate that the slowdown is sufficient to make 3 streams much less aggressive than two.

Vectorizing

Now there’s sufficient of a framework to make use of all of the components of the AVX vector, not simply the one lane that the earlier code has been utilizing. AVX vectors are 256 bits large, which permits SIMD on 4 doubles on the identical time. I’m going to borrow from ARM terminology and name every of those 4 parts a lane.

There are circumstances which the compiler can auto-vectorize code, however I’ve not explored whether or not it’s attainable on this case, so I exploit the Intel Intrinsics right here.

The check for whether or not a stream has completed wants some effort to make quick – it could be costly to check every element of the vector one after the other, so I need one check for every vector. One among my favourite intrinsics is _mm256_movemask_pd which makes use of the vmovmskpd instruction. _mm256_movemask_pd takes the highest bit of every 64 bit element of a vector and packs them into the low 4 bits of an integer. The aim is to get the expression (rely!=0 && mag_squared=<4.0) to set the excessive bit to 1 if it succeeds, 0 in any other case.

Firstly, I modified the rely loop a bit of by decreasing the worth by one, so as a substitute of checking rely!=0 the verify is modified to rely>=0, which can set the excessive bit when it’s time to exit. The _mm256_cmp_pd intrinsic will handle the mag_squared=<4.0. Then or-ing these collectively will make the check prepared for _mm256_movemask_pd.

As soon as the top of a stream is detected, the consequence must be written and get_next_point invoked. Each x86-64 AVX2 and ARM NEON help extracting a price from one lane of a vector, however neither of them help the lane chosen to be a variable, so my answer is to retailer all of the vectors to reminiscence, do the setup for the subsequent level from there, then load the reminiscence again into the registers. That is costly, however for deep zooms will occur occasionally in comparison with the case the place a traditional iteration occurs.

ARM notes: ARM NEON solely has 128 bit vectors so there are solely two lanes per vector, the intrinsics documentation is here, and there’s no direct vmovmskpd instruction on ARM, however the equal might be carried out with a couple of directions. Other than that, the ARM NEON code can have fairly related steps.

void render_avx_sheeprace2(iterations_t* p, double cen_x, double cen_y, 
		double zoom, int width, int top, iterations_t iterations) {
	double xs = cen_x - 0.5 / zoom;
	double ys = cen_y + 0.5 * top / (zoom * width);
	double inc = 1.0 / (zoom * width);
	MandelbrotStateEngineSimple mse(p, width, top, xs, ys, inc, 7);
	iterations--;
	__m256d s1_x = _mm256_setzero_pd();
	__m256d s1_y = _mm256_setzero_pd();
	__m256d s1_x_squared = _mm256_setzero_pd();
	__m256d s1_y_squared = _mm256_setzero_pd();
	__m256d s1_mag_squared = _mm256_setzero_pd();
	__m256i s1_count = _mm256_set1_epi64x(iterations + 2);
	__m256d s2_x = _mm256_setzero_pd();
	__m256d s2_y = _mm256_setzero_pd();
	__m256d s2_x_squared = _mm256_setzero_pd();
	__m256d s2_y_squared = _mm256_setzero_pd();
	__m256d s2_mag_squared = _mm256_setzero_pd();
	__m256i s2_count = _mm256_set1_epi64x(iterations + 2);
	__m256i one_int64 = _mm256_set1_epi64x(1);
	__m256d 4 = _mm256_set1_pd(4.0);
	int s1_dx[4] = {}, s1_dy[4] = {};
	int s2_dx[4] = {}, s2_dy[4] = {};
	double cx_mem[4] __attribute__((__aligned__(32)));
	double cy_mem[4] __attribute__((__aligned__(32)));
	mse.get_next_point(cx_mem[0], cy_mem[0], s1_dx[0], s1_dy[0]);
	mse.get_next_point(cx_mem[1], cy_mem[1], s1_dx[1], s1_dy[1]);
	mse.get_next_point(cx_mem[2], cy_mem[2], s1_dx[2], s1_dy[2]);
	mse.get_next_point(cx_mem[3], cy_mem[3], s1_dx[3], s1_dy[3]);
	__m256d s1_cx = *(__m256d*)cx_mem;
	__m256d s1_cy = *(__m256d*)cy_mem;
	mse.get_next_point(cx_mem[0], cy_mem[0], s2_dx[0], s2_dy[0]);
	mse.get_next_point(cx_mem[1], cy_mem[1], s2_dx[1], s2_dy[1]);
	mse.get_next_point(cx_mem[2], cy_mem[2], s2_dx[2], s2_dy[2]);
	mse.get_next_point(cx_mem[3], cy_mem[3], s2_dx[3], s2_dy[3]);
	__m256d s2_cx = *(__m256d*)cx_mem;
	__m256d s2_cy = *(__m256d*)cy_mem;
	whereas(true) {
		__m256d cmp2;
		cmp2 = _mm256_or_pd((__m256d)s1_count, _mm256_cmp_pd(s1_mag_squared, 4, _CMP_GT_OS));
		if(!_mm256_testz_pd(cmp2, cmp2)) {
			int masks = _mm256_movemask_pd(cmp2);
			_mm256_store_pd(cx_mem, s1_cx);
			_mm256_store_pd(cy_mem, s1_cy);
			double s1_x_mem[4] __attribute__((__aligned__(32)));
			double s1_y_mem[4] __attribute__((__aligned__(32)));
			uint64_t s1_count_mem[4] __attribute__((__aligned__(32)));
			double s1_x_squared_mem[4] __attribute__((__aligned__(32)));
			double s1_y_squared_mem[4] __attribute__((__aligned__(32)));
			_mm256_store_pd(s1_x_mem, s1_x);
			_mm256_store_pd(s1_y_mem, s1_y);
			_mm256_store_si256((__m256i*)s1_count_mem, s1_count);
			_mm256_store_pd(s1_x_squared_mem, s1_x_squared);
			_mm256_store_pd(s1_y_squared_mem, s1_y_squared);
			whereas(masks != 0) {
				int b = __builtin_ctz(masks);
				p[s1_dx[b] + s1_dy[b] * width] = iterations - s1_count_mem[b];
				if(!mse.get_next_point(cx_mem[b], cy_mem[b], s1_dx[b], s1_dy[b])) {
					return;
				}
				s1_count_mem[b] = iterations + 2;
				s1_x_mem[b] = s1_y_mem[b] = s1_x_squared_mem[b] = s1_y_squared_mem[b] = 0;
				masks = masks&~(1U<<b);
			}
			s1_x = _mm256_load_pd(s1_x_mem);
			s1_y = _mm256_load_pd(s1_y_mem);
			s1_count = _mm256_load_si256((__m256i*)s1_count_mem);
			s1_x_squared = _mm256_load_pd(s1_x_squared_mem);
			s1_y_squared = _mm256_load_pd(s1_y_squared_mem);
			s1_cx = _mm256_load_pd(cx_mem);
			s1_cy = _mm256_load_pd(cy_mem);
		}
		cmp2 = _mm256_or_pd((__m256d)s2_count, _mm256_cmp_pd(s2_mag_squared, 4, _CMP_GT_OS));
		if(!_mm256_testz_pd(cmp2, cmp2)) {
			int masks = _mm256_movemask_pd(cmp2);
			_mm256_store_pd(cx_mem, s2_cx);
			_mm256_store_pd(cy_mem, s2_cy);
			double s2_x_mem[4] __attribute__((__aligned__(32)));
			double s2_y_mem[4] __attribute__((__aligned__(32)));
			uint64_t s2_count_mem[4] __attribute__((__aligned__(32)));
			double s2_x_squared_mem[4] __attribute__((__aligned__(32)));
			double s2_y_squared_mem[4] __attribute__((__aligned__(32)));
			_mm256_store_pd(s2_x_mem, s2_x);
			_mm256_store_pd(s2_y_mem, s2_y);
			_mm256_store_si256((__m256i*)s2_count_mem, s2_count);
			_mm256_store_pd(s2_x_squared_mem, s2_x_squared);
			_mm256_store_pd(s2_y_squared_mem, s2_y_squared);
			whereas(masks != 0) {
				int b = __builtin_ctz(masks);
				p[s2_dx[b] + s2_dy[b] * width] = iterations - s2_count_mem[b];
				if(!mse.get_next_point(cx_mem[b], cy_mem[b], s2_dx[b], s2_dy[b])) {
					return;
				}
				s2_count_mem[b] = iterations + 2;
				s2_x_mem[b] = s2_y_mem[b] = s2_x_squared_mem[b] = s2_y_squared_mem[b] = 0;
				masks = masks&~(1U<<b);
			}
			s2_x = _mm256_load_pd(s2_x_mem);
			s2_y = _mm256_load_pd(s2_y_mem);
			s2_count = _mm256_load_si256((__m256i*)s2_count_mem);
			s2_x_squared = _mm256_load_pd(s2_x_squared_mem);
			s2_y_squared = _mm256_load_pd(s2_y_squared_mem);
			s2_cx = _mm256_load_pd(cx_mem);
			s2_cy = _mm256_load_pd(cy_mem);
		}
		// s1_mag_squared = s1_x_squared + s1_y_squared;
		s1_mag_squared = _mm256_add_pd(s1_x_squared, s1_y_squared);
		s2_mag_squared = _mm256_add_pd(s2_x_squared, s2_y_squared);
		// Decrement counter
		s1_count = _mm256_sub_epi64(s1_count, one_int64);
		s2_count = _mm256_sub_epi64(s2_count, one_int64);
		// s1_y_squared = s1_y * s1_y;
		s1_y_squared = _mm256_mul_pd(s1_y, s1_y);
		s2_y_squared = _mm256_mul_pd(s2_y, s2_y);
		// s1_x_squared = s1_x * s1_x;
		s1_x_squared = _mm256_mul_pd(s1_x, s1_x);
		s2_x_squared = _mm256_mul_pd(s2_x, s2_x);
		// s1_y = 2 * s1_y * s1_x + s1_cy;
		s1_y = _mm256_fmadd_pd(_mm256_add_pd(s1_x, s1_x), s1_y, s1_cy);
		s2_y = _mm256_fmadd_pd(_mm256_add_pd(s2_x, s2_x), s2_y, s2_cy);
		// s1_x = s1_x_squared - s1_y_squared + s1_cx;
		s1_x = _mm256_sub_pd(s1_x_squared, _mm256_sub_pd(s1_y_squared, s1_cx));
		s2_x = _mm256_sub_pd(s2_x_squared, _mm256_sub_pd(s2_y_squared, s2_cx));
	}
}

	Check A	Check B	Check C	Check D
render1_avx	215.6	15.6	9.5	786.7
render1_avx_sheeprace	207.9	15.1	9.2	750.0
render2_avx	126.3	9.3	5.8	458.3
render2_avx_sheeprace	119.7	8.8	5.4	436.8
render3_avx	131.6	9.7	6.0	479.2
render3_avx_sheeprace	131.2	9.8	6.0	480.6

render2_avx_sheeprace is the clear winner right here, which isn’t an important shock given the sooner outcomes.

Conclusion

render2_avx_sheeprace with check D is working ~2.3 billion iterations/second on a single core. On a 2.9 GHz processor, that’s a mean of 10.13 clock cycles/iteration for every of 8 streams.

Evaluating this with the naive implementation proper in the beginning:

	Check A	Check B	Check C	Check D
naive	973.4	68.7	41.2	3550.7
render2_avx_sheeprace	119.7	8.8	5.4	436.8
velocity multiplier	8.1	7.8	7.6	8.1

So utilizing a trying fastidiously on the drawback, realizing concerning the machine structure and utilizing number of optimization strategies, I’ve managed to get a velocity enchancment of about ~8x.

Here’s a tough breakdown of of how a lot every approach contributes to the most effective consequence:

	multiplier
-ffast-math	1.04
Sheep race	1.09
Interleaving (2)	2.0
Vectorization	3.6

Other than vectorization and interleaving utilizing the identical framework with pipelining and a state engine, these are largely unbiased of one another.

I’m positive that this isn’t the ultimate enchancment that may very well be made on this space – I welcome some other recommendations, and please be happy to contact me.

There shall be a component 2 to this – in that, I’ll focus on alternate options to the MandelbrotStateEngineSimple state engine that may prune the variety of factors to calculate, and make some observations about multithreading. It will lead to one other substantial velocity enchancment on prime of what I’ve already offered (attending to ~100x in some circumstances) utilizing these strategies.

Any feedback, corrections, observations, or optimizations that I missed? Please be happy to contact me