mathr / blog / #

In a previous post I wrote about some symbolic algebra code that extracts the most parallelism possible for Mandelbrot set series approximation iterations. Unfortunately its time complexity was around \(O(n^{4.5})\) making it prohibitively slow to go beyond order 64. I was enlightened by discussions with hapf from fractalforums.com.

First it turns out that computing the series approximation coefficients to full precision is unnecessary, you can get away with 53 bits as provided by machine double, provided the exponent is extended, because the values get huge. This means the potential gains from parallelism are greatly reduced, because the overhead remains the same while the amount of work needed drops hugely.

Secondly, and more importantly, there is a simple formula for the series approximation coefficient iterations, which means that all the cleverstupid symbolic algebra can be done away with and replaced with a couple of nested loops, with the bonus that the order can be changed at runtime without needing to have generated code beforehand.

For the series approximation coefficients (using the notation from the previous post, extended with \(b\) for the series approximation of the derivative for distance estimation):

\[\begin{aligned} \left<\left< z_n \right>\right> &= \sum_{k=1}^\infty{a_{k,n} {\left<\left< c \right>\right>}^k} \\ \left<\left< z'_n \right>\right> &= \sum_{k=1}^\infty{b_{k,n} {\left<\left< c \right>\right>}^k} \end{aligned}\]

the iterations become:

\[\begin{aligned} a_{1,n+1} &= 2 z_n a_{1,n} + 1 \\ a_{k,n+1} &= 2 z_n a_{k,n} + \sum_{j=1}^{k-1} a_{j,n} a_{k-j,n} \\ b_{k,n+1} &= 2 \left(z_n b_{k,n} + z'_n a_{k,n} + \sum_{j=1}^{k-1} a_{j,n} b_{k-j,n}\right) \end{aligned}\]

The sum for the \(a\) coefficients has some redundancy, with terms multiplied in both orders - this can be optimized using case analysis for odd and even coefficient indices (left as an exercise).

Finally, hapf reassured me that it's normal for the number of per-pixel iterations that series approximation can skip reaches a plateau when zooming deeper towards a particular reference point, where doubling the order of the approximation gets you an extra period's (of the reference) worth of skipping.

I implemented this in mandelbrot-perturbator but there are some issues with glitches with high orders at low zoom levels, so I think I'll make the order be determined automatically by the size of the reference minibrot.