The Mandelbrot Set

What is the Mandelbrot Set?

Let's define the following sequence in the complex plane:

z_n+1  = z_n² + c  with  z₀ = (0,0)

This sequence depends on parameter c. If ||z_n|| remains bounded for all n for a given c, then c belongs to the Mandelbrot set.

There are a few properties useful for computing the Mandelbrot set. First, any given point c such that ||c|| is greater than 2 does not belong to the set. The proof is easy and relies only on the classical property :

 | ||a||-||b|| | <= ||a+b||

which comes directly from the triangular inequality. So we have:

||z_n+1|| = ||z_n²+c|| >= | ||z_n²||-||c||

Let's suppose that ||c||=2+a, with a>0 and that we have the following property for a given n>0

 ||z_n||>2+na

Then we can write:

||z_n+1|| >= | ||z_n²||-||c|| | > (2+na)²-2-a > 2+3na+(na)² > 2+(n+1)a

So, if the property is true for n, it is true for n+1. As the property is true for n=1, it is true for all n, and ||z_n|| is not bound if ||c||>2.

For the same reason, it is easy to show that if, for a given c, k exists such as ||z_k||>2, then c is not in the Mandelbrot set.

So to find points inside the Mandelbrot set, we have only to check for points inside the circle (O,R=2), and we can stop calculus as soon as we find some k such as ||z_k||>2.

Strictly speaking, the Mandelbrot set should be black and white. But soon, the first k value such as ||z_k||>2 was used as a color index. The problem to have good looking Mandelbrot set is two folded. On the one side, you have to design a good colormap. On the other side you have to choose an upper bound for the value to test: the property of belonging to the mandelbrot set is in the general case an indecidable problem. If you find a k such as ||z_k||>2 for a given c then c is not in the set. However, if you can prove for some points that they do belong to the set (i.e there isn't any k such as the upper property is true), it is impossible in the general case, and you can't keep on computing for ever. The upper bound U chosen for the values of k checked gives the level of detail of the set. The higher it is, the more detail you have, but the longer it takes to compute the set, especially when you are close to the border.

When you are zooming inside the set, choosing U is a major concern. However, there are other limitations, such as the numerical precision of your processor.

I started working on Mandelbrot sets in 1985 during a training period at the Ecole Polytechnique under the direction of Jean-François Colonna. I wrote the assembly language routines to compute the set on RIDGE computers, which were one of the first RISC computers available. It is now a hobby and I still take the time to rewrite a few lines of assembly language when new processors or new instructions are available. The set below is computed in assembly language using the EMMX instruction set available on the Pentium-III computers. It is an extension of the Pentium MMX SIMD instructions set to floating point calculations. The EMMX instruction set operates on simple precision (24 bits of mantissa, 8 bits of exponent) float vectors of size 4. The assembly code routine is here .

It is quite interesting to try to estimate the real performance of the processor. On a P-3 800Mhz, a 1280x1024 image with U=256 is computed in 0.77 second. As the assembly routine returns the complement of the number of loops, it is easy to compute the total number of loops executed. Each loop contains 15 vector floating point operations (VFlop) and two integer instructions. 21147490 loops are executed, so 317212350 VFlop are executed. This gives a rate of 0.41 GVFlops, soit 1.6 GFlops. It must be noted that on a P-4 2.53Ghz, the same program takes 0.29 s, which gives 1.1GVflops or 4.4GFlops. On an Athlon 2100+ (real frequency 1.7Ghz), it takes 0.35s, for a performance of 3.6GFlops. (At the same frequency, the Athlon is slightly faster than the P-3 which is itself faster than the P-4; these results are due to differences in the internal structures of the processors).

With the Pentium-4, you can also use 2 elements vectors in double precision, instead of 4 elements in simple precision. It is slower, but you have a much better precision. The assembly code routine is here . On the same P-4 at 2.53Ghz, the program runs 41588101 loops in 0.52s, for a performance of 2.4GFlops.