mathr / blog / #

Start with the iteration formula for the Burning Ship escape time fractal:

\[ X := X^2 - Y^2 + A \quad\quad Y := \left|2 X Y\right| + B \]

Perturb this iteration, replacing \(X\) by \(X+x\) (etc) where \(x\) is a small deviation:

\[ x := (2 X + x) x - (2 Y + y) y + a \quad\quad y := 2 \operatorname{diffabs}(X Y, X y + x Y + x y) + b \]

Here \(\operatorname{diffabs}(c,d)\) is laser blaster's formula for evaluating \(|c + d| - |c|\) without catastrophic cancellation; see my post about perturbation algebra for more details.

Now replace \(x\) and \(y\) by bivariate power series in \(a\) and \(b\):

\[ x = \sum_{i,j} s_{i,j} a^i b^j \quad\quad y = \sum_{i,j} t_{i,j} a^i b^j \]

To implement this practically (without lazy evaluation) I pick an order \(o\) and limit the sum to \(i + j \le o\). Substituting these series into the perturbation iterations, and collecting terms, gives iteration formulae for the series coefficients \(s_{i,j}\) and \(t_{i,j}\):

\[ s_{i,j} := 2 X s_{i,j} - 2 Y t_{i,j} + \sum_{k=0,l=0}^{k=i,l=j} \left( s_{k,l} s_{i-k,j-l} - t_{k,l} t_{i-k,j-l} \right) + \mathbb{1}_{i=1,j=0} \]

The formula for \(t\) requires knowing which branch of the \(\operatorname{diffabs}(X Y, 0)\) was taken, which turns out has a nice reduction to \(\operatorname{sgn}(XY)\):

\[ t_{i,j} := 2 \operatorname{sgn}(X Y) \left( X t_{i,j} + s_{i,j} Y + \sum_{k=0,l=0}^{k=i,l=j} \left( s_{k,l} t_{i-k,j-l} \right) \right) + \mathbb{1}_{i=0,j=1} \]

\(\mathbb{1}_F\) is the indicator function, \(1\) when \(F\) is true, \(0\) otherwise. For distance estimation, the series of the derivatives are just the derivatives of the series, which can be computed very easily.

The series is only valid in a region that doesn't intersect an axis, at which point the next iteration will fold the region in a way that a series can't represent. Moreover, the series loses accuracy the further from the reference point, so there needs to be a way to check that the series is ok to use. One approach is to iterate points on the boundary of the region using perturbation, and compare the relative error against the same points calculated with the series. Only if all these probe points are accurate, is it safe to try the next iteration.

This usually means that the series will skip at most \(P\) iterations near a central miniship reference of period \(P\). But one can do better, by subdividing the region into parts that are folded in different ways. The parts that are folded the same way as the reference can continue with series approximation, with probe points at the boundary of each part, with the remainder switching to regular perturbation initialized by the series.

It may even be possible to use knighty's techniques like Taylor shift, which is a way to rebase a series to a new reference point (for example, one in the "other side" of the fold) to split the region into two or more separate parts each with their own series approximation. The Horner shift algorithm is not too complicated, and I think it can be extended to bivariate series by shifting along each variable in succession:

// Horner shift code from FLINT2
for (i = n - 2; i >= 0; i--)
  for (j = i; j < n - 1; j++)
    poly[j] += poly[j + 1] * c

Untested Haskell idea for bivariate shift:

shift :: Num v => v -> Map Int v -> Map Int v
shift = undefined -- left as an exercise

shift2 :: Num v => v -> v -> Map (Int, Int) v -> Map (Int, Int) v
shift2 x y
  = twiddle . pull . fmap (shift x) . push
  . twiddle . pull . fmap (shift y) . push

push :: (Ord a, Ord b) => Map (a, b) v -> Map a (Map b v)
push m = M.fromListWith M.union [ (i, M.singleton j e) | ((i,j), e) <- M.assocs m ]

pull :: (Ord a, Ord b) => Map a (Map b v) -> Map (a, b) v
pull m = M.fromList [((i,j),e) | (i, mj) <- M.assocs m, (j, e) <- M.assocs mj ]

twiddle :: (Ord a, Ord b) => Map (a, b) v -> Map (b, a) v
twiddle = M.mapKeys (\(i,j) -> (j,i))

Exciting developments! I hope to release a new version of Kalles Fraktaler 2 + containing at least some of these algorithms soon...