# Fractals/Iterations in the complex plane/demm

This algorithm has 2 versions :

• exterior
• interior

Compare it with version for dynamic plane and Julia set : DEM/J

# Exterior DEM/M

Distance Estimation Method for Mandelbrots set ( DEM/M ) estimates distance from point ${\displaystyle c\,}$ ( in the exterior of Mandelbrot set ) to nearest point in Mandelbrot set.[1][2]

Names

• the (directional) distance estimate formula

It can be used to create B&W images from BOF :[3]

• map 41 on page 84
• map 43 on page 85
• an unnumbered plot on page 188

## types

• analytic DE = DEM = true DE
• Claude DE = fake DE ( FDE): FDE=1/(g*log(2)) when g>0 and undefined when g=0 (i.e., in interior)[4]

## Math formule

Math formula :[5][6] [7]

${\displaystyle distance(c_{0},\sigma M)=\lim _{n\to \infty }2{\frac {|z_{n}|}{|z'_{n}|}}\ln |z_{n}|}$

where :

${\displaystyle z_{n}=f_{c}^{n}(z=0.0)}$
${\displaystyle z_{n}'={\frac {d}{dc}}f_{c}^{n}(z_{0}).}$

is a first derivative of ${\displaystyle f_{c}^{n}(z_{0})}$ with respect to c which can be found by iteration starting with

${\displaystyle z_{0}'={\frac {d}{dc}}f_{c}^{0}(z_{0})=1}$

and then replacing at every consecutive step

${\displaystyle z_{n+1}'={\frac {d}{dc}}f_{c}^{n+1}(z_{0})=2\cdot {}f_{c}^{n}(z)\cdot {\frac {d}{dc}}f_{c}^{n}(z_{0})+1=2\cdot z_{n}\cdot z_{n}'+1.}$

• "you need to change "If mag > 2" to increase the escape radius. Attached is a comparison, showing that the "straps" reduce when the radius is increased, no harm going huge like R = 1e6 or so because once it escapes 2 it grows very fast. The bigger the radius, the more accurate the distance estimate." Claude [8]
• Basically in testing magnitude squared against bailout
• if your bailout is 10^4 then your DE estimates are accurate to 2 decimal figures i.e. if your DE is 0.01 then it's only accurate to +/-0.0001,
• if your bailout is 10^6 then your DE estimates are accurate to 3 decimal figures (David Makin in fractalforums.com)

"For exterior distance estimation, you need a large escape radius, eg 100100. Iteration count limit is arbitrary, with a finite limit some pixels will always be classified as "unknown"." Claude Heiland-Allen[9]

## Comparison

• The Beauty of Fractals
• The Science of Fractal Images, page 198,
• Distance Estimator by Robert P. Munafo[10]

Pseudocode by Claude Heiland-Allen[11]

foreach pixel c
while not escaped and iteration limit not reached
dc := 2 * z * dc + 1
z := z^2 + c
de := 2 * |z| * log(|z|) / |dz|
d := de / pixel_spacing
plot_grey(tanh(d))


Code [12]

      double _Complex m_exterior_distance(int N, double R, double _Complex c)
{
double _Complex z = 0;
double _Complex dc = 0;
for (int n = 0; n < N; ++n)
{
if (cabs(z) > R)
return 2 * z * log(cabs(z)) / dc;
dc = 2 * z * dc + 1;
z = z * z + c;
}
return -1;
}


### Analysis of code from BOF

"The Beauty of Fractals" gives an almost correct computer program for the distance estimation shown in the right image. A possible reason that that method did not gain ground is that the procedure in this program is seriously flawed: The calculation of ${\displaystyle z_{k}}$ is performed (and completed) before the calculation of ${\displaystyle z'_{k}}$, and not after as it ought to be (${\displaystyle z'_{k+1}}$ uses ${\displaystyle z_{k}}$, not ${\displaystyle z_{k+1}}$). For the successive calculation of ${\displaystyle z'_{k}}$, we must know ${\displaystyle f'(z_{k})}$ (which in this case is 2${\displaystyle z_{k}}$). In order to avoid the calculation of ${\displaystyle z_{k}}$ (k = 0, 1, 2, ...) again, this sequence is saved in an array. Using this array, ${\displaystyle z'_{k+1}=2z_{k}z'_{k}+1}$ is calculated up to the last iteration number, and it is stated that overflow can occur. If overflow occurs the point is regarded as belonging to the boundary (the bail-out condition). If overflow does not occur, the calculation of the distance can be performed. Apart from it being untrue that overflow can occur, the method makes use of an unnecessary storing and repetition of the iteration, making it unnecessarily slower and less attractive. The following remark in the book is nor inviting either: "It turns out that the images depend very sensitively on the various choices" (bail-out radius, maximum iteration number, overflow, thickness of the boundary and blow-up factor). Is it this nonsense that has got people to lose all desire for using and generalizing the method? " Gertbuschmann

So :

• in the loop, derivative should be calculated before the new z

"I'm finally generating DEM (Distance Estimate Method) based images that I'm happy with. It turns out that my code had a couple of bugs in it. The new code runs reasonably fast, even with all the extra math to compute distance estimates.

The algorithm in S of F I ("The Science of Fractal Images" uses a sharp cut-off from white to black when pixels get close enough to the set, but I find this makes for jagged looking plots. I've implemented a non-linear color gradient that works pretty well.

For pixels where their DE value is threshold>=DE>=0, I scale the value of DE to 1>=DE>=0, then do the calculation: gradient_value = 1 - (1-DE)^n, ("^n" means to the n power. I wish I could use superscripts!)

"n" is an adjustment factor that lets me change the shape of the gradient curve. The value "gradient_value" determines the color of the pixel. If it's 0, the pixel is colored at the starting color (white, for a B&W plot.) If "gradient_value" is 1, the pixel gets the end color (e.g.

black) For values of n>1, this gives a rapid change in color for pixels that are far from the set, and a slower change in value as DE approaches 0. For values of n<1, it gives a slow color change for pixels far from the set, and rapid change close to the set. For n=1, I get a linear

color gradient. The non-linear gradient lets me use the DE value to anti-alias my plots. By adjusting my threshold value and my adjustment value, I can get good looking results for a variety of images." Duncan C [13]

"All our DE formulas above are only approximations – valid in the limit n→∞, and some also only for point close to the fractal boundary. This becomes very apparent when you start rendering these structures – you will often encounter noise and artifacts.

Multiplying the DE estimates by a number smaller than 1, may be used to reduce noise (this is the Fudge Factor in Fragmentarium).

Another common approach is to oversample – or render images at large sizes and downscale." Mikael Hvidtfeldt Christensen [14]

## Example code

### Two algorithm in two loops

Here is Java function. It computes in one loop both : iteration ${\displaystyle z_{n}\,}$ and derivative ${\displaystyle z'_{n}\,}$. If (on dynamic plane) critical point :

• does not escapes to infinity ( bailouts), then it belongs to interior of Mandelbrot set and has colour maxColor,
• escapes then it maybe in exterior or near boundary. Its colour is "proportional" to distance between the point c and the nearest point in the Mandelbrot set. It uses also colour cycling ( (int)R % maxColor ).
// Java function by Evgeny Demidov from http://www.ibiblio.org/e-notes/MSet/DEstim.htm
// based on code by Robert P. Munafo from http://www.mrob.com/pub/muency/distanceestimator.html
public int iterate(double cr, double ci, double K, double f) {
double Cr,Ci, I=0, R=0,  I2=I*I, R2=R*R, Dr=0,Di=0,D;   int n=0;
if (f == 0){ Cr = cr; Ci = ci;}
else{ Cr = ci; Ci = cr;}
do  {
D = 2*(R*Dr - I*Di) +1;  Di = 2*(R*Di + I*Dr);  Dr = D;
I=(R+R)*I+Ci;  R=R2-I2+Cr;  R2=R*R;  I2=I*I;  n++;
} while ((R2+I2 < 100.) && (n < MaxIt) );
if (n == MaxIt) return maxColor; // interior
else{ // boundary and exterior
R = -K*Math.log(Math.log(R2+I2)*Math.sqrt((R2+I2)/(Dr*Dr+Di*Di))); // compute distance
if (R < 0) R=0;
return (int)R % maxColor; };
}


Here is cpp function. It gives integer index of color as an output.

/*
this function is  from program mandel ver 5.10 by Wolf Jung
see file mndlbrot.cpp
http://www.mndynamics.com/indexp.html

It is checked first (in pixcolor)
that the point escapes before calling this function.
So we do not compute the derivative and the logarithm
when it is not escaping.
This would not be a good idea if most points escape anyway.

*/
int mndlbrot::dist(double a, double b, double x, double y)
{  uint j;
double xp = 1, yp = 0; // zp = xp+yp*i =  1 ; derivative
double nz, nzp;

// iteration
for (j = 1; j <= maxiter; j++)
{  // zp
nz = 2*(x*xp - y*yp);
yp = 2*(x*yp + y*xp);
xp = nz; //zp = 2*z*zp; on the dynamic plane
// if sign is positive it is parameter plane, if negative it is dynamic plane.
if (sign > 0) xp++; //zp = 2*z*zp + 1 ; on the parameter plane
//z = z*z + c;
nz = x*x - y*y + a;
y = 2*x*y + b;
x = nz;
//
nz = x*x + y*y;
nzp = xp*xp + yp*yp;
// 2 conditions for stopping the iterations
if (nzp > 1e40 || nz > bailout) break;
}

// 5 types of points but 3 colors
/*  not escaping, rare
Is it not simply interior ???
This should not be necessary.  I do not know why I included it,
because in this case  pixcolor should not call dist.  If you do not
have pixcolor before,  you should return 0 here. */
if (nz < bailout) return 1; // not escaping, rare
/*  If The Julia set is disconnected and the orbit of z comes close to
z = 0 before escaping,  nzp will be small  */
if (nzp < nz) return 10; // includes escaping through 0

// compute estimated distance  = x
x = log(nz);
x = x*x*nz / (temp[1]*nzp); //4*square of dist/pixelwidth
if (x < 0.04) return 1; // exterior but very close to the boundary
if (x < 0.24) return 9; // exterior but close to the boundary
return 10; //exterior : escaping and not close to the boundary
} //dist


### Two algorithms in one loop

• distance rendering for fractals by ińigo quilez in c++[15]
Function GiveDistance(xy2,eDx,eDy:extended):extended;

begin
result:=2*log2(sqrt(xy2))*sqrt(xy2)/sqrt(sqr(eDx)+sqr(eDy));
end;

//------------------------------------

Function PointIsInCardioid (Cx,Cy:extended):boolean;
//Hugh Allen
// http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip
var DeltaX,DeltaY:extended;
//
PDeltyX,PDeltyY:extended;
//
ZFixedX,ZFixedY:extended;

begin
result:=false;
// cardioid checkig - thx to Hugh Allen
//sprawdzamy Czy punkt  C jest w głównej kardioidzie
//Cardioid in squere :[-0.75,0.4] x [ -0.65,0.65]
if InRange(Cx,-0.75,0.4)and InRange(Cy,-0.65,065) then
begin
// M1= all C for which Fc(z) has  attractive( = stable) fixed point
// znajdyjemy punkt staly z: z=z*z+c
// czyli rozwiazujemy rownanie kwadratowe
// zmiennej zespolonej o wspolczynnikach zespolonych
// Z*Z - Z + C = 0
//Delta:=1-4*a*c; Delta i C sa liczbami zespolonymi
DeltaX:=1-4*Cx;
DeltaY:=-4*Cy;
// Pierwiastek zespolony z delty
CmplxSqrt(DeltaX,DeltaY,PDeltyX,PDeltyY);
// obliczmy punkt staly jeden z dwóch, ten jest prawdopodobnie przycišgajšcy
ZFixedX:=1.0-PDeltyX; //0.5-0.5*PDeltyX;
ZFixedY:=PDeltyY; //-0.5*PDeltyY;
// jesli punkt stały jest przycišgajšcy
// to należy do M1
If (ZfixedX*ZFixedX + ZFixedY*ZFixedY)<1.0
then result:=true;

// ominięcie iteracji M1 przyspiesza 3500 do 1500 msek
end; // if InRange(Cx ...

end;
//------------------------------------
Function PointIsInComponent (Cx,Cy:extended):boolean;
//Hugh Allen
// http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip

var Dx:extended;
begin
result:=false;
// czy punkt C nalezy do koła na lewo od kardioidy
// circle: center = -1.0  and radius 1/4
dx:=Cx+1.0;
if (Dx*Dx+Cy*Cy) < 0.0625
then result:=true;

end;

//------------------------------
Procedure DrawDEM_DazibaoTrueColor;
// draws Mandelbrot set in black and its complement in true color
// see   http://ibiblio.org/e-notes/MSet/DEstim.htm
// by  Evgeny Demidov
//
//http://www.mandelbrot-dazibao.com/Mset/Mset.htm
// translation ( with modification) of Q-Basic program:
//  http://www.mandelbrot-dazibao.com/Mset/Mdb3.bas
//

var iter:integer;
iY,iX:integer;
eCy ,eCx:extended; // C:=eCx + eCy*i
eX,eY:extended;    // Zn:=eX+eY*i
eTempX,eTempY:extended;
eX2,eY2:extended;  //x2:=eX*eX;  y2:=eY*eY;
eXY2:extended;  //  xy2:=x2+y2;
eXY4:extended;
eTempDx,eTempDy:extended;
eDx,eDy:extended; // derivative
distance:extended;
color:TColor;

begin
//compute bitmap
for iY:= iYmin to iYMax do
begin
eCy:=Convert_iY_to_eY(iY);
for iX:= iXmin to iXmax do
begin
eCx:=Convert_iX_to_eX(iX);
If not PointIsInCardioid (eCx,eCy) and not PointIsInComponent(eCx,eCy)
then
begin
//  Z0:=0+0*i
eX:=0;
eY:=0;
eTempX:=0;
eTempY:=0;
//
eX2:=0;
eY2:=0;
eXY2:=0;
//
eDx:=0;
eDy:=0;
eTempDx:=0;
eTempDy:=0;
//
iter:=0;
// iteration of Z ; Z= Z*z +c
while ((iter<IterationMax) and (eXY2<=BailOut2)) do
begin
inc(iter);
//
eTempY:=2*eX*eY + eCy;
eTempX:=eX2-eY2 + eCx;
//
eX2:=eTempX*eTempX;
eY2:=eTempY*eTempY;
//
eTempDx:=1+2*(eX*eDx-eY*eDy);
eTempDy:=2*(eX*eDY+eY*eDx);
//
eXY2:=eX2+eY2;
//
eX:=ETempX;
eY:=eTempY;
//
eDx:=eTempDx;
eDy:=eTempDy;
end; // while

// drawing procedure
if (iter<IterationMax)
then
begin
distance:= GiveDistance(eXY2,eDx,eDy);
color:=Rainbow(1,500,Abs(Round(100*Log10(distance)) mod 500));
with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
begin
B := GetBValue(color);
G := GetGValue(color);
R := GetRValue(color);
//A := 0;
end; // with FirstLine[Y*LineLength+X]

end // if (iter<IterationMax) then
else  with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
begin
B := 0;
G := 0;
R := 0;
//A := 0;
end;
//--- end of drawing procedure ---
end //  If not PointIsInCardioid ... then
else with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
begin
B := 0;
G := 0;
R := 0;
//A := 0;
end;
//--- If not PointIsInCardioid ...
end; // for iX
end; // for iY

end;
//------------------------------------------------------


## Modifications / optimisations

Creating DEM images can be improved by use of :

### equalisation

Distance Estimate Equalisation[19]

Program exrdeeq[20]

• reads DEX and DEY channels
• outputs Y channel
• DE less than 0 becomes 0
• DE greater than 1 becomes 1
• otherwise Y is the index of DE in the sorted array, scaled to the range 0 to 1
exrdeeq in.exr out.exr


## Examples

• B of F map 43 DEM by Duncan Champney[21]
• This plot is centered on -0.7436636773584340, 0.1318632143960234i at a width of about 6.25E-11
• Smily Kaleidoscope by Duncan Champney [22]
height=800
width=800
max_iterations=10000
center_r=-1.74920463346
center_i=-0.000286846601561
r_width=3.19E-10


# Interior distance estimation

Estimates distance from limitly periodic point ${\displaystyle c\,}$ ( in the interior of Mandelbrot set ) to the boundary of Mandelbrot set.

Descriptions :[23]

// https://mathr.co.uk/mandelbrot/book-draft/#interior-distance
// Claude Heiland-Allen
#include <complex.h>
#include <math.h>

double cnorm(double _Complex z)
{
return creal(z) * creal(z) + cimag(z) * cimag(z);
}

double m_interior_distance
(double _Complex z0, double _Complex c, int p)
{
double _Complex z = z0;
double _Complex dz = 1;
double _Complex dzdz = 0;
double _Complex dc = 0;
double _Complex dcdz = 0;
for (int m = 0; m < p; ++m)
{
dcdz = 2 * (z * dcdz + dz * dc);
dc = 2 * z * dc + 1;
dzdz = 2 * (dz * dz + z * dzdz);
dz = 2 * z * dz;
z = z * z + c;
}
return (1 - cnorm(dz))
/ cabs(dcdz + dzdz * dc / (1 - dz));
}

double m_distance(int N, double R, double _Complex c)
{
double _Complex dc = 0;
double _Complex z = 0;
double m = 1.0 / 0.0;
int p = 0;
for (int n = 1; n <= N; ++n)
{
dc = 2 * z * dc + 1;
z = z * z + c;
if (cabs(z) > R)
return 2 * cabs(z) * log(cabs(z)) / cabs(dc);
if (cabs(z) < m)
{
m = cabs(z);
p = n;
double _Complex z0 = m_attractor(z, c, p);
double _Complex w = z0;
double _Complex dw = 1;
for (int k = 0; k < p; ++k)
{
dw = 2 * w * dw;
w = w * w + c;
}
if (cabs(dw) <= 1)
return m_interior_distance(z0, c, p);
}
}
return 0;
}


## Math description

The estimate is given by :

${\displaystyle distance(c,\sigma M)=d(c,z_{0},p)={\frac {1-\mid {{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}\mid ^{2}}{\mid {{\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})+{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0}){\frac {{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})}{1-{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}}}\mid }}}$

where

${\displaystyle p}$ is the period = length of the periodic orbit

${\displaystyle c}$ is the point from parameter plane to be estimated

${\displaystyle f_{c}(z)}$ is the complex quadratic polynomial ${\displaystyle f_{c}(z)=z^{2}+c}$

${\displaystyle f_{c}^{p}(z_{0})}$ is the ${\displaystyle p}$-fold iteration of ${\displaystyle f_{c}(z)}$

${\displaystyle z_{0}}$ is any of the ${\displaystyle p}$ points that make the periodic orbit ( limit cycle ) ${\displaystyle \{z_{0},\dots ,z_{p-1}\}}$

4 derivatives of ${\displaystyle f_{c}^{p}}$ evaluated at ${\displaystyle z_{0},c,p}$ :

${\displaystyle {\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}$

${\displaystyle {\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}$

${\displaystyle {\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})}$

${\displaystyle {\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}$

### First partial derivative with respect to z

It must be computed recursively by applying the chain rule :

Starting points : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=1\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}}$

First iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=2z_{0}\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}}$

Second iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=4z_{0}^{3}+4cz_{0}\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}}$

