 Download an exciting AlmondBread technology demo for DOS that gives you a preview of what speed you can expect from the upcoming AlmondBread 0.3 which will include support for Windows95/NT. Here's the README.

AlmondBread is a program for the X Window System that generates images of the Mandelbrot set. Here are some of its features:

• Fast algorithms, among them an implementation of Simultaneous Orbit Iteration.
• Tcl/Tk user interface.
• Save/Load Fractint-compatible GIF files.
• Load Fractint-style colormap files (`.MAP` files).
• Shared or own colormap using an allocation scheme similar to that of xv.
• Easy-to-use zoombox function.

News in version 0.22 include:

• Boundary Tracing.
• Improved Unroll Iteration Method.
• experimental Multithreading support (using POSIX 1003.1c threads).
• several bug fixes (e.g. 16bit Visuals work properly now).

AlmondBread is known to work on Linux using XFree86, Solaris, and IRIX. It requires X11 and the aforementioned Tcl/Tk (at least version 7.4/4.0), plus the jstools-package by Jay Sekora.

View a screenshot of AlmondBread in action, then download AlmondBread (400k gzipped tar-file from Berlin, Germany). Outside of Europe, you might want to use the Canadian mirror instead.
By popular demand, it is also available as a zip file.

The remainder of this page is an adapted version of the AlmondBread help file.

### The Mandelbrot Set

This section describes some of the basic mathematical concepts behind the Mandelbrot set and its generation.

The Mandelbrot set - or M - was discovered in 1980 by mathematician Benoit B. Mandelbrot. It is considered to be the most complex object mathematics has ever seen .

In  M is defined as follows: "The Mandelbrot set is the set of all complex numbers c such that iterating z->z²+c does not go to infinity (starting with z=0)".

It can be shown that the sequence diverges if the magnitude of z exceeds 2.

#### The Algorithm

The basic algorithm to generate an image of M and its basin of attraction for a given region R in the complex plane goes like this: For each pixel p representing a point c in R do the following:
1. z := 0, n := 0.
2. z := z² + c.
3. n := n + 1.
4. If |z|<2 and n<N go to 2.

N is the maximum number of iterations you allow (smaller values mean faster execution, larger ones better accuracy). If the iteration loop terminates because n>=N, there is a good chance c lies in M, so you can color p accordingly. If, however, |z|>=2 then p lies outside of M and you can color p according to the number of iterations when z "bailed out", i.e. according to n.

### Drawing Methods

A drawing method refers to the sequence of pixels that are chosen to generate the image of M for a given region R. The goal, of course, is to compute as few pixels as possible to speed up the imaging process. Certain properties of M allow us not to compute each and every pixel of the image.

#### Scan

This is the simplest method of generating an image of M. As its name suggests, this algorithm simply scans the chosen region from top to bottom, doing a complete iteration for each pixel on the screen (well, at least half the pixels, since for most coloring schemes the symmetrical properties of the generated image can be taken advantage of). The complexity of the algorithm is in O(n²), where n is proportional to the side length of the image. For implementational details check out the file `scan.c` that comes with AlmondBread.

#### Interleave

This algorithm is significantly faster than the Scan algorithm. However, it can be inaccurate to some extent.

The idea is the following: For each i consecutive pixels, compute only one. If its color matches that of the previous computed pixel, assume all i are of same color and go on. If not, go back one pixel at a time until it matches the "saved" color (plotting each of those different colored pixels along the way, of course). Now, plot the whole line of unique color defined by the first point of "saved" color and the last.

Obviously, this is very fast for images containing large areas of a single color, but may err in images with fine filaments extending through areas of few colors. Thus, it follows the same policy as Fractint's solid-guessing logic. Its complexity is O(c·n²), with 1/i<=c<=1. For further details, check out `interleave.c`.

#### Tesseral

This method is an imitation of Fractint's mode of same name. It is also known as the Mariani-Silver Algorithm, the Quadtree Algorithm, or the Rectangular Subdivision Algorithm. Although in most cases it is not as fast as Interleave, it has a nice theoretical approach and is fun to look at.

