The algorithms are presented using a C-like pseudocode that you should be able to adapt to the language of your choice. Even Python.

The code, as presented, treats complex numbers as individual (x,y) components
and performs primitive arithmetic operations on those. Some bits might say
simply "calculate |z|", but you should know how to compute the modulus from the
individual components without having it explicitly spelled out. If you're using
a language that has complex numbers as a primitive type, such as C99, then you
should be able to implement many of these algorithms much more simply - you
can divide complex numbers using the '/' operator and the compiler will generate
all the boilerplate and presumably do so very efficiently. If, however, you're
using a language where complex number support is via a library that gives you
Complex objects with arithmetic methods called "add", "conjugate", "multiply",
etc. that return new instances of the Complex class, steer well clear. Fractal
generation typically requires billions of floating point calculations. This is
not a problem - floating point performance these days is measured in gigaflops.
However, creating billions of objects on the heap and then freeing them when
they go out of scope *is* and will make your runtime hours rather than
seconds. Using a built-in list type should be avoided as well; even if you're not
reallocating Array objects, storing a value in an Array is orders of magnitude
slower than a primitive arithmetic operation. Inside the iteration loops, you
should also avoid function calls where possible, since the overhead of tens of
millions of function calls is noticeable. That said, be sensible: it is much
faster to calculate the exponent via a library call than to calculate the
Taylor series approximation using basic arithmetic.

Things that are in ALL_CAPS are assumed to defined outside the functions (e.g. as object properties).

Smooth colouring requires the ability to obtain a colour in between two colours picked from a discrete palette:

function: getIntermediateColourValue parameters: c1(R,G,B), c2(R',G',B'), m: 0 <= m <= 1 return: [R + m * (R' - R), G + m * (G' - G), B + m * (B' - B)]

Fractal images are typically coloured according to the iteration count after some bailout condition has been met. It is important not to use too few colours. Otherwise images of the fractal boundary will have a different, unrelated colour for every pixel, since the iteration counts tend to vary wildly. This will produce colour noise which will destroy detail and result in an ugly image.

define: paletteObject properties: controlColours: Array(0..2N); RGB start points and end points of N colour channels weightingPerGroup: Array(0..N); relative importance of each colour group; only used for the first cycle coloursPerChannel: integer; number of colours in a group when the weighting is zero

Choose your control colours. If you're a colour theory expert, this should be a doddle. If not, don't use arithmetical relationships between R, G and B values since this will produce a really bad palette. What I do is use a decent colour picker that shows colours by their perceptual value. I then trace a diagonal line from a dark tone to a light tone and use the two end values as a start point and end point of a colour group. This provides lots of intermediate, perceptually related tones.

The weighting array is used to tune the palette for the fractal image. Typically, you'll want the first group transition to occur at the fractal boundary. You can use a negative weighting value to create a smaller colour group. You'll also probably want to assign a higher weighting to green tones than red / orange tones, since the eye is more sensitive to greens than reds and oranges.

A paletteObject is passed to a palette expansion function. The expansion function ensures two things:

- That there is one palette entry per iteration to avoid colour noise
- That colours don't change when the user increases the iteration count to reveal more detail in the image.

function: makePalette parameters: paletteDef: paletteObject; maxColours: integer >= 0 to limit the palette size with 0 signifying "unlimited" let: p = [], toGenerate = maxColours > 0 : min(MAX_ITER, maxColours) : MAX_ITER, controlColours = paletteDef.controlColours, weightingPerGroup = paletteDef.weightingPerGroup, coloursPerChannel = paletteDef.coloursPerChannel, totalBias = 0, weighting2 = [], bias2 = 0, groupCount = min(weightingPerGroup.length, controlColours.length / 2) generated = 0, thisRound = min(toGenerate, DEF_MAX_ITER); // Calculate total bias and set up even weighting for subsequent rounds for (i = 0..groupCount - 1) { totalBias += weightingPerGroup[i]; weighting2[i] = 1; ++bias2; } while (generated < toGenerate) { let: coloursPerPointOfBias = ceil((thisRound - groupCount * coloursPerGroup) / totalBias); for (i = 0..groupCount - 1) { let: idx = 2 * i, c1 = controlColours[idx], c2 = controlColours[idx + 1], weighting = weightingPerGroup[i], // The number of colours to create in this channel colours = coloursPerGroup + coloursPerPointOfBias * weighting, fInterval = 1 / colours; // Add start point push(p, c1); // Intermediate points for (j = 1..colours - 1) { push(p, getIntermediateColourValue(c1, c2, j * fInterval); } // End point push(p, c2); } generated += thisRound; // Subsequent rounds use even weighting weightingPerGroup = weighting2; totalBias = bias2; // Each subsequent round generates twice as many colours as the last thisRound *= 2; } // White and black for set interior points push(p, WHITE, BLACK); return p;

