How to Write a Mandelbrot Viewer

The Mandelbrot equation is easy to write, but there are some interesting challenges in programming a graphical viewer. Mathematically, the set is defined by this sequence:

Z = Z² + C    where  C = x + iy, the point to test. The sequence is iterated to infinity with an initial value of Z = 0. If the magnitude of Z remains finite, then the point is in the set. Otherwise, if the magnitude grows infinitely large (diverges), then the point is not in the set.

The naive implementation is always a good place to start. This pseudo-code takes a graph point (x, y) and computes a grayscale color there by measuring how quickly the magnitude diverges when it does. (For the sake of clarity, the complex number arithmetic is written explicitly as real and imaginary parts.)

Var Cr = x Var Ci = y Var Zr = 0 Var Zi = 0 Loop N = 0 to 64 Zr = Zr*Zr - Zi*Zi + Cr Zi = 2.0*Zr*Zi + Ci If sqrt(Zr*Zr + Zi*Zi) < 2 Then Exit Loop Do Loop PixelColor = N*4

To draw an image, you must run this code for each pixel. (The pixel coordinates must be transformed to graph space depending on the viewer's current pan and zoom state. This is a standard problem of any graphical viewer and can be easily solved by using a library, like my open source IceView library.)

That simple Mandelbrot code generates images like this:

Center: (x = -0.377, y = 0.676), Zoom: (1.4 x 10⁵)

It works! This is a lovely image of the Mandelbrot set, which can be colored using any function you want. I use a color function of (Red, Green*1.5, Blue*10) in my Mandelbrot Viewer to make it look like ice crystals!

Speed Optimizations

If you clicked the link to my viewer, you've probably noticed that it works in real time, while the code example takes several seconds to render the image. Speed optimizations and GPU acceleration are the key. There are two simple optimizations that make a significant improvement.

1. Remove the sqrt() by squaring both sides of the comparison. The expression sqrt(Zr*Zr + Zi*Zi) < 2 is equivalent to Zr*Zr + Zi*Zi < 4, but this is much faster!
2. Save the result of Zr*Zr and Zi*Zi from each previous loop, thus removing two more multiplications.

The optimized code looks like this:

Var Cr = x Var Ci = y Var Zr = 0 Var Zi = 0 Var Zr_sq = 0 Var Zi_sq = 0 Loop N = 0 to 64 Zr = Zr_sq - Zi_sq + Cr Zi = 2.0*Zr*Zi + Ci Zr_sq = Zr*Zr Zi_sq = Zi*Zi If Zr_sq + Zi_sq < 4 Then Exit Loop Do Loop PixelColor = N*4

And the image drawn looks identical.

GPU Acceleration

However, the greatest speed improvement comes courtesy of the GPU. Notice that every pixel in the image is calculated independently of the others. They can all be drawn in parallel!

However, if you simply make a pixel shader with the above code, you will have a nasty surprise. The number precision is cut in half (32 bits instead of 64), so you can hardly zoom in at all! Even if you have a fancy GPU with 64-bit support, this is not available in WebGL.

However, there is a solution: Use a pair of singles to approximate double-precision.

My WebGL viewer uses pairs of numbers to enhance the precision to approximately 50 bits (100,000 times more zoom). It's not as good as true 64-bit double precision, but the difficulty in approximating higher precision goes up exponentially, so 50 bits is a compromise that keeps the viewer fast enough to run on most GPUs.

Mathematically, given a number x, you can express x as the sum a + b. There are an infinite number of choices for a and b, but we choose them such that half of x's digits are in each one.

x = 0.101010101111100101010011001 a = 0.101010101111000000000000000 b = 0.000000000000100101010011001

When we pass x to the GPU half of it's digits are lost, but when we pass a and b, then already have half, so nothing is lost! This is the first step, because once information is lost, it's game over.

Since x is equal to a + b, we can use treat a and b as unknowns and carry them through all math operations using algebra. Then we separate the lesser and greater parts to prevent them from being added together.

For example:

x*x = a*a + 2*a*b + b*b a_result = a*a + 2*a*b b_result = b*b

However, after a few operations, the digits begin to drift, leaving a and b either too close or too far to each other. This can be fixed by normalizing them when needed:

sum = a + b // The value of x, but without extra precision. b = b - (sum - a) // Computes the residual error, which becomes the new b. a = sum

Mathematically, the normalization does nothing. However, because the numbers are truncated to their precision after each operation, the result reveals exactly where the significant bits stop and become all zeros, thus enabling us to start b's digits at that point to extend the precision of a.

If you apply these principles to the simple Mandelbrot code to create an extended-precision GPU shader, the GLSL code looks like this:

highp float zrh = 0.0; highp float zrl = 0.0; highp float zih = 0.0; highp float zil = 0.0; highp float zrsq_h = 0.0; highp float zrsq_l = 0.0; highp float zisq_h = 0.0; highp float zisq_l = 0.0; for (int n = 0; n < 2048; n++) { zil = (zrh*zil + zrl*(zih + zil))*2.0 + yl; zih = zrh*zih*2.0 + yh; zrh = zrsq_h - zisq_h + xh; zrl = zrsq_l - zisq_l + xl; if (zrh > 2.0 || zrh < -2.0) { highp float totalr = zrh + zrl; zrl = zrl - (totalr - zrh); zrh = totalr; } if (zih > 2.0 || zih < -2.0) { highp float totali = zih + zil; zil = zil - (totali - zih); zih = totali; } zrsq_h = zrh*zrh; zrsq_l = (zrh*2.0 + zrl)*zrl; zisq_h = zih*zih; zisq_l = (zih*2.0 + zil)*zil; m = zrsq_h + zrsq_l + zisq_h + zisq_l; if (m > 4.0) {en = n; break;} if (n == nmax) break; }

Because each operation loses a little bit of precision, the result is slightly less than 64-bit precision. It gets worse the further you are from (0,0). The precision is usually between 48 and 50 bits of precision, but that is a significant improvement over 32 bits.

You can experiment with a GPU implementation by downloading this HTML file (fractalviewer.html) and opening it in your browser to use it. Then open the file in a text editor, and scroll down until you see the pixel shader, which has the simple Mandelbrot code. Try uncommenting a line to make it Z = Z⁴ + C and see what that does! You can modify the code to draw nearly any iterative 2D fractal, and experiment with the colors too.