It makes use of the Mandelbrot set's connectedness property, which has as a consequence that there cannot be a region of one color inside a rectangle of another color, i.e. once you find that the sides of a rectangular region are of a single color, it is safe to assume the whole rectangle to be of said color (with some margin of error, of course, since strictly speaking one pixel defines a region, not just one point). The procedure used to generate the image is this:

1. Compute the border of the region, keeping track of the number of colors used.
2. If only one color was used, assume the whole region to be of that color and "blit" it in.
3. If more than one color was used, subdivide the image into 2 smaller ones and calculate them recursively (stopping at a certain minimum side length, beyond which use scanning to finish up the image).
The following image shows the behavior of the tesseral algorithm. In the top half, no "blitting" was done, so you can see the outline of the rectangles made up of one color.

There are some details I think are worth elaborating on. First, AlmondBread's algorithm subdivides at least once to avoid "blitting" in the initial image. Secondly, it skips computing those edges of the image that were already done in the previous stage (which is either 3 or all 4 using the "DIN" approach). Thirdly, AlmondBread uses the "DIN" approach to subdividing the image (I call it that because of the way paper sizes in the DIN system are related, e.g. DIN A5 is half the size of DIN A4), which is also known as hv-splitting. This means simply that subdividing in horizontal or vertical direction is done in alternating fashion, which has proven to be very efficient over other methods like 4-way splitting. The complexity of the Tesseral algorithm is in O(n·log(n)). The code can be found in `tesseral.c`.

#### Boundary Tracing

The Boundary Tracing method takes the idea that led to the Tesseral algorithm one step further. The fact that the Mandelbrot set is connected, implies that there won't be areas of one color "inside" areas of a different color. This means, that if you trace the boundary of a region of a single color, i.e. determine the connected set of pixels of one color whose members have at least one neighboring pixel of a different color, all the pixels inside that area will have that one color also. Color in this context refers to the iteration value of a single pixel. The algorithm is this:

We scan the whole screen pixel by pixel, keeping an offscreen map of iteration values for each pixel (initially set to all -1s to signal "not done"). If we hit a pixel not done, we start tracing its border going right and keeping the "wall" to our left until we get back to the starting point. For each pixel, the next pixel visited will be computed thus:

1. check pixel to left (looking in the current direction, i.e. if the current direction is right, the pixel to the left is the one above it). If it's part of the wall (i.e. offscreen or a different color than the one we're tracking), go to 2. If it's not, it is the next pixel.
2. do same for pixel straight ahead. Go there if it's not part of the wall.
3. ditto for pixel to the right. If this pixel is also part of the wall, go back one pixel.

When we're back at the starting point we trace the boundary once more. Now, whenever we move up, we plot pixels to the right of the current one until we hit the wall.

This method was inspired by a post on sci.fractals by Maarten Egmond (cp@stack.urc.tue.nl).

The complexity of the Boundary Tracing algorithm is in O(n2). The code can be found in `boundary.c`.

#### SOI

SOI stands for Simultaneous Orbit Iteration. It was "invented" by Steven Stoft and introduced in his program Fractal Witchcraft. Most of the information I had about the method before starting out came from a response to my inquiry about Fractal Witchcraft's speed from Robert P. Munafo on sci.fractals.

As you may recall from section The Mandelbrot Set, calculating M involves an iterative process, which means on successive iterations one point is mapped onto another one. Now, if we look at not only one point but a whole set of points, say a rectangle, these points get mapped onto a different set of points in the complex plane. If this (iterated) set of points is still a rectangle, we can make the reasonable assumption that in order to obtain the iterated value for a point inside the original rectangle all we have to do is a linear interpolation, i.e. determine where the point in the original rectangle would be in the new, iterated rectangle if we had merely "resized" it and moved it to the new location (which in effect is what we have done by iterating). To put it in other words, if point x is in the middle of the rectangle before iterating it will be in the middle of the iterated rectangle also.

