KF's uses floatexp with renormalization of the exponent every arithmetic operation, for zooms deeper than e9800. What Pauldelbrot describes is potentially very much more efficient (renormalization only every 500 iterations, no additional squaring etc), but harder to get right.

KF also uses "rescaled (long) double" at medium zooms (e300-e600, e4900-e9800), by considering an unevaluated product of two doubles, which is a lot faster than floatexp due to hardware support.  Probably the method Pauldelbrot describes is faster there too, but the details seem a bit tricky.

I don't know if Pauldelbrot has ever published an implementation.  MandelMachine might use similar techniques, if you can decipher the assembler.

I suppose any of these methods could be combined with series approximation without any major issues.

Thanks for the link.

I've implemented this method and got it working for at least some cases (without series approximation) - the performance is really good, basically like normal double precision however there are some bugs. I think this methods is very dependent on the reference orbit; like pauldelbrot states I use a regular double precision to record the reference (it is calculated with arbitrary precision with precision relative to window radius). Sometimes, the values recorded for the reference will hit zero (or near). I'm not sure if this is a bug on my end but when a zero is recorded in the double precision reference, the formula breaks down, specifically when the delta_0 value has been scaled such that it rounds to 0. Once this happens, the values of the orbit zero out and the same value is present in the entire image.

Edit: I think the solution to this is to have the reference orbit in floatexp form; that should allow for it to work well - the benefit of this method over the normal floatexp method is that the re-scales only happen every 100-400 iterations and that you don't need to consider the squared term in the perturbation iterations until the delta is more than ~1e-150.

Yes I think the assumptions of the formula break down if the reference is 0.
In those cases you need to do a full iteration with the delta^2 part included.
Probably you can analyze the reference orbit up front, so the main per pixel loop is more like
for (int i = 0; ; ++i)
{
  int iterations_until_next_reference_zero = iterations_until_next_reference_zero_array[i];
  for (int j = 0, k = 0; j < iterations_until_next_reference_zero; ++j, =+k)
  {
    if (k = 500)
    {
      k = 0;
      // rescale delta_scaled back down to be near 1.0, and rescale delta_scaled_0 by the same amount
     }
    delta_scaled_n+1 = delta_scaled_n * twice_reference_orbit_n + delta_scaled_0 // pauldelbrot optimization
  }
  // rescale delta_scaled to be near 1.0, and rescale delta_scaled_0 by the same amount
  // do one iteration assuming twice_reference_orbit is 0
  delta_scaled_n+1 = ldexp(delta_scaled_n * delta_scaled_n, scaling_exponent)  + delta_scaled_0
}
you could probably even omit the addition of delta_scaled_0 if delta_scaled gets big enough (equivalently if delta_scaled_0 gets small enough): if z escapes the minibrot of the reference then it probably escapes everywhere (other non-tuned minibrot interior nearby to a minibrot is vanishingly small once you zoom deep enough).

I'm not sure what happens if the scaling_exponent exceeds the range of double in the last line - probably you have to go full floatexp for one iteration.

and maybe there will be cases where the reference orbit underflows to 0 in double when it's not exactly 0, so perhaps storing (only) those iterations in a floatexp structure could work out ok - duplicating the whole reference orbit as floatexp would seem to be overkill when only a few points will be problematic.

I've now optimized my code a little more and done some benchmarking, comparing a double precision version (with series approximation in double precision), to one using floatexp for series approximation and the extended algorithm for the inner loop. The performance of double vs. this method is very comparable (look at the numbers for iteration). I haven't got distance estimation working at the moment which could add some significant performance issues as it probably will require floatexp. Are the numerical methods for D.E. looking at the escape iteration/radius of the pixels and doing a numerical approximation (and how good do these look compared to the analytical inner loop method)?