Certain colouring operations - basins of attraction, distance shading - need to be done in the HSV domain. Therefore, we need to be able to convert RGB to HSV and HSV to RGB.

function rgb2hsv parameters: r, g, b: integer in the range 0 - 255 // Convert to values between 0 and 1 such that 255 => 1.0 r = (r & 0xFF) / 255; g = (g & 0xFF) / 255; b = (b & 0xFF) / 255; let: Min = min(r, g, b), Max = max(r, g, b); if (Max == Min) { // Greyscale return [0, 0, Max]; } let: h = 0, d = Max - Min, s = d / Max, v = Max; // Get the sector in which the hue falls; r g b are 120° apart when Max is r: h = (g - b) / d; g: h = (b - r) / d + 2; b: h = (r - g) / d + 4; // Normalize h so that 0 <= h < 1 if (h < 0) { h = h + 6; } return [h / 6, s, v];

function hsv2rgb parameters: hue, saturation, value: real: normalized such that 0 <= v <= 1 if (hue == 1) { // Hue is an angle around a colour wheel hue = 0; } // Colour wheel is divided into six sectors hue = hue * 6; let: C = saturation * value, inter = hue % 2 - 1, m = value - C; if (inter < 0) { inter = -inter; } let: X = C * (1 - inter), _r = 0, _g = 0, _b = 0, sector = floor(hue); when sector is 0: _r = C; _g = X; 1: _r = X; _g = C; 2: _g = C; _b = X; 3: _g = X; _b = C; 4: _r = X; _b = C; 5: _r = C; _b = X; return [(_r + m) * 255, (_g + m) * 255, (_b + m) * 255];

If we are able to calculate a distance estimate to the fractal boundary we can use it to adjust the v-value in HSV colour space so that pixels closer to the boundary are darker. This allows us to pick out the boundary without having explicitly to colour pixels black. We use the boundaryFraction parameter to emphasize (small values) or de-emphasize the boundary.

function colourDistance parameters: h, s, v: real: components of the HSV colour space; distance: real: the distance estimate calculated during iteration; pixelWidth: real: the width represented by a pixel on the screen; boundaryFraction: real: value > 0 used to control prominence of the boundary let: dScale = log_{2}(distance / pixelWidth), factor = 0; if (dScale > 0) { factor = 1; } else if (dScale > -boundaryFraction) { factor = (boundaryFraction + dScale) / boundaryFraction; } // Darken v when factor < 1 return hsv2rgb(h, s, v * factor);

The iteration function calculates an index into a palette which is generally real-valued allowing for continuous colouring and a distance estimation from the set boundary. The distance estimate may be null, such as when the Newton iteration function fails to converge. Given a palette (generated using the algorithm outlined above) we can use these two values to obtain a colour value for the pixel.

function lookupColour parameters: iter: real: the iteration count before bailout (negative value signifies set interior); distance: real: the distance estimate calculated during iteration; p: array: RGB triplets let: // Last two colours are used for set interior maxC = length(p) - 3, r = 0, g = 0, b = 0, idx1 = 0, idx2 = 0, frac = 0; if (iter >= 0) { idx1 = floor(iter) % maxC, idx2 = idx1 + 1; if (idx2 > maxC) { idx2 = 0; } frac = iter - floor(iter); } else { // Set interior idx1 = idx2 = length(p) + iter; } (r, g, b) = getIntermediateColourValue(p[idx1], p[idx2], frac); if (distance) { (h, s, v) = rgb2hsv(r, g, b); (r, g, b) = colourDistance(h, s, v, distance, PIXEL_WIDTH, BOUNDARY_FRACTION); } // One of the colour options for a dark interior / light boundary if (REVERSED_COLOURS) { r = 255 - r; g = 255 - g; b = 255 - b; } return [r, g, b];

This takes an integer iteration value and turns it into a real value that can be fed into the lookupColour routine. It takes the escape value into account when the bailout condition has been reached. Identify the order of the polynomial, P. If P <= 0 we can't use this function because of the need to take a logarithm; in this case we simply use the integer value.

function: normalizeIterationCount_rp parameters: r: the modulus of z_{n}; logp: the natural logarithm of P; n: the iteration count at bailout let: inter = log_{2}(r); if (inter <= 0) { return n; } let: nu = 1 - log(inter) / logp; return max(n + nu, 0);

To use it, calculate the iteration count, distance and z_{n} as normal, then:

let: r = sqrt(x_{n}* x_{n}+ y_{n}* y_{n}), niter = normalizeIterationCount_rp(r, log(P), iter), return lookupColour(niter, distance, palette);

This is used as part of the basin colouring algorithm (see below). It is strictly correct only when the point converged on is a distinct root of the function. It can be used for general-purpose colouring but I personally prefer the results from the exponential smoothing algorithm.

Calculate the difference between two iteration values when the convergence criterion has been met. Then:

let: q = log(log(diff) / log(CONVERGENCE_RADIUS)) / log(2)

This produces a value between 0 and 1 that should be subtracted from the iteration count.

This is the algorithm used for convergent fractals. The colour value needs to be calculated as part of the iteration loop:

let: s = 0, diff = 1, distance = 0; while (ITERATING) { s += exp(-1 / diff); Calculate diff = |(z_{k}- z_{k-1})|; Calculate distance estimate } return lookupColour(s, distance, PALETTE);

An option for Newtonian fractals is colouring according to which root of the function a point in the complex plane converges to. This algorithm allows for an unlimited number of basins to be coloured, although there is no guarantee that each basin will have a perceptually distinct colour. Colours are distinct with 6th order polynomials and distinguishable up to 10th order.

function makeBasinColour parameters: rootIdx: the index of the root found by convergence; rootCount: the number of distinct roots of the function; iter: the number of iterations required for convergence; diff: the absolute difference between the last two iterated values; distance: the distance estimate calculated during iteration let: // Add to the hue angle to make the basin for the principal root an attractive colour hueAngle = rootIdx / rootCount + 0.75; // Ensure that adjacent roots don't have similar colours if (rootIdx % 2) { hueAngle = hueAngle + 0.5; } if (hueAngle >= 1) { hueAngle = hueAngle - 1; } // Base value on normalized iteration count for a super attracting fixed point let: q = log(log(diff) / log(CONVERGENCE_RADIUS)) / log(2), v = 1 - log(iter - q) / log(MAX_ITER); return colourDistance(hueAngle, 1, v, distance, PIXEL_WIDTH, BOUNDARY_FRACTION);

The roots of the function are assumed to be sorted according to the complex number sorting algorithm given below.

We choose a fudge factor that makes the principal root (the one closest to the positive real axis) a colour that we find pleasing and go on from there. 0.75 corresponds to an angle of 3 * π / 2 and results in a nice shade of purple. A fudge factor of 0 gives pure red. Adjust to taste. Since the angles are evenly spaced around the colour wheel, a large number of roots will result in perceptually indistinct colouring for adjacent basins (since they'll all end up in the same π / 3 sector). Therefore, we make the angle of alternating basins opposite by adding π.

We only use two dimensions of a three-dimensional colour space as I've not been able to figure out a meaningful saturation calculation that produces good results. The value is based on iteration count - the fewer the brighter which makes the roots stand out as distinct bright spots. This also means the more the darker, so zooming into a singularity will make things look distinctly dim. Then again, zooming into a singularity is a fool's errand since there's nothing to see.

Passing the HSV co-ordinates and the distance estimate to the colourDistance routine picks out the boundary very nicely indeed.

Distance estimation is a method to calculate the distance of a point in the complex plane from the boundary of the fractal. We use this distance estimate to darken pixels close to the border. Boundary tracing controls noise, picks out detail and increases the subjective image quality.

We will only consider *exterior* distance estimation. Interior distance
estimation for filled sets is possible, but has a prerequisite of identifying the period length
which can be computationally prohibitive for points very close to the set
boundary. Plus some seriously scary maths once we've done so.

This algorithm is for divergent fractals where infinity is a super-attracting
fixed point and requires that we know the differential of the
iteration function *f*.
Since we're dealing with simple powers of *z*, this is straightforward. For
powers of z >= 2, the differential of *f* is:

2z, 3z^{2}, 4z^{3}, 5z^{4}, etc.

For powers of z <= -2, we have:

-2z^{-3}, -3z^{-4}, -4z^{-5}, etc.

For intermediate powers of z, use the general rule:

d(z^{P}) ----- = Pz^{P-1}dz

The differential has to be calculated as part of the iteration loop. For Julia-style
iterations, we are calculating the partial differential with respect to *z*
while for Mandelbrot-style iterations the calculation is with respect to *c*.
Fortunately, we can handle this difference with a simple change in initial values
and a constant of addition. For illustration, I'll show the calculation for the
z^{2} + c iteration function. Important:
we have to calculate the differential using the results from the last iteration,
otherwise we'll get things badly wrong.

Parameters: x, y, cX, cY let: x_ = x, y_ = y, iter = 0, eV = 256, rsquared = (x_ * x_ + y_ * y_), dzx = ISJULIA ? 1 : 0, dzy = 0, dzPlus = ISMANDEL ? 1 : 0, // Initialize distance to the null value distance = null; while (iter < MAX_ITER AND rsquared < eV) { // Update dz tmp = 2 * (dzx * x_ - dzy * y_) + dzplus; dzy = 2 * (dzx * y_ + dzy * x_); dzx = tmp; // Update z tmp = x_ * x_ - y_ * y_ + cX; y_ = 2 * x_ * y_ + cY; x_ = tmp; rsquared = x_ * x_ + y_ * y_; ++iter; }

At the end of the loop, we have calculated the differential. Now we can calculate the distance estimation:

if (iter < MAX_ITER) { Let: r = sqrt(rsquared), dz = sqrt(dzx * dzx + dzy * dzy); // This will be fed into lookupColour distance = log(r * r) * r / dz; }

Newtonian fractals are generated by applying the generalized Newton-Raphson root finding method to points in the complex plane:

f(z_{k}) z_{k+1}= z_{k}- R·------ fʹ(z_{k})

The algorithm outlined here assumes that *f(z) = 0*
has *n* distinct roots where *n* is the order of *f*. This
makes the roots of *f* super-attracting. If this is not so, the world won't
end but the distance estimation will be off. For example:

zonly has three distinct roots, despite being fourth order. It is the product of:^{4}- (1+i)z^{3}- (1-i)z^{2}- (1-i)z -i

(z + 1)(z +so that +1 is a root with a multiplicity of 2, making it attracting but not super-attracting (and makes Newton's method slow to boot).i)(z - 1)^{2}

For shorthand, let's call Newton's method *g(z)*. To differentiate it
we need to apply the quotient rule for differentiating a function of the form *p/q*:

pʹq - pqʹ --------- q^{2}

Noting that q is *fʹ(z)* and that qʹ is *fʹʹ(z)* (the second differential
of *f*), we can obtain the formula for differentiating *g(z)*:

fʹ(z)^{2}- f(z)fʹʹ(z) g'(z) = 1 - R·------------------- fʹ(z)^{2}

function dgz parameters: fzx, fzy: the components of f(z); dfzx, dfzy: the components of fʹ(z); d2fzx, d2fzy: the components of fʹʹ(z); rX, rY: the components of the relaxation parameter let: denomX = dfzx * dfzx - dfzy * dfzy, denomY = 2 * dfzx * dfzy, nomX = denomX - (fzx * d2fzx - fzy * d2fzy), nomY = denomY - (fzx * d2fzy + fzy * d2fzx); // Multiply by R tmp = nomX * rX - nomY * rY; nomY = nomY * rX + nomX * rY; nomX = tmp; // Multiply by conjugate of denominator tmp = nomX * denomX + nomY * denomY; nomY = nomY * denomX - nomX * denomY; nomX = tmp; denomX = denomX * denomX + denomY * denomY; return [1 - nomX / denomX, nomY / denomX];

The cumulative differential of *f(z)* is calculated during the main iteration
loop. Note that you have to keep track of the *last* value of the cumulative
differential as well as the current one. Most references will indicate |dz| in
the final distance estimate calculation, but this is several orders
of magnitude out. The correct factor is |dz - dzlast| (credit: David Makin).

Parameters: x, y Constants: rX, rY: the components of the relaxation parameter let: x_ = x, y_ = y, diff = 1, dzx = 1, dzy = 0, dzxl = 0, dzyl = 0, eV = 1e-5, iter = 0, distance = null, s = 0; while (iter < MAX_ITER AND diff > eV) { // Update colour index s = s + exp(-1 / diff); Compute [nomX, nomY] = f(z); Compute [denomX, denomY] = fʹ(z); Compute [d2zx, d2zy] = fʹʹ(z); // Update differential [_dzx, _dzy] = dgz(nomX, nomY, denomX, denomY, d2zx, d2zy, rX, rY); dzxl = dzx; dzyl = dzy; tmp = dzx * _dzx - dzy * _dzy; dzy = dzx * _dzy + dzy * _dzx; dzx = tmp; // Multiply by R tmp = nomX * rX - nomY * rY; nomY = nomY * rX + nomX * rY; nomX = tmp; // Multiply top and bottom by conjugate of denominator tmp = nomX * denomX + nomY * denomY; nomY = nomY * denomX - nomX * denomY; nomX = tmp; denomX = denomX * denomX + denomY * denomY; if (!denomX) { return BLACK; } // Compute g(z) _x = x_ - nomX / denomX; _y = y_ - nomY / denomX; // Get the differences from last round diffX = x_ - _x; diffY = y_ - _y; diff = sqrt(diffX * diffX + diffY * diffY); // Set up the next round x_ = _x; y_ = _y; ++iter; } if (diff < eV) { // Calculate distance estimate let: dz = sqrt((dzx - dzxl) * (dzx - dzxl) + (dzy - dzyl) * (dzy - dzyl)); distance = diff * -log(diff) / dz; } return lookupColour(s, distance, PALETTE);

Implementation note: for the convergent formulae, if the iteration loop
ends without the escape value being reached, we don't automatically return the set interior
colour (which basically means "no idea"). Injudicious values of *R* can
prevent Newton's method converging though over/under relaxation. Nonetheless,
the exponential smoothing algorithm still provides a usable colour value so we
can produce (sometimes striking) images even when the iteration count is
maximum. Assuming that we haven't screwed up Newton's method (*R* is
sensible) and we know the roots that the iteration *should* converge to,
we may still find that some starting values fail to converge on a root. There
are certain combinations of coefficients and starting values where Newton's method
ends up in a cycle. In this case, returning the set interior colour is the correct thing
to do.

When finding roots of a function, we have no expectation of strict equality, since the last value of an iteration is only an approximation of the actual root within some limit of tolerance. Therefore, we need to be able to compare complex numbers for approximate equality.

function complexCompare parameters: x1, y1, x2, y2: the values to compare; tolerance: the acceptable radius within which values are deemed to be equal let: _x = x1 - x2, _y = y1 - y2, delta = sqrt(_x * _x + _y * _y); return (delta < tolerance);

This section assumes a parameterized sort method on a list where one of the parameters is a pointer or reference to a comparison function. It also assumes a list of points [x, y] representing values in the complex plane and that the semantics of comparing two values a, b are:

< 0: a sorts before b > 0: b sorts before a 0: a and b are equal

It also assumes an *atan2* function that returns an angle *theta*:

Upper two quadrants of the complex plane: 0 <= theta <= π Lower two quadrants: 0 > theta > -π

The sort order for complex numbers is undefined. We cannot say which of
-1-3*i* and -2-1*i* is more negative than the other. So we'll
define our own sort order. Here's mine:

The basin colouring algorithm for the Newtonian fractal uses roots sorted using the following algorithm:

let: theta1 = atan2(a.y, a.x), theta2 = atan2(b.y, b.x), r1 = sqrt(a.x * a.x + a.y * a.y), r2 = sqrt(b.x * b.x + b.y * b.y); if (theta1 < 0) { theta1 += 2 * PI; } if (theta2 < 0) { theta2 += 2 * PI; } let: ret = theta1 - theta2; if (0 == ret) { ret = r1 - r2; } return ret;

Implementation note: compute the angles and moduli *before* sorting
and sort objects with *r* and *theta* properties.

Basin colouring for Newtonian fractals requires us to know the roots of the function before we generate the fractal. We can also anticipate a future requirement of identifying critical points of a function for generating Mandelbrot sets. Therefore, we need a generic root finding algorithm. This algorithm uses - what else? - Newton's method.

It assumes the existence of a random function that returns a real value 0 <= x < 1.

definition rootObject properties: x: the real component of z; y: the imaginary component of z; r: the modulus of z; theta: the angle of z

function findRoots parameters: f: a pointer to a function that takes a complex number, z, and computes f(z); fprime: a pointer to a function that takes a complex number, z, and computes fʹ(z) constants: XMAX, XMIN, YMAX, YMIN: the limits of the complex plane to be tested let: eV = 1e-8, // For checking suspiciously small real / imaginary components delta = eV * 100, // For comparing roots found tolerance = delta * 100, rSpan = XMAX - XMIN, iSpan = YMAX - YMIN, // If we don't find a root after n iterations, give up and move on iterations = 100, // We'll test this many points in the x- and y- axes points = 30, rInt = rSpan / points, iInt = iSpan / points, roots = []; // Test a reasonably large sample of points in the plane; assuming that we only have to // do this when we regenerate a fractal, we can use a conservatively large value to minimize // the chances of failing to find distinct roots that are close to each other for (i = 0..points) { for(j = 0..points) { let: x = XMIN + i * rInt, y = YMIN + j * iInt, diff = 1, iter = 0; while (iter < iterations AND diff > eV) { // Get the numerator and denominator [nomX, nomY] = f(x, y); [denomX, denomY] = fprime(x, y); // Divide f(z) by f'(z) - multiply through by conjugate of denominator tmp = nomX * denomX + nomY * denomY; nomY = nomY * denomX - nomX * denomY; nomX = tmp; denomX = denomX * denomX + denomY * denomY; if (!denomX) { // Try another number break; } // Now we can divide _x = x - nomX / denomX; _y = y - nomY / denomX; // Get the differences from last round diffX = x - _x; diffY = y - _y; diff = sqrt(diffX * diffX + diffY * diffY); // Set up the next guess x = _x; y = _y; ++iter; } if (diff < eV) { let: newRoot = TRUE; // Zero suspiciously small real / imaginary components so that the roots can // be sorted deterministically if (abs(x) < delta) { x = 0; } if (abs(y) < delta) { y = 0; } for (k = 0..LENGTH(roots)) { if (complexCompare(x, y, roots[k].x, roots[k].y, tolerance)) { newRoot = FALSE; break; } } if (newRoot) { push(roots, rootObject{ x : x, y : y, r : sqrt(x * x + y * y), theta : atan2(y, x) }); } } } } SORT roots; return roots;

The idea behind the algorithm is that multiple independent threads of execution will find the same roots. The algorithm will treat distinct roots that differ by less than the tolerance threshold as multiple identical roots. Since any function for which this was true would be both stupid and deliberately constructed to have very close roots, we'll not care.

Back to the index