In 2012, I made a Mandelbrot renderer in JavaScript. Here I will explain how you can do it yourself, including how to implement smooth coloring and making it fast.

The famous Mandelbrot set is
a set of points in the complex plane. In essence, what we want to find out
is if the iterative function C below will *converge* to some constant or *diverge*
to infinity.

The definition is simply

z

_{n+1}= z_{n}^{2}+ c

with the initial condition simply formed by taking the coordinates in the complex plane,

z

_{0}= c = x + iy

Pretend for a moment that you've never seen the Mandelbrot plot before. By only looking at the definition above, can you guess how it would look like when rendered? First we need to look at a few expansions of z expressed with x and y.

z

_{0}= x + iyz

_{1}= (x + iy)^{2}+ (x + iy) = (x^{2}- y^{2}+ x) + i(2xy + y)z

_{2}= z_{1}^{2}+ (x + iy) = ...

So, as n → ∞, for which values of x and
y will z_{n} converge, and for which will it
diverge?

Well, for large values of |x| and small for |y|
(points close to origo), the sequence might diverge, because
x^{2} would be dominating over y^{2}.
Diverged points are painted black, so we can guess that the plot will be black
in all directions some distance from origo.

But for points close to origo — say, |x| and
|y| less than 1 — we would expect it to converge, because
the product of two numbers less than one will always be smaller than each of
its parts, giving small values for x^{2} - y^{2}. So in a
circle around origo with a radius of one, we'd expect to see colored pixels.

But what if the signs of x and y differ? Then things quickly get more complicated. In fact, if we plot using a computer, we won't get a nice colored disc centered at origo. What we get is an infinitely complex and fractured plot.

We can even zoom endlessly into the plot and it will *still* be as
non-uniform and complex as before. While a disc would have a dimension of two,
a fractal has a so-called so-called Hausdorff dimension
which is something in between a line and a plane. We'd get a non-integer
dimension; a *fractal*.

Calculating the Mandelbrot set numerically is easy.

Given the equations above, take any point z_{0} = (x, y)
and then calculate z_{1} = (x+iy)^{2} + (x+iy) and
continue doing this. For practical purposes, let's decide on a *threshold*
value. If the magnitude of z — the distance to origo, or
√(x^2+y^2)— ever becomes larger than this threshold
value we will assume that it will diverge into infinity. If so, stop the
calculation and plot a *black dot* at the current location.

*
For the Mandelbrot set, it can be shown that the threshold value is exactly two, i.e.
any sequence with a |z _{n}| > 2 will diverge.
Read more about this.
*

If |z| has not exceeded the threshold value after a predecided number of iterations (which we choose at will), we will assume that the current parameters makes the function converge. In this case, plot a non-black dot at the current location.

I said above that if the function diverges, one should plot a non-black dot.
One could simply paint a white dot here. But instead, maybe we want to get
an idea of *how fast* the function is diverging to infinity at this point.

To do this, just take the current value of the number of steps performed
and *map* that against a color spectrum, and paint that color.

So, functions diverging quickly will get about the same color.

If you use the number of iterations to pick a color, you'll get ugly color bands in the plot. There is a really cool trick to get smooth, gradual color changes.

So, you *basically* calculate `Z = Z^2`

until it diverges and make a note of
the iteration count. What we really want, though, is a *fractional*
iteration count, so we can multiply that with a color value to get smooth
colors.

The trick is to note that when you calculate `Z = Z^2`