The problem is that for a given rectangle the iterated point set is very rarely a perfect rectangle. The approach used by AlmondBread involves a polynomial interpolation of 2nd degree (using Newton interpolation). For every rectangle we have 9 key points plus four test points. All points are iterated and the test points are interpolated. If the ratio of iterated to interpolated test points is greater than an allowed margin, we subdivide, proceeding in a fashion similar to the Tesseral algorithm (using a 4-way split, though). Otherwise, we do the iterating and interpolating over until we have to split or maxiter is reached, in which case we simply fill the rectangle and exit. This may sound trivial at first, but the whole interpolating stuff can get pretty confusing after a while, particularly when you have to interpolate certain points and don't know which key points to use (just take a look at the subdividing part in the code), but this is mainly due to my quick'n'dirty implementation tactics.

SOI does a very good job on images that have a high number of average iterations, typically very deep zooms that take a long time using other algorithms such as Interleave. The following image is a very deep zoom that took over 2 hours using Fractint under DOS on a 486DX-33. It took about 4 minutes using AlmondBread under Linux on the same machine.

Load the image into Fractint and check out yourself how long this image takes to calculate (Download the Parameter file if you just want to calculate the image).

SOI does not do as good a job on low-magnification images, because unlike Interleave, for example, SOI has a large organizational overhead. There are a few tricks, however, which can make SOI as fast or even faster on low-magnificaction images. The main problem with low-magnification images is that the number of simultaneous iterations is very small and you have to subdivide after very few iterations. Therefore, AlmondBread has a minprogress variable, which controls the minimum number of iterations required to allow subdividing. If this number is not reached, Interleave scanning is employed.

A word of warning is necessary at this point. The SOI in AlmondBread is far from perfect. With certain images you may notice severe distortions caused by a failure of the tolerance code. You can then lower the value of the variable tolerance (possibly from 0.1 to 0.05). Of course, this is not acceptable, but up till now, I have not found a fool-proof way of detecting distortion (I would greatly appreciate any input you might have on this or any other aspect of SOI. Plus, let me know if you have implemented SOI or a similar method yourself). The complexity of SOI cannot be easily compared to previous methods, since it does not depend on the side length of the image. The code for SOI can be found in `soi.c`.

### Iteration Methods

An iteration method determines if a pixel belongs to M or not. If not, it also determines the iteration count, i.e. the number of iterations it took before |z| "bailed out". The goal here is to minimize the computational complexity of the iteration loop, because this is where the most time is spent. For example, for a 1000x1000 image that has an average iteration count of 1000 (not uncommon), the iteration loop is executed a billion times when using simple scanning.

#### Conventional

As you might have guessed, this is the most basic algorithm to perform the iteration of a single point. It has no optimizations (except for obvious ones). It is a straightforward implementation of the method discussed in The Mandelbrot Set. You can find the code in `conventional.h`.

#### Optimized

This routine has two major optimizations, both of which concern themselves with the iteration of points inside M:
Periodicity checking
Recall that a point is in M if during iteration its magnitude does not exceed bailout. We know that this is the case for a certain point z0:=c if we discover that zi+p=zi, i.e. the iteration has entered a cycle. This sounds easy enough, but implementing it proves rather cumbersome since we do not know in advance what p is and saving every zi and testing it against all its predecessors would take too long. Therefore AlmondBread saves only every "maxiter/8"th z, exiting if |z-"saved z"|<pixelsize. This is nowhere near optimal, but suffices to catch a lot of "offenders". If you have a better method, please let me know.
Cardioid checking
If you take a look a the Mandelbrot set, you will notice that most of it is made up of a number of discs that are attached to an object that looks like a sliced apple, which is called a cardioid (presumably because it looks a little like a heart, too). With a couple of operations we can check whether a given point lies inside one of those objects and exit if it does (at least for the main cardioid, and the 3 biggest discs attached to it). For example, if z=(r,i) satisfies (1/4)²<=(r+1)²+i², it lies inside the disc at (-1,0) with radius 1/4. Most of the cardioid checking code in AlmondBread was taken from a post by Todd S. Lehman (toddl@county.lmt.mn.org) on sci.fractals. Also, check out "The Beauty of Fractals" (p.56).
The code can be found in `optimized.h`.

#### Unroll

