Mathematical background

For a basic introduction to the Mandelbrot set we recommend this wiki page by Arnaud Chéritat.

This section will only present a brief outline of the algorithms used in fractalshades for deep Mandelbrot zooms. Do not expect rigorous demonstrations. For a more in-depth coverage of these points you should also have a look at Claude Heiland-Allen blog, particularly these entries:

Zoom level and native floating-point precision

We recall the basic iterations fo the Mandelbrot set:

\[\begin{split}\left\{ \begin{array}{rcl} z_0 &= &0 \\ z_{n+1} &= &z_{n}^2 + c \end{array} \right.\end{split}\]

Here, \(c\) represents a pixel in the image. For all interesting areas in the Mandelbrot set, the magnitude of \(|c|\) is close to \(1.0\) - somewhere between \(0.5\) and \(2.0\).

On most computers the highest precision floating-point format available natively is double, with a significand precision of 53 bits. Its means that, at machine precision we have

> 1. + 2 ** -53 == 1.
True

So for 2 pixels separated by less than \(2^{-53}\), a calulation at the highest available machine precision will result in exactly the same sequence. This is the maximal resolution for doubles, roughly \(10^{-16}\). If the image is 1000 pixels width, issues begins to occur for a width of \(10^{-13}\).

To go deeper one need to switch to an arbitrary-precision library (like gmpy2).

Unfortunately calculations will become a lot slower, and moreover this time will increase quickly with the zoom factor. Algorithms have been developped to avoid running all pixels in full precision (in fact, one orbit calculated at full precision and stored at standard precision if often enough)

Perturbation Theory

This technique is based on an original idea from K. I. Martin (see [Mar13])

Let’s assume we know a reference sequence \(Z_n\) (for instance, the center of the image iterated with arbitrary precision):

\[\begin{split}\left\{ \begin{array}{rcl} Z_0 &= &0 \\ Z_{n+1} &= &Z_{n}^2 + C \end{array} \right.\end{split}\]

For any point \(c\) in the image we can express it as \(c = C + \delta c\) where \(\delta c\) is small. Similarly for its iterated sequence we will write: \(z_n = Z_n + \delta z_n\).

Some elementary manipulations yield the recurrence formula for \(\delta z_n\) :

(1)\[\begin{split}\left\{ \begin{array}{rcl} \delta z_0 &= &0 \\ \delta z_{n+1} &= & \delta z_n (2 Z_n + \delta z_n) + \delta c \end{array} \right.\end{split}\]

Similar formulas can be found for the derivatives (useful for instance for distance approximation, or early interior points detection) \(\frac {\partial \delta z_n}{\partial c}\) and \(\frac {\partial \delta z_n}{\partial z_1}\).

The good and a bit unexpected news is that \(\delta z_n\) can be iterated at standard precision and, with a few caveat, still give similar results to a calculation with full precision! A rigorous demonstation is probably quite involving and in any case beyond the scope of this section…

Avoiding loss of precision

When iterating (1) numerical loss of precision may occur, leading to invalid results. This is often evidenced by so-called “glitches” in the image: large zones with uniform color or noisy rendering.

A deeper investigation shows that these issues happen when the calculated orbit \(z_n\) goes close to zero [1], under the influence of a nearby minibrot of period \(n\).

Several methods have been investigated to avoid such glitches.

We present here the one used in fractalshades which is based on an original idea from Zhuoran in 2021 (see [Zhu21]). Whenever \(|z_n| <= |Z_n|\) we do a rebase i.e. basically restart the reference orbit at index 0:

(2)\[\begin{split}\left\{ \begin{array}{rcll} \delta z_{n+1} &\leftarrow &z_n & \\ Z_{n+1+i} &\leftarrow &Z_i &\forall i \in \mathbb{N} \end{array} \right.\end{split}\]

Of course from an implementation point of view it is easier to just reset the index \(n\) in \(Z_n\) to 0 - rather than making the substitution above.

The same strategy can also be applied whenever the reference orbit diverges.

Billinear approximations

When \(|\delta z_n| \ll |2 Z_n|\), (1) can be approximated without loss of precision by a billinear expression in \(\delta z_n\) and \(\delta c\) :

\[\delta z_{n+1} = 2 Z_n \delta z_n + \delta c = A_n \delta z_n + B_n \delta c\]

Or, with more precise notation highlighting that this set of billinear coefficients map \(\delta z_n\) to \(\delta z_{n+1}\):

(3)\[\delta z_{n+1} = A_{n, n+1} \cdot \delta z_n + B_{n, n+1} \cdot \delta c \quad \forall \, \delta z_n, \, |\delta z_n| \ll |2 Z_n| = r_{n, n+1}\]

Now, the composition of 2 bilinear functions is still a billinear function. Combining 2 by 2 gives:

(4)\[\begin{split}\left\{ \begin{array}{rcll} A_{k, m} &= & A_{k, l} \cdot A_{l, m}& \\ B_{k, m} &= & B_{k, l} \cdot A_{l, m} + A_{l, m} \end{array} \right.\end{split}\]

Valid if:

\[\begin{split}\begin{array}{rcl} |\delta z_n| &\leq & \min \left( r_{k, l}\, , \: \max \left( 0\, , \: r_{l, m} - \frac{|A_{k, l}| |\delta c|}{|B_{k, l}|} \right) \right) \\ &\leq & r_{k, m} \end{array}\end{split}\]

It means, we can combine 2 by 2 such approximations (and their validity bound \(r_{k, m}\)) and pre-compute the results - for instance in a binary tree structure. This gives a set of precomputed shortcuts allowing to skip many iterations, when the validity conditions are met.

Note: The same approach can also be applied to the derivatives, they just becomes a simple linear relationship.

Extended range floating points

Another limitation of native double format is met when the size of a pixel drops below approx \(10^{-308}\). This time, it is the range of the exponent.

For IEEE 754 double, exponent has 11 bits and its valid range is \(-1022\) to \(+1023\). Numbers below \(2^{-1022}\) are first stored by compromising precision (subnormal numbers), down to about \(5 \cdot 10^{-324}\) and then they vanish:

> 10. ** -324 == 0.
True

In short, manipulating numbers below \(10^{-308}\) in double format shall be avoided. fractalshades switches to a ad-hoc datatype based on a couple (float64, int32) and the main operations have been implemented for this datatype.

Finding minibrots

In order to find a minibrot center and size, the process is the following:

  • find the period first by iterating a small ball until it contains the origin (using ball arithmetic),

  • then use Newton’s method to find the location of the minibrot with the computed period

  • then compute the size (and orientation) using a renormalization formula (as the one found in [BRH97])

We will not describe in details here these techniques which have already been better explained elsewhere. The interested reader might refer to:

  • Van den Doel [vdD18]

  • Heiland-Allen [HA16].

A generalisation of these techniques for the Burning Ship is described in Heiland-Allen [HA19].

An implementation of these algorithms is available from the GUI (right-clik on the image) for the fractals supporting arbitrary-precision.

Going further: a few references

[BRH97]

Edward Ott Brian R. Hunt. Structure in the parameter dependence of order and chaos for the quadratic map. J. Phys. A: Math. Gen. 30 (1997) 7067–7076, 1997. https://terpconnect.umd.edu/ bhunt/papers/wind.html.

[HA13]

Claude Heiland-Allen. Perturbation techniques applied to the mandelbrot set. 2013. https://mathr.co.uk/mandelbrot/perturbation.pdf.

[HA16]

Claude Heiland-Allen. Deriving the size estimate. 2016. https://mathr.co.uk/blog/2016-12-24_deriving_the_size_estimate.html.

[HA18]

Claude Heiland-Allen. Newton method for periodic points. 2018. https://mathr.co.uk/blog/2018-11-17_newtons_method_for_periodic_points.html.

[HA19]

Claude Heiland-Allen. At the helm of the burning ship. In Proceedings of EVA London 2019. 2019. http://dx.doi.org/10.14236/ewic/EVA2019.74.

[HA21]

Claude Heiland-Allen. Deep zoom theory and practise. 2021. https://mathr.co.uk/blog/2021-05-14_deep_zoom_theory_and_practice.html.

[HA22]

Claude Heiland-Allen. Deep zoom theory and practise again. 2022. https://mathr.co.uk/blog/2022-02-21_deep_zoom_theory_and_practice_again.html.

[Har07]

Jussi Härkönen. On smooth fractal coloring techniques. 2007. https://www.scribd.com/document/91355075/On-Fractal-Coloring-res.

[Mar13]

K.I. Martin. Superfractalthing maths. 2013. http://www.science.eclipse.co.uk/sft_maths.pdf.

[vdD18]

Kees van den Doel. Notes on mandelbrot set. 2018. http://persianney.com/fractal/fractalNotes.pdf.