Perturbation theory to speed up Julia fractal drawing

Question

I have a really underpowered platform here and want to draw a Julia (and possibly Mandelbrot and Burning ship) fractals using a 8.24 fixed point class. I use iteration count for coloring and need to speed up rendering considerably. I read about perturbation theory (see this) and am trying to do the following:

• Iteratively calculate one image pixel $X$ using $X_{n+1}=X_n+X_0$ until the iteration exits due to iteration limit or distance estimate ($|X_n| >= 2$). I calculate and store all values of $X_n$, $A_n$, $B_n$, $C_n$ on the way.
• Then for the next 3 pixels I use $\Delta_n = A_n\delta_i+B_n\delta_i²+C_n\delta_i³$, where $\delta_i = pixel\space index *{horizontal\space range \over nrOfPixels}$ (the "horizontal step value"). Then supposedly the (estimated) value for the new pixel should be $Y_n = X_n + \Delta_n$. This is where the speedup should come from.
• Then I use $Y_n$ and decide whether to iterate further if $|Y_n| < 2$, or if $|Y_n| >= 2$ "step backwards" (calculate $Y_{n-1}$ like above), to see if $Y_n$ would have exited in an earlier iteration (until $|Y_{n-1}| < 2$).

This should make the 3 pixels "I cheat on" much cheaper to calculate (if this works I'd use 4x4 pixel chunks). Question is: Should this work in theory? It is not really working for me and I'm trying to find out what I'm doing wrong. "Regular" per-pixel Julia rendering is working ok for 16.16 and 8.24 fixed-point formats, so I suspect it's not calculation errors piling up... I can post C++ code if that helps.

Answer 1

As far as I know, the perturbation techniques when applied to deep zooming have only been shown to work reliably for "Mandelbrot" iterations of the form $z := F(z) + c$ where $c$ is the pixel coordinates. I don't know if they can be expected to work for "Julia" iterations where it is the initial $z$ that is the pixel coordinates and $c$ is a constant.

Your first bullet point is fine, that's generally how it is done.

Your second bullet point is mostly ok, for Mandelbrot type iterations you have per pixel $\delta = c_\text{pixel} - c_\text{ref}$, for Julia sets I suppose you replace $c$ with $z_0$. Re-deriving the update formulae for $A_n,B_n,C_n$ because they will be different to the Mandelbrot case:

$$\begin{aligned} Z_{n+1} &= Z_n^2 + C \\ Z_{n+1}+z_{n+1} &:= (Z_n+z_n)^2 + C \\ z_{n+1} &:= 2 Z_n z_n + z_n^2 \\ \delta &= z_0 \\ z_n &= A_n \delta + B_n \delta^2 + C_n \delta^3 \\ \end{aligned}$$ Substituting the last line into the third line and equating coefficients of $\delta^k$ gives $$\begin{aligned} A_0 &= 1 & A_{n+1} &= 2 Z_n A_n \\ B_0 &= 0 & B_{n+1} &= 2 Z_n B_n + A_n^2 \\ C_0 &= 0 & C_{n+1} &= 2 Z_n C_n + 2 A_n B_n\\ \end{aligned}$$ which you can see breaks down if ever $Z_n = 0$.

Your third bullet point is very problematic - the series is only approximate, and beyond a certain iteration count (dependent on location) it goes wild and all bets are off. You can use probe points on the boundary of the image (iterated in another way, and compared with the approximation) to detect when the series is starting to diverge from reality. This will typically be significantly before the pixel escapes to infinity.

One issue with fixed point is underflow and ensuing loss of precision. If $\delta < \frac{1}{10}$ then $\delta^3 < \frac{1}{1000} \approx 2^{-10}$ so you have already lost almost half of your precision. On the other side, the coefficients $A_n, B_n, C_n$ get increasingly large, so you will probably overflow your limited range.

One possible solution is to rescale all your numbers, so instead of multiplying big $C_n$ with small $\delta^3$ you take out a constant factor $K$ to make $\epsilon = K \delta$ have magnitude approximately equal to $1$. It becomes:

$$ \frac{A_n}{K} \times K \delta + \frac{B_n}{K^2} \times K^2 \delta^2 + \frac{C_n}{K^3} \times K^3 \delta^3 = a_n \epsilon + b_n \epsilon^2 + c_n \epsilon^3$$ with all the variables on the right hand side having magnitude near $1$, so no risk of overflow or underflow.

It's a matter of careful algebra to work out the update iterations for $a_n, b_n, c_n$, you will need to involve $K$ at certain points when multiplying scaled values, perhaps something like this: $$\frac{X}{K} \times_1 \frac{Y}{K} = \frac{X\times_1Y/_1K}{K} \equiv x \times_K y = x\times_1y/_1K$$

I suggest instrumenting your program to tell you loudly when overflow occurs, underflow is harder to deal with I suppose.