you'll get values ```
Z,
Z^2, Z^4, Z^8
```

and so on. If you take the logarithm of this, you'll get the
values 1, 2, 4, 8 etc. If you take the logarithm one more time, you'll get
1, 2, 3, 4, 5 etc. So to get a fractional number of iterations, just do:

This is all explained over at linas.org.

In my code, I originally used the following smoothing equation:

With some elementary logarithm rules, we can simplify this to

which is faster. The constant `5`

is another little trick, which should
be explained in the code itself.

Finally, when you calculate the color value of a single pixel, it is in
reality just the color of a single point in the Mandelbrot set that is
situated somewhere *inside* that pixel.

What I'm saying is that you'll basically get pixel artifacts in the image, especially in dense areas where the color changes (near the black set, for instance).

So what I do is to use random sampling: Just sample a given number of random points inside the pixel and average the sum of the color values. This is equivalent to rendering the plot at a higher resolution and scaling down.

There are many supersampling techniques to use, and the random sampling was chosen because of its simplicity. The problem is that the resulting picture will look a bit blurry (there are ways around this as well).

Calculating the Mandelbrot set is quite slow, but there are a lot of tricks to speed it up.

When speeding up any code, the first step (after making the code *correct*,
of course) is to look at the algorithm and try to use one with a simpler
complexity class. Unfortunately, for the Mandelbrot set, this isn't really
possible. So the tricks mentioned here are all cases of
*micro-optimizations*. Nevertheless, they will improve the running time
quite a lot.

We also have to remember that we're using Javascript here, which is a relatively slow language because of its dynamic nature. What's interesting in this regard, though, is to identify performance hot spots in the typical javascript engines. It's interesting to test the code on different browsers.

First, let's look at the inner loop. It continually calculates the magnitude of the complex number C, and compares this with a threshold value. Observe that it takes the square root in doing so:

If we just square the treshold value, we should be able to do away with the square root operation:

You've probably noticed that the plot is reflected vertically over the line y = 0. You can take advantage of this mirroring to halve the computation time. I don't, because you'll mostly render plots that are massively zoomed in.

The main equation is

z

_{n+1}= z_{n}^{2}+ c

Setting `C = z`

and `Cr = Re(z)`

and `Ci = Im(z)`

, we get

giving us

If we introduce two variables `Tr = Cr^2`

and `Ti = Ci^2`

, we get

So we have now replaced some multiplications with additions, which is
normally faster in most CPUs. But, again, this is javascript, and
javascript has quite a different performance profile. The code above indeed
does *not* give us any **significant** speedup --- for a 640x480 image, we
only save a hundred milliseconds, or so.

To plot individual pixels in HTML5 canvas, you get an array and you have to calculate the array offset for a given coordinate pair.

I.e., given RGBA pixel format (four positions), an (x, y) coordinate pair and a width and height, you calculate it by

so that you can now set the RGBA values as

There are several ways of optimizing this. For instance, we can simply
multiply the whole offset by four, which is the same as shifting all bits
left two positions. However, javascript works in mysterious ways, so the
customary shift operations may not be as fast as in other languages like C
and C++. The reason *probably* has to do with the fact that javascript only
has *one* data type for numbers, and my guess is that it's some kind of
float.

Anyway, we now have

Another trick I'd like to mention. Say that the width and height are fixed to, say 640 and 480, respectively. And old trick to multiply y by 640 would be notice that 640 = 512 + 128 = 2^9 + 2^7, giving us

So now we've converted one multiplication into two shifts and an add. In your commodity language and hardware, this might be quite fast in tight innerloops.

Anyway, we still want to be able to use arbitrary heights and widths, so let's skip that one.

By far, the fastest way of accessing the array is by doing it sequentially.

That is, instead of doing

a *much* faster way would be to do

So now we've basically saved the work of doing `2*width*height`

multiplications, or 600 thousand of them, assuming a 640x480 image.

To draw in the canvas, you request an array, update it and copy it back to the canvas.

Of course, you want to reduce the number of such operations. Because we want an animation showing each line as it is drawn, we'll do this:

- Get an image data array
- For each line: Update the array
- For each line: Copy the array back to the canvas

The trick here, though is to *not* use `getImageData`

. You're going to
overwrite all existing image data, so you can use the same buffer for every
line. So instead, we'll use these operations:

- Get a line buffer by calling
`createImageData(canvas.width, 1)`

- For each line: Update the line buffer array
- For each line: Call
`putImageData(linebuffer, 0, y_position)`

to copy only*one*line

This ensures that we only copy *one* line per frame update.

Since the algorithm above is referentially transparent, meaning that it
always produces the same output for the same input (where input is defined
as `x, y, steps, threshold`

), you could in theory calculate all points in
parallel.

Such algorithms are colloquially called embarrassingly parallel.

Now, JavaScript is inherently single-threaded: You can only use so-called green threads, meaning that the javascript engine will yield control between them.

However, there is a new HTML5 APi called web workers that you can use to create real, OS-level threads. That should make it easy to split up plotting into several threads (preferrably one per logical core).

The algorithm is very well suited for vector operations. Most modern computers come with hardware optimizations for such operations (SSE or the GPU, etc.) However, we are again limited to what the javascript engines will do for us.

Take a look at the optimizations done to the Mandelbrot set in The Computer Language Benchmarks Game

There are a lot of cool tricks going on there. Most of *those* use SSE
parallelism for hardware speedup or offloads to the GPU.