There are many programs and algorithms used to plot the Mandelbrot set and other fractals, some of which are described in fractal-generating software. These programs use a variety of algorithms to determine the color of individual pixels efficiently.
Escape time algorithm
The simplest algorithm for generating a representation of the Mandelbrot set is known as the "escape time" algorithm. A repeating calculation is performed for each x, y point in the plot area and based on the behavior of that calculation, a color is chosen for that pixel.
Unoptimized naïve escape time algorithm
In both the unoptimized and optimized escape time algorithms, the x and y locations of each point are used as starting values in a repeating, or iterating calculation (described in detail below). The result of each iteration is used as the starting values for the next. The values are checked during each iteration to see whether they have reached a critical "escape" condition, or "bailout". If that condition is reached, the calculation is stopped, the pixel is drawn, and the next x, y point is examined. For some starting values, escape occurs quickly, after only a small number of iterations. For starting values very close to but not in the set, it may take hundreds or thousands of iterations to escape. For values within the Mandelbrot set, escape will never occur. The programmer or user must choose how many iterations–or how much "depth"–they wish to examine. The higher the maximal number of iterations, the more detail and subtlety emerge in the final image, but the longer time it will take to calculate the fractal image.
Escape conditions can be simple or complex. Because no complex number with a real or imaginary part greater than 2 can be part of the set, a common bailout is to escape when either coefficient exceeds 2. A more computationally complex method that detects escapes sooner, is to compute distance from the origin using the Pythagorean theorem, i.e., to determine the absolute value, or modulus, of the complex number. If this value exceeds 2, or equivalently, when the sum of the squares of the real and imaginary parts exceed 4, the point has reached escape. More computationally intensive rendering variations include the Buddhabrot method, which finds escaping points and plots their iterated coordinates.
The color of each point represents how quickly the values reached the escape point. Often black is used to show values that fail to escape before the iteration limit, and gradually brighter colors are used for points that escape. This gives a visual representation of how many cycles were required before reaching the escape condition.
To render such an image, the region of the complex plane we are considering is subdivided into a certain number of pixels. To color any such pixel, let be the midpoint of that pixel. We now iterate the critical point 0 under , checking at each step whether the orbit point has modulus larger than 2. When this is the case, we know that does not belong to the Mandelbrot set, and we color our pixel according to the number of iterations used to find out. Otherwise, we keep iterating up to a fixed number of steps, after which we decide that our parameter is "probably" in the Mandelbrot set, or at least very close to it, and color the pixel black.
In pseudocode, this algorithm would look as follows. The algorithm does not use complex numbers and manually simulates complex-number operations using two real numbers, for those who do not have a complex data type. The program may be simplified if the programming language includes complex-data-type operations.
for each pixel (Px, Py) on the screen do x0 := scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2.00, 0.47)) y0 := scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-1.12, 1.12)) x := 0.0 y := 0.0 iteration := 0 max_iteration := 1000 while (x*x + y*y ≤ 2*2 AND iteration < max_iteration) do xtemp := x*x - y*y + x0 y := 2*x*y + y0 x := xtemp iteration := iteration + 1
color := palette[iteration] plot(Px, Py, color)
Here, relating the pseudocode to , and :
and so, as can be seen in the pseudocode in the computation of x and y: and
To get colorful images of the set, the assignment of a color to each value of the number of executed iterations can be made using one of a variety of functions (linear, exponential, etc.). One practical way, without slowing down calculations, is to use the number of executed iterations as an entry to a palette initialized at startup. If the color table has, for instance, 500 entries, then the color selection is n mod 500, where n is the number of iterations.
Optimized escape time algorithms The code in the previous section uses an unoptimized inner while loop for clarity. In the unoptimized version, one must perform five multiplications per iteration. To reduce the number of multiplications the following code for the inner while loop may be used instead:
x2:= 0 y2:= 0 w:= 0
while (x2 + y2 ≤ 4 and iteration < max_iteration) do x:= x2 - y2 + x0 y:= w - x2 - y2 + y0 x2:= x * x y2:= y * y w:= (x + y) * (x + y) iteration:= iteration + 1
The above code works via some algebraic simplification of the complex multiplication:
Using the above identity, the number of multiplications can be reduced to three instead of five.
The above inner while loop can be further optimized by expanding w to
Substituting w into yields
and hence calculating w is no longer needed.
The further optimized pseudocode for the above is:
x:= 0 y:= 0 x2:= 0 y2:= 0
while (x2 + y2 ≤ 4 and iteration < max_iteration) do x2:= x * x y2:= y * y y:= 2 * x * y + y0 x:= x2 - y2 + x0 iteration:= iteration + 1
Note that in the above pseudocode, seems to increase the number of multiplications by 1, but since 2 is the multiplier the code can be optimized via .
Coloring algorithms
In addition to plotting the set, a variety of algorithms have been developed to efficiently color the set in an aesthetically pleasing way show structures of the data (scientific visualisation)
Histogram coloring
A more complex coloring method involves using a histogram which pairs each pixel with said pixel's maximum iteration count before escape/bailout. This method will equally distribute colors to the same overall area, and, importantly, is independent of the maximum number of iterations chosen.
This algorithm has four passes. The first pass involves calculating the iteration counts associated with each pixel (but without any pixels being plotted). These are stored in an array IterationCounts[x][y], where x and y are the x and y coordinates of said pixel on the screen respectively.
The top row is a series of plots using the escape time algorithm for 10000, 1000 and 100 maximum iterations per pixel respectively. The bottom row uses the same maximum iteration values but utilizes the histogram coloring method. Notice how little the coloring changes per different maximum iteration counts for the histogram coloring method plots.
The first step of the second pass is to create an array NumIterationsPerPixel[n], where the array size n is the maximum iteration count. Next, one must iterate over the array of pixel-iteration count pairs IterationCounts[x][y], and retrieve each pixel's saved iteration count, i, via e.g. i = IterationCounts[x][y]. After each pixel's iteration count i is retrieved, it is necessary to index the NumIterationsPerPixel array at i and increment the indexed value (which is initially zero) -- e.g. NumIterationsPerPixel[i] = NumIterationsPerPixel[i] + 1.
for (x = 0; x < width; x++) do for (y = 0; y < height; y++) do i:= IterationCounts[x][y] NumIterationsPerPixel[i]++
The third pass iterates through the NumIterationsPerPixel array and adds up all the stored values, saving them in total. The array index represents the number of pixels that reached that iteration count before bailout.
total: = 0 for (i = 0; i < max_iterations; i++) do total += NumIterationsPerPixel[i]
After this, the fourth pass begins and all the values in the IterationCounts array are indexed, and, for each iteration count i, associated with each pixel, the count is added to a global sum of all the iteration counts from 1 to i in the NumIterationsPerPixel array . This value is then normalized by dividing the sum by the total value computed earlier.
hue[][]:= 0.0 for (x = 0; x < width; x++) do for (y = 0; y < height; y++) do iteration:= IterationCounts[x][y] for (i = 0; i <= iteration; i++) do hue[x][y] += NumIterationsPerPixel[i] / total /* Must be floating-point division. */
...
color = palette[hue[m, n]]
...
Finally, the computed value is used, e.g. as an index to a color palette.
This method may be combined with the smooth coloring method below for more aesthetically pleasing images.
Continuous (smooth) coloring
The escape time algorithm is popular for its simplicity. However, it creates bands of color, which, as a type of aliasing, can detract from an image's aesthetic value. This can be improved using an algorithm known as "normalized iteration count", which provides a smooth transition of colors between iterations. The algorithm associates a real number with each value of z by using the connection of the iteration number with the potential function. This function is given by
where zn is the value after n iterations and P is the power for which z is raised to in the Mandelbrot set equation (zn+1 = znP + c, P is generally 2).
If we choose a large bailout radius N (e.g., 10100), we have that
for some real number , and this is
and as n is the first iteration number such that |zn| > N, the number we subtract from n is in the interval [0, 1).
For the coloring we must have a cyclic scale of colors (constructed mathematically, for instance) and containing H colors numbered from 0 to H − 1 (H = 500, for instance). We multiply the real number by a fixed real number determining the density of the colors in the picture, take the integral part of this number modulo H, and use it to look up the corresponding color in the color table.
For example, modifying the above pseudocode and also using the concept of linear interpolation would yield for each pixel (Px, Py) on the screen do x0:= scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2.5, 1)) y0:= scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-1, 1)) x:= 0.0 y:= 0.0 iteration:= 0 max_iteration:= 1000 // Here N = 2^8 is chosen as a reasonable bailout radius.
while x*x + y*y ≤ (1 << 16) and iteration < max_iteration do xtemp:= x*x - y*y + x0 y:= 2*x*y + y0 x:= xtemp iteration:= iteration + 1
// Used to avoid floating point issues with points inside the set. if iteration < max_iteration then // sqrt of inner term removed using log simplification rules. log_zn:= log(x*x + y*y) / 2 nu:= log(log_zn / log(2)) / log(2) // Rearranging the potential function. // Dividing log_zn by log(2) instead of log(N = 1<<8) // because we want the entire palette to range from the // center to radius 2, NOT our bailout radius. iteration:= iteration + 1 - nu
color1:= palette[floor(iteration)] color2:= palette[floor(iteration) + 1] // iteration % 1 = fractional part of iteration. color:= linear_interpolate(color1, color2, iteration % 1) plot(Px, Py, color)
Exponentially mapped and cyclic iterations
Typically when we render a fractal, the range of where colors from a given palette appear along the fractal is static. If we desire to offset the location from the border of the fractal, or adjust their palette to cycle in a specific way, there are a few simple changes we can make when taking the final iteration count before passing it along to choose an item from our palette.