Besides the optimizations mentioned in the previous paragraph, this method includes a technique that is geared toward reducing the number of instructions used to perform one iteration. The "normal" iteration loop looks something like this:
``````
do
{
im=(im+im)*re+cim;
re=rq-iq+cre;

rq=re*re;
iq=im*im;

/* maybe some periodicity checking here */

} while (i++<=maxiter && rq+iq<bailout);
```
```
As you can see, this process involves 3 floating point multiplies. We would be able to eliminate one of them if we didn't have to check for bailout on each iteration. And that's exactly what we do: We unroll the loop and check for the bound only every n-th time.
``````
do
{
/* do this n times */
i1 = (im+im)*re+cim;
r1 = (re+im)*(re-im)+cre;

/* now check the bound once */
rq = rn*rn;
iq = in*in;

if(rq+iq>bailout) { /* determine correct iter and return */ }

i+=n;
} while (i<maxiter);
```
```

In addition to saving one multiplication (which in fact isn't such a big savings on modern processors such as the Pentium) the comparison and looping overhead can be avoided. This turns out to be an even greater advantage, especially for Intel processors where floating point comparisons are very costly.

There are a couple of problems, though. First of all, we are running the risk of doing too many iterations if we don't check the bound each time, but if maxiter is large compared to the number of unrolled iterations the time-savings outweigh the overhead.

Secondly, the magnitude of z can get very large after it surpasses bailout, generating a floating-point exception. Therefore, AlmondBread unrolls only to a depth of eight, which seems to be safe for IEEE 754 doubles.

Thirdly, it has to be determined when exactly the bailout value was exceeded in order to color the corresponding pixel correctly. Nick Haines (nickh@cmu.edu) suggested doing an "exponent correction hack", which uses the exponent of the magnitude of z to determine the correct iteration number. It is based on the observation that after |z| surpasses bailout, its exponent is doubled on each consecutive iteration. However, this works only for bailout values larger than 4. The code is in `unroll.h`.

### Coloring Schemes

Once the iteration count for a certain pixel is established it can be plotted using a color that is a function of the iteration count. Various coloring schemes can be used to describe this mapping.

#### Normal

This coloring scheme will color a pixel according to the level of iteration where the magnitude of z got bigger than bailout (therefore sometimes termed Level Set Method).

#### Potential

This method yields smoother contours of the color bands than Normal at the cost of slower computation. The equipotential lines computed by this method are lines of equal escape time towards infinity (as are the bands computed by Normal, albeit on a much cruder basis). Consult  and  for mathematical details. The potential is computed as

`pot:=maxcolor-slope*log(|zn|)/2n`,

where `maxcolor` is the maximum value for `pot`. The variable `slope` affects the width of the color bands. zn is the value of z at the end of iteration, i.e. the value it holds when bailout is reached, after n iterations. Values less than 1 are "truncated" to 1.

#### Decomposition

This method determines the color for pixel p on the basis of the angle of zn in polar coordinates, where zn is the value of z at the end of iteration. The normal representation of complex numbers is as a pair of real numbers (r,i) where r describes the real part and i the imaginary one. You can also describe a complex number completely by its radius (magnitude) and angle on a plane with real and imaginary axes. The angle can be computed as follows:

`angle:=(atan2(i,r)/(2*pi)+0.5)*arity`.

This will yield a value in the range 0-`arity`. The `arity` parameter determines how many segments the circle of possible angle values is divided into, i.e. if `arity` is 2, p will be colored according to whether the imaginary part of zn is positive or negative. Again, try  or  for an in-depth explanation. `arity` 256 and a nicely shaded colormap gives beautiful pictures.

### References


Peitgen, H.-O., Saupe, D., et al., The Science of Fractal Images. Springer Verlag, Berlin, 1988.

Taylor, Michael C., et al., The sci.fractals FAQ.

Peitgen, H.-O., Richter, P.H., The Beauty of Fractals. Springer Verlag, Berlin, 1986.

Giffin, Noel, The Spanky Fractal Database.

### Further Information

Mu-Ency - The Encyclopedia of the Mandelbrot Set by Robert Munafo.
Although a little chaotic (but that's what fractals are about, right), this has some very good information on the Mandelbrot set and how to speed up its generation.
Fractal Geometry of the Mandelbrot Set by Robert L. Devaney.
This explains some of the properties of the Mandelbrot set in detail.

Michael R. Ganss (rms@cs.tu-berlin.de)
Last modified: Thu Apr 18 15:01:44 MEST 2002