General rule for p iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}}$

### First partial derivative with respect to c

It must be computed recursively by applying the chain rule :

Starting points : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})=0\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}}$

First iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{0}(z_{0})+1=1\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}}$

Second iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{1}(z_{0})+1=2z_{0}^{2}+2c+1\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}}$

General rule for p iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})+1\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}}$

### Second order partial derivative with respect to z

It must be computed :

• together with :${\displaystyle f_{c}^{n}}$ and ${\displaystyle {\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}$
• recursively by applying the chain rule

Starting points : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=1\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=0\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}}$

First iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=2z_{0}\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2\{({\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0}))^{2}+f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})\}=2\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}}$

Second iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=4z_{0}^{3}+4cz_{0}\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}}$

General rule for p iteration : ${\displaystyle {\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2\{({\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0}))^{2}+f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\}\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}}$

### Second order mixed partial derivative

${\displaystyle {\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2({\frac {\partial }{\partial {c}}}f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})+f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0}))}$

## Algorithm

### Steps

• choose c
• check if critical point for given c is periodic ( interior ) or not ( exterior or near boundary or to big period )
• if not periodic do not use this algorithm ( use exterior version)
• if periodic compute period and periodic point ${\displaystyle z_{0}}$
• using ${\displaystyle c,z_{0},p}$ compute distance using p iteration

## Code

• Java code by Albert Lobo Cusidó[24]
• Haskell code by Claude[25] and image [26]
• processing code[27] and description[28] by tit_toinou

## Problems

There are two practical problems with the interior distance estimate:

• first, we need to find ${\displaystyle z_{0}}$ precisely,
• second, we need to find ${\displaystyle p}$ precisely.

The problem with ${\displaystyle z_{0}}$ is that the convergence to ${\displaystyle z_{0}}$ by iterating ${\displaystyle P_{c}(z)}$ requires, theoretically, an infinite number of operations.

The problem with period is that, sometimes, due to rounding errors, a period is falsely identified to be an integer multiple of the real period (e.g., a period of 86 is detected, while the real period is only 43=86/2). In such case, the distance is overestimated, i.e., the reported radius could contain points outside the Mandelbrot set.

For interior distance estimation, you need the period, then a number (maybe 10 or so should usually be enough) Newton steps to find the limit cycle.

the interior checking code absolutely requires the reference to be at the nucleus of the island (not any of its child discs, and certainly not some random point)

# References

1. A Cheritat wiki-draw: Mandelbrot_set
2. fractalforums : m-set-distance-estimation
3. fractalforums discussion : How are B&W images from "Beauty of Fractals" created?
4. fractalforums.org : m-brot-distance-estimation-versus-claudes-fake-de
5. Milnor (in Dynamics in One Complex Variable, appendix G)
6. Heinz-Otto Peitgen (Editor, Contributor), Dietmar Saupe (Editor, Contributor), Yuval Fisher (Contributor), Michael McGuire (Contributor), Richard F. Voss (Contributor), Michael F. Barnsley (Contributor), Robert L. Devaney (Contributor), Benoit B. Mandelbrot (Contributor) : The Science of Fractal Images. Springer; 1 edition (July 19, 1988), page 199
7. distance rendering for fractals by ińigo quilez
8. fractalforums : mandelbrot-distance-estimation-problem
9. fractalforums.org : list-of-numbers-with-fixed-digit-length-in-the-mandelbrot-set
10. Distance Estimator by Robert P. Munafo
11. math.stackexchange question: how-to-draw-a-mandelbrot-set-with-the-connecting-filaments-visible
12. Mandelbrot book
13. Fractal forum : How are B&W images from "Beauty of Fractals" created?
14. Distance Estimated 3D Fractals (V): The Mandelbulb & Different DE Approximations by Mikael Hvidtfeldt Christensen
15. distance rendering for fractals by ińigo quilez
16. distance estimation by Claude Heiland-Allen
17. distance rendering for fractals by ińigo quilez
18. fractalforums gallery by Pauldelbrot
19. fractalforums.org : histogram-de-coloring
20. Distance-Estimate-Equalisation by Claude Heiland-Allen
21. B of F map 43 DEM by Duncan Champney
22. Smily Kaleidoscope by Duncan Champney
23. Interior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó
24. Interior and exterior distance bounds for the Mandelbrot by Albert Lobo Cusidó
25. Interior coordinates in the Mandelbrot set by Claude
26. Fractal forum : Coloring points inside the Mandelbrot set
27. MandelbrotDE in processing by ttoinou
28. fractalforums : classic-mandelbrot-with-distance-and-gradient-for-coloring