# Fractals/Iterations in the complex plane/Mandelbrot set

This book shows how to code different algorithms for drawing parameter plane [1] ( Mandelbrot set [2] ) for complex quadratic polynomial [3].

One can find different types of points / sets on parameter plane [4]

# Exterior of Mandelbrot set

Colouring of exterior of Mandelbrot set can be :

• boolean ( bETM/M )
• discrete ( LSM/M = iETM/M)
• Smooth  :
• Real Escape Time ( rETM/M )
• Distance estimation ( DEM/M )
• Complex potential ( CPM/M )
• "triangle inequality"[5]

One can also draw curves :

• external rays
• equipotential lines ( closed curves - quasi circles)

## Escape time

Here for given point c on parameter plane one checks how critical point $z = 0.0$ behaves on dynamical plane under forward iteration. If you change initial point ypu will get different result [6]

To draw given plane one needs to check/scan (all) its points. See here for more details ( optimisation) Read definitions first.

### Boolean escape time

Here complex plane consists of 2 sets : Mandelbrot set $M\,$ and its complement $M^c\,$ :

$\mathbb{C}= M \cup M^c$

#### ASCI graphic ( on screen)

ASCI graphic : Boolean escape time in text mode

-- Haskell code by Ochronus
import Data.Complex
mandelbrot a = iterate (\z -> z^2 + a) a !! 500
main = mapM_ putStrLn [[if magnitude (mandelbrot (x :+ y)) < 2 then '*' else ' '
| x <- [-2, -1.9685 .. 0.5]]
| y <- [1, 0.95 .. -1]]


; common lisp
(loop for y from -1.5 to 1.5 by 0.05 do
(loop for x from -2.5 to 0.5 by 0.025 do
(let* ((c (complex x y)) ; parameter
(z (complex 0 0))
(iMax 20) ; maximal number of iterations
(i 0)) ; iteration number

(loop   while (< i iMax ) do
(setq z (+ (* z z) c)) ; iteration
(incf i)
(when (> (abs z) 2) (return i)))
; color of pixel
(if (= i iMax) (princ (code-char 42)) ; inside M
(princ (code-char 32))))) ; outside M
(format t "~%")) ; new line


Comparison programs in various languages [7]

#### Graphic file ( PPM )

Here are various programs for creating pbm file [8]

##### C

This is complete code of C one file program.

• It makes a ppm file which consists an image. To see the file (image) use external application ( graphic viewer).
• Program consists of 3 loops:
• iY and iX, which are used to scan rectangle area of parameter plane
• iterations.

For each point of screen (iX,iY) it's complex value is computed c=cx+cy*i.

For each point c is computed iterations of critical point $z_0 = z_{cr} = 0 \,$

It uses some speed_improvement. Instead of checking :

sqrt(Zx2+Zy2)<ER


it checks :

(Zx2+Zy2)<ER2 // ER2 = ER*ER


It gives the same result but is faster.

 /*
c program:
--------------------------------
1. draws Mandelbrot set for Fc(z)=z*z +c
using Mandelbrot algorithm ( boolean escape time )
-------------------------------
2. technique of creating ppm file is  based on the code of Claudio Rocchini
http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
create 24 bit color graphic file ,  portable pixmap file = PPM
see http://en.wikipedia.org/wiki/Portable_pixmap
to see the file use external application ( graphic viewer)
*/
#include <stdio.h>
#include <math.h>
int main()
{
/* screen ( integer) coordinate */
int iX,iY;
const int iXmax = 800;
const int iYmax = 800;
/* world ( double) coordinate = parameter plane*/
double Cx,Cy;
const double CxMin=-2.5;
const double CxMax=1.5;
const double CyMin=-2.0;
const double CyMax=2.0;
/* */
double PixelWidth=(CxMax-CxMin)/iXmax;
double PixelHeight=(CyMax-CyMin)/iYmax;
/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
const int MaxColorComponentValue=255;
FILE * fp;
char *filename="new1.ppm";
static unsigned char color[3];
/* Z=Zx+Zy*i  ;   Z0 = 0 */
double Zx, Zy;
double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
/*  */
int Iteration;
const int IterationMax=200;
/* bail-out value , radius of circle ;  */
/*create new file,give it a name and open it in binary mode  */
fp= fopen(filename,"wb"); /* b -  binary mode */
/*write ASCII header to the file*/
fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
/* compute and write image data bytes to the file*/
for(iY=0;iY<iYmax;iY++)
{
Cy=CyMin + iY*PixelHeight;
if (fabs(Cy)< PixelHeight/2) Cy=0.0; /* Main antenna */
for(iX=0;iX<iXmax;iX++)
{
Cx=CxMin + iX*PixelWidth;
/* initial value of orbit = critical point Z= 0 */
Zx=0.0;
Zy=0.0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
/* */
for (Iteration=0;Iteration<IterationMax && ((Zx2+Zy2)<ER2);Iteration++)
{
Zy=2*Zx*Zy + Cy;
Zx=Zx2-Zy2 +Cx;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
/* compute  pixel color (24 bit = 3 bytes) */
if (Iteration==IterationMax)
{ /*  interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
else
{ /* exterior of Mandelbrot set = white */
color[0]=255; /* Red*/
color[1]=255;  /* Green */
color[2]=255;/* Blue */
};
/*write color to the file*/
fwrite(color,1,3,fp);
}
}
fclose(fp);
return 0;
}


### Integer escape time = LSM/M

Here color is proportional to last iteration.

This is also called Level Set Method ( LSM )

$L_n= \{ c : z_n \in T ~~\mbox{and} ~~ z_k \notin T ~~\mbox{where}~~ k

#### C

LSM/M image with full code in C

Difference between Mandelbrot algorithm and LSM/M is in only in part instruction, which computes pixel color of exterior of Mandelbrot set. In LSM/M is :

 if (Iteration==IterationMax)
{ /* interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
/* exterior of Mandelbrot set = LSM */
else if ((Iteration%2)==0)
{ /* even number = black */
color[0]=0; /* Red */
color[1]=0; /* Green */
color[2]=0; /* Blue */
}
else
{/* odd number =  white */
color[0]=255; /* Red */
color[1]=255; /* Green */
color[2]=255; /* Blue */
};


Here is faster C function which can be used instead of above C++ function:

int GiveEscapeTime(double C_x, double C_y, int iMax, double _ER2)
{
int i;
double Zx, Zy;
double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */

Zx=0.0; /* initial value of orbit = critical point Z= 0 */
Zy=0.0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;

for (i=0;i<iMax && ((Zx2+Zy2)<_ER2);i++)
{
Zy=2*Zx*Zy + C_y;
Zx=Zx2-Zy2 +C_x;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
return i;
}


#### C++

Here is C++ function which can be used to draw LSM/M :

 int iterate_mandel(complex C , int imax, int bailout)
{
int i;
std::complex Z(0,0); // initial value for iteration Z0

for(i=0;i<=imax-1;i++)
{
if(abs(Z)>bailout)break;
}
return i;
}


I think that it can't be coded simpler (it looks better than pseudocode), but it can be coded in other way which can be executed faster .

Here is faster code :

// based on cpp code by Geek3
inline int fractal(double cx, double cy, int max_iters)
// gives last iteration
{
double zx = 0, zy = 0;
if (zx * zx + zy * zy > 4) return(0); // it=0
for (int it = 1; it < max_iters; it++)
{       double zx_old = zx;
zx = zx * zx - zy * zy;
zy = 2 * zx_old * zy;
zx += cx;
zy += cy;
if (zx * zx + zy * zy > 4.0) return(it);
}
return(max_iters);
}


#### Java

//Java code by Josef Jelinek
// http://java.rubikscube.info/

int mandel(double cx, double cy) {
double zx = 0.0, zy = 0.0;
double zx2 = 0.0, zy2 = 0.0;
int iter = 0;
while (iter < iterMax && zx2 + zy2 < 4.0) {
zy = 2.0 * zx * zy + cy;
zx = zx2 - zy2 + cx;
zx2 = zx * zx;
zy2 = zy * zy;
iter++;
}
return iter;
}


#### Java Script

Here is JavaScript function which does not give last iteration but LastIteration modulo maxCol. It makes colour cycling ( if maxCol < maxIt ).

function iterate(Cr,Ci) {
// JavaScript function by Evgeny Demidov
// http://www.ibiblio.org/e-notes/html5/fractals/mandelbrot.htm
var I=0, R=0,  I2=0, R2=0,   n=0;
if (R2+I2 > max) return 0;
do  {  I=(R+R)*I+Ci;  R=R2-I2+Cr;  R2=R*R;  I2=I*I;  n++;
} while ((R2+I2 < max) && (n < maxIt) );
if (n == maxIt) return maxCol;  else return n % maxCol;
}


Above functions do not use explicit definition of complex number.

#### Lisp program

Whole Lisp program making ASCI graphic based on code by Frank Buss [9] [10]

; common lisp
(loop for y from -1.5 to 1.5 by 0.1 do
(loop for x from -2.5 to 0.5 by 0.04 do
(let* ((i 0)
(z (complex x y))
(c z))
(loop while (< (abs
(setq z (+ (* z z) c)))
2)
while (< (incf i) 32))
(princ (code-char (+ i 32))))) ; ASCII chars <= 32 contains non-printing characters
(format t "~%"))


#### MathMap plugin for Gimp

filter mandelbrot (gradient coloration)
c=ri:(xy/xy:[X,X]*1.5-xy:[0.5,0]);
z=ri:[0,0]; # initial value z0 = 0
# iteration of z
iter=0;
while abs(z)<2 && iter<31
do
z=z*z+c;  # z(n+1) = fc(zn)
iter=iter+1
end;
coloration(iter/32) # color of pixel
end


#### Pov-Ray

Pov-Ray has a built-in function mandel[11]

#### Matemathica

Here is code by Paul Nylander

### Level Curves of escape time Method = LCM/M

edge detection of Level sets
Lemniscates of Mandelbrot set

Lemniscates are boundaries of Level Sets of escape time ( LSM/M ). They can be drawn using 2 methods:

• edge detection of Level sets.
• Algorithm described in paper by M. Romera et al.[12] This method is fast and allows looking for high iterations.
• boundary trace[13]
• drawing curves $L_n(T)=\{c: abs(z_n)=ER \}\,$, see explanation and source code. This method is very complicated for iterations > 5.

### Decomposition of exterior of Mandelbrot set

Decomposition is modification of escape time algorithm.

The target set is divided into parts (2 or more). Very large escape radius is used, for example ER = 12.

#### Binary decomposition of LSM/M

binary decomposition: image with full code in C

Here target set $T\,$ on dynamic plane is divided into 2 parts (binary decomposition = 2-decomposition ):

• upper half ( white) $T_u = \{ z : |z| > ER ~~\mbox{and}~~ Im (z) > 0\}\,$
• lower half (black) $T_l = \{ z : |z| > ER ~~\mbox{and}~~ Im (z) \le 0 \}\,$

Division of target set induces decomposition of level sets $L_n\,$ into $2^{n+1}\,$ parts:

• $L_{n,u} =\{ c: |z_n| > ER ~~\mbox{and}~~ Im (z_n) > 0 \}\,$ which is colored white,
• $L_{n,l} = \{ c : |z_n| > ER ~~\mbox{and}~~ Im (z_n) \le 0 \}\,$ which is colored black.

External rays of angles (measured in turns):

$angle = (k / 2^n ) ~~\mbox{mod }~1\,$

can be seen.

Difference between binary decomposition algorithm and Mandel or LSM/M is in only in part of instruction , which computes pixel color of exterior of Mandelbrot set. In binary decomposition is :

 if (Iteration==IterationMax)
{ /* interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
/* exterior of Mandelbrot set = LSM */
else if (Zy>0)
{
color[0]=0; /* Red */
color[1]=0; /* Green */
color[2]=0; /* Blue */
}
else
{
color[0]=255; /* Red */
color[1]=255; /* Green */
color[2]=255; /* Blue */
};


Point c is plotting white or black if imaginary value of last iteration ( Zy) is positive or negative.[14]

#### nth-decomposition

This method is extension of binary decomposition.

The target set T = { z : |zn| > R } with a very large escape radius ( for example R = 12 ) is divided into more then 2 parts ( for example 8).[15]

### Real Escape Time

Other names of this method/algorithm are :

• Normalized Iteration Count Algorithm
• Continuous coloring
• fractional iterations

Here color of exterior of Mandelbrot set is proportional not to Last Iteration ( which is integer number) but to real number :

$\nu(z) = \lim_{i \to \infty} (i - \log_2 \log_2 |z_i|)\,$

Other methods and speedups

Colouring formula in Ultrafractal :[16]

smooth iter = iter + 1 + ( log(log(bailout)-log(log(cabs(z))) )/log(2)


where :

• log(log(bailout) can be precaclulated

#### C

To use log2 function add :

#include <math.h>

at the beginning of program.

if (Iteration==IterationMax)
{ /*  interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
/* exterior of Mandelbrot set  */
else GiveRainbowColor((double)(Iteration- log2(log2(sqrt(Zx2+Zy2))))/IterationMax,color);


where :

• Zx2 = Zx*Zx
• Zy2 = Zy*Zy

Here is another version by Tony Finch[17]

while (n++ < max &&
x2+y2 < inf) {
y = 2*x*y + b;
x = x2-y2 + a;
y2 = y*y;
x2 = x*x;
}
nu = n - log(log(x2+y2)/2)/ log(2);


based on equation [18]

$\nu(z) = n - \log_2 \log (z_n)\,$

#### C++

// based on cpp code by Geek3 from http://en.wikibooks.org/wiki/File:Mandelbrot_set_rainbow_colors.png
sqrxy = x * x + y * y;
double m = LastIteration + 1.5 - log2(log2(sqrxy));


#### Matemathica

Here is code by Paul Nylander. It uses different formula :

$cet = n + log_2ln(R) - log_2ln|z|$

#### Python

Python code using mpmath library[19]

def mandelbrot(z):
c = z
for i in xrange(ITERATIONS):
zprev = z
z = z*z + c
return ctx.exp(1j*(i + 1 - ctx.log(ctx.log(abs(z)))/ctx.log(2)))
return 0


## Distance estimation DEM/M

Exterior DEM/M
simple boundary with DEM/M
Boundary with DEM/M and Sobel filter

Estimates distance from point $c\,$ ( in the exterior of Mandelbrot set ) to nearest point in Mandelbrot set.

Math formula : [20] [21]

$distance(z_0,\sigma M) =\lim_{n \to \infty} 2 \frac{|z_n|}{|z'_n|}\ln|z_n|$

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

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

Creating DEM images can be improved by use of :

### 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 $z_k$ is performed (and completed) before the calculation of $z'_k$, and not after as it ought to be ($z'_{k+1}$ uses $z_k$, not $z_{k+1}$). For the successive calculation of $z'_k$, we must know $f'(z_k)$ (which in this case is 2$z_k$). In order to avoid the calculation of $z_k$ (k = 0, 1, 2, ...) again, this sequence is saved in an array. Using this array, $z'_{k+1} = 2z_kz'_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

### Example code

#### Two algorithm in two loops

It is easy to create such loop, but it is slow method

#### Two algorithms in one loop

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

// this function is  from program mandel ver 5.3 by Wolf Jung
// http://www.iram.rwth-aachen.de/~jung/indexp.html
// it needs escape time before

int mndlbrot::dist(double a, double b, double x, double y)
{  uint j; double xp = 1, yp = 0, nz, nzp; //zp = 1
for (j = 1; j <= maxiter; j++)
{  nz = 2*(x*xp - y*yp); yp = 2*(x*yp + y*xp); xp = nz; //zp = 2*z*zp;
if (sign > 0) xp++; //zp = 2*z*zp + 1
nz = x*x - y*y + a; y = 2*x*y + b; x = nz; //z = z*z + c;
nz = x*x + y*y; nzp = xp*xp + yp*yp;
if (nzp > 1e60 || nz > 1e60) break;
}
if (nz < bailout) return 1; //not escaping
if (nzp < 0.01) return 11; //escaping through 0
x = log(nz); x = x*x*nz / (temp[1]*nzp); //square of dist/pixelwidth
if (x < 0.08) return 1;
if (x < 0.15) return 9;
if (x < 1.0) return 3;
return 11;
} //dist


Here is Java function. It computes in one loop both : iteration $z_n\,$ and derivative $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; };
}


## Complex potential

Uniformization of complement of Mandelbrot set

Complex potential is a complex number, so it has a real (real potential) and an imaginary part(external angle). One can take its curl, divergence or its absolute value. So on one image one can use more than one variable to color image.[26]

### Real potential = CPM/M

In Fractint :

potential =  log(modulus)/2^iterations


One can use real potential to:

• smooth (continuous) coloring[27]
• discrete coloring ( level sets of potential)
• 3D view

Here is Delphi function which gives level of potential :

 Function GiveLevelOfPotential(potential:extended):integer;
var r:extended;
begin
r:= log2(abs(potential));
result:=ceil(r);
end;


### External angle and external ( parameter) ray

First find angle of last iteration. It is easy to compute and shows some external rays as a borders of level sets.

Then one go futher :

• Wolf Jung explanation[28]
• Tomoki Kawahira paper [29]

#### Tests

The Wolf Jung test : The external parameter rays for angles (in turns)

• 1/7 (period 3)
• 321685687669320/2251799813685247 (period 51)
• 321685687669322/2251799813685247 ( period 51 )

Angles differ by about $10^{-15}$, but the landing points of the corresponding parameter rays are about 0.035 apart. [30]

The test by G. Pastor and Miguel Romera The external parameter rays for angles (in turns)

• 6871947673/34359738367 ( period 35 )
• 9162596898/34359738367 ( period 35 )

the central babies Mandelbrot sets of the cauliflowers located at -0.153756141 + 1.030383223i

(not that 34359738367 = 2^35 - 1)

# Boundary of Mandelbrot set

Boundary of Mandelbrot set as an image of unit circle under $\Psi_M\,$

## Mini Mandelbrot sets

" a saddle-node (parabolic) periodic point in a complex dynamical system often admits homoclinic points, and in the case that these homoclinic points are nondegenerate, this is accompanied by the existence of infinitely many baby Mandelbrot sets converging to the saddle-node parameter value in the corresponding parameter plane." Devaney [31]

## Drawing boundaries

Methods used to draw boundary of Mandelbrot set :[32]

## Topological models of Mandelbrot set

Topological model of Mandelbrot set = Mandelbrot cactus
Arcs up to period 12
• Thurston model ( abstract Mandelbrot set )[37]
• Cactus model [40]( go to the image, where the src code and description of algorithm is )
• Fake Mandelbrot Set = M-set without hairs, filaments and primitive hyberbolic componennts[41]

## Computing Misiurewicz points of complex quadratic mapping

"... we do not know how to compute (...) Misiurewicz parameters (with high (pre)periods) for the family of quadratic rational maps. One might need to and a non-rigorous method to and Misiurewicz parameter in a reasonable time like Biham-Wenzel's method." HIROYUKI INOU [42]

### First method

Misiurewicz points [43] are special boundary points.

Define polynomial in Maxima CAS :

P(n):=if n=0 then 0 else P(n-1)^2+c;


Define a Maxima CAS function whose roots are Misiurewicz points, and find them.

M(preperiod,period):=allroots(%i*P(preperiod+period)-%i*P(preperiod));


Examples of use :

(%i6) M(2,1);
(%o6) [c=-2.0,c=0.0]
(%i7) M(2,2);
(%o7) [c=-1.0*%i,c=%i,c=-2.0,c=-1.0,c=0.0]


### Second Method

" factorizing the polynomials that determine Misiurewicz points. I believe that you should start with

  ( f^(p+k-1) (c) + f^(k-1) (c) ) / c


This should already have exact preperiod k , but the period is any divisor of p . So it should be factorized further for the periods.

Example: For preperiod k = 1 and period p = 2 we have

  c^3 + 2c^2 + c + 2 .


This is factorized as

(c + 2)*(c^2 + 1)


for periods 1 and 2 . I guess that these factors appear exactly once and that there are no other factors, but I do not know."[44]

# Interior of Mandelbrot set - hyperbolic components

## The Ljapunov Exponent

• image and description by janthor [45]
• image by Anders Sandberg [46]

## Interior distance estimation

Interior distance estimation

Description of method [47]

## absolute value of the orbit

# Hypercomputing the Mandelbrot Set? by Petrus H. Potgieter February 1, 2008
n=1000; # For an nxn grid
m=50; # Number of iterations
c=meshgrid(linspace(-2,2,n))\ # Set up grid
+i*meshgrid(linspace(2,-2,n))’;
x=zeros(n,n); # Initial value on grid
for i=1:m
x=x.^2+c; # Iterate the mapping
endfor
imagesc(min(abs(x),2.1)) # Plot monochrome, absolute
# value of 2.1 is escape


### internal level sets

Color of point  :

• is proportional to the value of z is at final iteration.
• shows internal level sets of periodic attractors.

### bof60

Image of bof60 in on page 60 in the book "the Beauty Of Fractals".Description of the method described on page 63 of bof. It is used only for interior points of the Mandelbrot set.

Color of point is proportional to :

• the smallest distance of its orbit from origin[48][49]
• the smallest value z gets during iteration [50]
• illuminating the closest approach the iterates of the origin (critical point) make to the origin inside the set

Level sets of distance are sets of points with the same distance[51]

if (Iteration==IterationMax)
/* interior of Mandelbrot set = color is proportional to modulus of last iteration */
else { /* exterior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}

• fragment of code : fractint.cfrm from Gnofract4d [52]
bof60 {
init:
float mag_of_closest_point = 1e100
loop:
float zmag = |z|
if zmag < mag_of_closest_point
mag_of_closest_point = zmag
endif
final:
#index = sqrt(mag_of_closest_point) * 75.0/256.0
}


### bof61

3D view of bof61
2D view of BOF61

This is the method described in the book "The Beauty of Fractals" on page 63, but the image in on page 61.

Color of point is proportional to :

• the time it takes z to reach its smallest value
• iterate of the critical point makes the closest approach
• Index (c) is the iteration when point of the orbit was closest to the origin. Since there may be more than one, index(c) is the least such.

This algorithms shows borders of domains with the same index(c)[53] [54].

Fragment of code : fractint.cfrm from Gnofract4d [55]

bof61 {
init:
int current_index = -1 ; -1 to match Fractint's notion of iter count
int index_of_closest_point = -1
float mag_of_closest_point = 1e100
loop:
current_index = current_index + 1
float zmag = |z|
if zmag < mag_of_closest_point
index_of_closest_point = current_index
mag_of_closest_point = zmag
endif
final:
#index = index_of_closest_point /256.0
}


Cpp function

//       function is based on function mclosetime
//       from mbrot.cpp
//       from program mandel by Wolf Jung
//      http:www.iram.rwth-aachen.de/~jung/indexp.html
////8 = iterate =  bof61
// bailout2=4
int mclosetime(std::complex<double> C , int iter_max,  int bailout2)
{ int j, cln = 0;
double x = C.real(), y = C.imag(), u, cld = 4;
for (j = 0; j <= iter_max; j++)
{ u = x*x + y*y;
if (u > bailout2) return j;
if (u < cld) {cld = u;cln = j; }

u = x*x - y*y + C.real();
y = 2*x*y + C.imag();
x = u;
}
return iter_max + cln % 15; //iterate, bof61

}


It can be used :

// compute escape time
int last_iter= mclosetime(C,iter_max,bailout2);
//  drawing code */
if (last_iter>=iter_max) { putpixel(ix,iy,last_iter - iter_max);} // interior
else putpixel(ix,iy,WHITE); // exterior


## Period of hyperbolic components

period of hyperbolic components

Algorithms :

• direct period detection from iterations
• the spider algorithm
• "quick and dirty" algorithm ( check if $abs(z_n ) < eps$ then colour c-point with colour n )

## Multiplier map

definition

### Internal rays

Internal and external rays

When $radius\,$ varies and $angle\,$ is constant then $c\,$ goes along internal ray. It is used as a path inside Mandelbrot set

/* find c in component of Mandelbrot set
uses complex type so #include <complex.h> and -lm
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html

*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
//0 <= InternalRay<= 1
//0 <= InternalAngleInTurns <=1
double t = InternalAngleInTurns *2*M_PI; // from turns to radians
double Cx, Cy; /* C = Cx+Cy*i */
switch ( period ) {
case 1: // main cardioid
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t) - 1.0;
break;
// for each period  there are 2^(period-1) roots.
default: // safe values
Cx = 0.0;
Cy = 0.0;
break; }

return Cx+ Cy*I;
}

// draws points to memmory array data
int DrawInternalRay(double InternalAngleInTurns , unsigned int period, int iMax, unsigned char data[])
{

complex double c;
int i; // number of point to draw

for(i=0;i<=iMax;++i){
DrawPoint(c,data);
}

return 0;
}


Example : internal ray of angle =1/6 of main cardioid.

Internal angle :

$angle = 1/6 \,$

$0 \le radius \le 1 \,$

Point of internal radius of unit circle :

$w = radius * e^{i * angle}\,$

Map point $w$ to parameter plane :

$c = \frac{w}{2} - \frac{w^2}{4} \,$

For $epsilon = 0 \,$ this is equation for main cardioid.

### Internal curve

When $radius\,$ is constant varies and $angle\,$ varies then $c\,$ goes along internal curve.

/* find c in component of Mandelbrot set
uses complex type so #include <complex.h> and -lm
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html

*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
//0 <= InternalRay<= 1
//0 <= InternalAngleInTurns <=1
double t = InternalAngleInTurns *2*M_PI; // from turns to radians
double Cx, Cy; /* C = Cx+Cy*i */
switch ( period ) {
case 1: // main cardioid
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t) - 1.0;
break;
// for each period  there are 2^(period-1) roots.
default: // safe values
Cx = 0.0;
Cy = 0.0;
break; }

return Cx+ Cy*I;
}

// draws points to memmory array data
int DrawInternalCurve(double InternalRadius , unsigned int period,  int iMax, unsigned char data[])
{

complex double c;
double InternalAngle; // in turns = from 0.0 to 1.0
double AngleStep;
int i;
// int iMax =100;

AngleStep = 1.0/iMax;

for(i=0;i<=iMax;++i){
InternalAngle = i * AngleStep;
DrawPoint(c,data);
}

return 0;
}


### Boundaries of hyperbolic components

Boundaries of hyperbolic components for periods 1-6 as a solutions of boundary equations

Methods of drawing boundaries:

• solving boundary equations :
• using method by Brown,Stephenson,Jung . It works only up to period 7-8 [56]
• using resultants [57]
• parametrisation of boundary with Newton method near centers of components[58] [59]. This methods needs centers, so it has some limitations,
• Boundary scanning  : detecting the edges after detecting period by checking every pixel. This method is slow but allows zooming. Draws "all" components
• Boundary tracing for given c value. Draws single component.
• Fake Mandelbrot set by Anne M. Burns : draws main cardioid and all its descendants. Do not draw mini Mandelbrot sets. [60]

"... to draw the boundaries of hyperbolic components using Newton's method. That is, take a point in the hyperbolic component that you are interested in (where there is an attracting cycle), and then find a curve along which the modulus of the multiplier tends to one. Then you will have found an indifferent parameter. Now you can similarly change the argument of the multiplier, again using Newton's method, and trace this curve. Some care is required near "cusps". " Lasse Rempe-Gillen[61]

#### solving boundary equations

##### System of 2 equations defining boundaries of period $n\,$ hyperbolic components
• first defines periodic point,
• second defines indifferent orbit ( multiplier of periodic point is equal to one ).

$\begin{cases} F(n,z,c) = z \\ abs(\lambda (z) )=1 \end{cases}$

Because stability index $abs(\lambda)\,$ is equal to radius of point of unit circle $abs(w)\,$:

$abs(\lambda) = abs(w)\,$

so one can change second equation to form [62] :

$\lambda = w \,$

It gives system of equations :

$\begin{cases} F(n,z,c) = z \\ \lambda = w \end{cases}$

It can be used for :

• drawing components ( boundaries, internal rays )
• finding indifferent parameters ( parabolic or for Siegel discs )

Above system of 2 equations has 3 variables : $z,c,w \,$ ( $n\,$ is constant and multiplier $\lambda\,$ is a function of $z,c \,$). One have to remove 1 variable to be able to solve it.

Boundaries are closed curves : cardioids or circles. One can parametrize points of boundaries with angle $t\,$ ( here measured in turns from 0 to 1 ).

After evaluation of $w = l(t) \,$ one can put it into above system, and get a system of 2 equations with 2 variables $z,c\,$.

Now it can be solved

For periods:

• 1-3 it can be solved with symbolical methods and give implicit ( boundary equation) $b_p(w,c)=0 \,$ and explicit function (inverse multiplier map)  : $c=\gamma_p(w)\,$
• 1-2 it is easy to solve [63]
• 3 it can be solve using "elementary algebra" ( Stephenson )
• >3 it can't be reduced to explicitly function but
• can be reduced to implicit solution ( boundary equation) $b_p(w,c)=0 \,$ and solved numerically
• can be solved numerically using Newton method for solving system of nonlinear equations

###### Solving system of equation for period 1

Here is Maxima code :

(%i4) p:1;
(%o4) 1
(%i5) e1:F(p,z,c)=z;
(%o5) z^2+c=z
(%i6) e2:m(p)=w;
(%o6) 2*z=w
(%i8) s:eliminate ([e1,e2], [z]);
(%o8) [w^2-2*w+4*c]
(%i12) s:solve([s[1]], [c]);
(%o12) [c=-(w^2-2*w)/4]
(%i13) define (m1(w),rhs(s[1]));
(%o13) m1(w):=-(w^2-2*w)/4


Second equation contains only one variable, one can eliminate this variable. Because boundary equation is simple so it is easy to get explicit solution

m1(w):=-(w^2-2*w)/4

###### Solving system of equation for period 2

Here is Maxima code using to_poly_solve package by Barton Willis:

(%i4) p:2;
(%o4) 2
(%i5) e1:F(p,z,c)=z;
(%o5) (z^2+c)^2+c=z
(%i6) e2:m(p)=w;
(%o6) 4*z*(z^2+c)=w
(%i7) e1:F(p,z,c)=z;
(%o7) (z^2+c)^2+c=z
to_poly_solve([e1, e2], [z, c]);
(%o10) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac
(%o11) [[z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4]]
(%i12) s:to_poly_solve([e1, e2], [z, c]);
(%o12) [[z=sqrt(w)/2,c=-(w-2*sqrt(w))/4],[z=-sqrt(w)/2,c=-(w+2*sqrt(w))/4],[z=(sqrt(1-w)-1)/2,c=(w-4)/4],[z=-(sqrt(1-w)+1)/2,c=(w-4)/4]]
(%i14) rhs(s[4][2]);
(%o14) (w-4)/4
(%i16) define (m2 (w),rhs(s[4][2]));
(%o16) m2(w):=(w-4)/4


explicit solution :

m2(w):=(w-4)/4

###### Solving system of equation for period 3

For period 3 ( and higher) previous method give no results (Maxima code) :

(%i14) p:3;
e1:z=F(p,z,c);
e2:m(p)=w;
to_poly_solve([e1, e2], [z, c]);
(%o14) 3
(%o15) z=((z^2+c)^2+c)^2+c
(%o16) 8*z*(z^2+c)*((z^2+c)^2+c)=w
(%i17)
(%o17) C:/PROGRA~1/MAXIMA~1.1/share/maxima/5.16.1/share/contrib/topoly_solver.mac
algsys' cannot solve - system too complicated.
#0: to_poly_solve(e=[z = ((z^2+c)^2+c)^2+c,8*z*(z^2+c)*((z^2+c)^2+c) = w],vars=[z,c])
-- an error.  To debug this try debugmode(true);


I use code by Robert P. Munafo[64] which is based on paper of Wolf Jung.

##### Boundary equation

The result of solving above system with respect to $c\,$ is boundary equation,

$b_p(w,c)=0 \,$

where $b_p(w,c)\,$ is boundary polynomial.

It defines exact coordinates of hyperbolic componenets for given period $p \,$.

It is implicit equation.

Period 1

One can easly compute boundary point c

$c = c_x + c_y*i$

of period 1 hyperbolic component ( main cardioid) for given internal angle ( rotation number) t using this code by Wolf Jung[65]

t *= (2*PI); // from turns to radians
cx = 0.5*cos(t) - 0.25*cos(2*t);
cy = 0.5*sin(t) - 0.25*sin(2*t);



Period 2

t *= (2*PI);  // from turns to radians
cx = 0.25*cos(t) - 1.0;
cy = 0.25*sin(t);


##### Solving boundary equations

Solving boundary equations for various angles gives list of boundary points.

### Centers of components

Methods of finding centers :

• all centers for given period $p\,$
• solving $F_p(0,c) = 0 \,$ or $G_p(0,c) = 0\,$
• showing results in graphical form
• using numerical methods[66]
• finding all roots of polynomial
• Myrberg's method [67][68]
• interval arithmetics and a divide and conquere strategy[69]
• one center $c_p\,$ for given period $p\,$ near given point $c\,$
• some steps of Newton method for solving $F_p(0,c_p) = 0 \,$

#### System of equations defining centers

Centers of hyperbolic components for period 10
Centers of hyperbolic components for period 1-10 computed using irreducible polynomials

Center of hyperbolic component is a point $c\,$ of parameter plain, for which periodic z-cycle is superattracting. It gives a system of 2 equations defining centers of period $p\,$ hyperbolic components :

• first defines periodic point $z_0\,$,
• second makes $z_0\,$ point superattracting.

$\begin{cases} F_p(z_0,c) = z_0 \\ \frac{df_c^{(p)}(z_0)}{dz}=0 \end{cases}$

see definitions for details.

#### Solving system of equations

Second equation can be solved when critical point $z_{cr} = 0 \,$ is in this cycle :

$z_0 = z_{cr} = 0 \,$

To solve system put $z_0\,$ into first equation.

##### Equation defining centers

One gets equation :

$F_p(0,c) = 0 \,$

Centers of components are roots of above equation.

Because $z = 0 \,$ one can remove z from these equations :

• for period 1 : z^2+c=z and z=0 so $c=0$
• for period 2 : (z^2+c)^2 +c =z and z=0 so $c^2 + c =0$
• for period 3 : ((z^2+c)^2 +c)^2 +c = z and z=0 so $(c^2 + c)^2 +c =0$

Here is Maxima functions for computing above functions :

P[n]:=if n=1 then c else P[n-1]^2+c;

##### Reduced equation defining centers

Set of roots of above equation contains centers for period p and its divisors. It is because :

$F_p(0,c) = \prod_{m|p} G_m(0,c) \,$

where :

• $G\,$ is irreducible divisors [70] of $F_p \,$ [71][1]
• capital Pi notation of iterated multiplication is used
• $m|p\,$ means here : for all divisors of $p\,$ ( from 1 to p ). See table of divisors.

For example :

$F_1 = c = G_1 \,$

$F_2 = c*(c+1) = G_1* G_2 \,$

$F_3 = c*(c^3+2*c^2+c+1) = G_1* G_3 \,$

$F_4 = c*(c+1)*(c^6+3*c^5+3*c^4+3*c^3+2*c^2+1) = G_1* G_2* G_4 \,$

So one can find irreducible polynomials using :

$G_p(0,c) = \frac{F_p(0,c) }{ \prod_{m|p,m

Here is Maxima function for computing $G \,$ :

GiveG[p]:=
block(
[f:divisors(p),
t:1], /* t is temporary variable = product of Gn for (divisors of p) other than p */
f:delete(p,f), /* delete p from list of divisors */
if p=1
then return(Fn(p,0,c)),
for i in f do t:t*GiveG[i],
g: Fn(p,0,c)/t,
return(ratsimp(g))
)$ Here is a table with degree of $F\,$ and $G\,$ for periods 1-10 and precision needed to compute roots of these functions. period $p \,$ $\deg(F_p(0,c)) \,$ $\deg(G_p(0,c)) \,$ $fpprec \,$ 1 1 1 16 2 2 1 16 3 4 3 16 4 8 6 16 5 16 15 16 6 32 27 16 7 64 63 32 8 128 120 64 9 256 252 128 10 512 495 300 11 1024 1023 600 Here is a table of greatest coefficients. period $p \,$ $log10(gc_p)\,$ $gc_p \,$ 1 0 1 2 0 1 3 1 2 4 1 3 5 3 116 6 4 5892 7 11 17999433372 8 21 106250977507865123996 9 44 16793767022807396063419059642469257036652437 10 86 86283684091087378792197424215901018542592662739248420412325877158692469321558575676264 11 179 307954157519416310480198744646044941023074672212201592336666825190665221680585174032224052483643672228833833882969097257874618885560816937172481027939807172551469349507844122611544 Precision can be estimated as bigger than size of binary form of greatest coefficient $log_2(gc)\,$ : $fpprec > log_2(gc)\,$  period $p \,$ $\deg(G_p(0,c)) \,$ $\log_2(gc_p) \,$ $fpprec \,$ 1 1 1 16 2 1 1 16 3 3 2 16 4 6 2 16 5 15 7 16 6 27 13 16 7 63 34 32 8 120 67 64 9 252 144 128 10 495 286 300 11 1023 596 600 Here is Maxima code for $gc\,$ period_Max:11; /* ----------------- definitions -------------------------------*/ /* basic function = monic and centered complex quadratic polynomial http://en.wikipedia.org/wiki/Complex_quadratic_polynomial */ f(z,c):=z*z+c$
/* iterated function */
fn(n, z, c) :=
if n=1         then f(z,c)
else f(fn(n-1, z, c),c) $/* roots of Fn are periodic point of fn function */ Fn(n,z,c):=fn(n, z, c)-z$
/* gives irreducible divisors of polynomial Fn[p,z=0,c] */
GiveG[p]:=
block(
[f:divisors(p),t:1],
g,
f:delete(p,f),
if p=1
then return(Fn(p,0,c)),
for i in f do t:t*GiveG[i],
g: Fn(p,0,c)/t,
return(ratsimp(g))
)\$
/* degree of function */
GiveDegree(_function,_var):=hipow(expand(_function),_var);
log10(x) := log(x)/log(10);
/* ------------------------------*/
file_output_append:true; /* to append new string to file */
grind:true;
for period:1 thru period_Max step 1 do
(
g[period]:GiveG[period], /* function g */
d[period]:GiveDegree(g[period],c), /* degree of function g */
cf[period]:makelist(coeff(g[period],c,d[period]-i),i,0,d[period]), /* list of coefficients */
cf_max[period]:apply(max, cf[period]), /* max coefficient */
disp(cf_max[period]," ",ceiling(log10(cf_max[period]))),
s:concat("period:",string(period),"  cf_max:",cf_max[period]),
stringout("max_coeff.txt",s)/* save output to file as Maxima expressions */
);


#### Graphical methods for finding centers

All these methods shows centers for period n and its divisors.

Animation showing colors proportional to absolute value of iterations 1-20
Absolute value of first iteration
Color shows in which quadrant 5-th iteration lands
##### color is proportional to magnitude of zn
• color is proportional to magnitude of zn[72]
• Parameter plane is scanned for points c for which orbit of critical point vanishes [73]
• youtube video : Mandelbrot Oscillations [74]
##### color shows in which quadrant zn lands
ninth-decomposition. Image and src code

This is radial nth-decomposition of exterior of Mandelbrot set ( compare it with n-th decomposition of LSM/M )

4 colors are used because there are 4 quadrants :

• re(z_n) > 0 and im(z_n) > 0 ( 1-st quadrant )
• re(z_n) < 0 and im(z_n) > 0 ( 2-nd quadrant )
• re(z_n) < 0 and im(z_n) < 0 ( 3-rd quadrant )
• re(z_n) > 0 and im(z_n) < 0 (4-th quadrant ).

"... when the four colors are meeting somewhere, you have a zero of q_n(c), i.e., a center of period dividing n. In addition, the light or dark color shows if c is in M." ( Wolf Jung )

Here is fragment of c code for computing 8-bit color for point c = cx + cy * i  :

/* initial value of orbit = critical point Z= 0 */
Zx=0.0;
Zy=0.0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
/* the same number of iterations for all points !!! */
for (Iteration=0;Iteration<IterationMax; Iteration++)
{
Zy=2*Zx*Zy + Cy;
Zx=Zx2-Zy2 +Cx;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
/* compute  pixel color (8 bit = 1 byte) */
if ((Zx>0) && (Zy>0)) color[0]=0;   /* 1-st quadrant */
if ((Zx<0) && (Zy>0)) color[0]=10; /* 2-nd quadrant */
if ((Zx<0) && (Zy<0)) color[0]=20; /* 3-rd quadrant */
if ((Zx>0) && (Zy<0)) color[0]=30; /* 4-th quadrant */


One can also find cpp code in function quadrantqn from class mndlbrot in mndlbrot.cpp file from source code of program Mandel by Wolf Jung[75]

The differences from standard loop are :

• there is no bailout test (for example for max value of abs(zn) ). It means that every point has the same number of iterations. For large n overflow is possible. Also one can't use test Iteration==IterationMax
• Iteration limit is relatively small ( start for example IterationMax = 8 )

# Video

One can make videos using :

• going along some paths on parameter plane ( for example internal and external rays )
• zoom into parameter plane[78][79][80] using automatic determination of Iteration Max number[81]
• changing coloring scheme ( for example color cycling - Fractint)
• changing some parameters of algorithm, for example :
• maximal iteration of escape time algorithm
• bailout value [82]

# Speed improvements - optimisation

## Symmetry

The Mandelbrot set is symmetric with respect to the x-axis in the plane :

colour(x,y) = colour(x,-y)


its intersection with the x-axis ( real slice of Mandelbrot set ) is an interval :

<-2 ; 1/4>


It can be used to speed up computations ( up to 2-times )[83]

## Bailout test

$\sqrt{Z_x^2 + Z_y^2} < ER$


compute ER2 once and check :

$Z_x^2 + Z_y^2 < ER^2$


It gives the same result and is faster.

## Period detection

"When rendering a Mandelbrot or Julia set, the most time-consuming parts of the image are the “black areas”. In these areas, the iterations never escape to infinity, so every pixel must be iterated to the maximum limit. Therefore, much time can be saved by using an algorithm to detect these areas in advance, so that they can be immediately coloured black, rather than rendering them in the normal way, pixel by pixel. FractalNet uses a recursive algorithm to split the image up into blocks, and tests each block to see whether it lies inside a “black area”. In this way, large areas of the image can be quickly coloured black, often saving a lot of rendering time. ... (some) blocks were detected as “black areas” and coloured black immediately, without having to be rendered. (Other) blocks were rendered in the normal way, pixel by pixel." Michael Hogg [84]

 // cpp code by Geek3
// http://commons.wikimedia.org/wiki/File:Mandelbrot_set_rainbow_colors.png
bool outcircle(double center_x, double center_y, double r, double x, double y)
{ // checks if (x,y) is outside the circle around (center_x,center_y) with radius r
x -= center_x;
y -= center_y;
if (x * x + y * y > r * r)
return(true);
return(false);

// skip values we know they are inside
if ((outcircle(-0.11, 0.0, 0.63, x0, y0) || x0 > 0.1)
&& outcircle(-1.0, 0.0, 0.25, x0, y0)
&& outcircle(-0.125, 0.744, 0.092, x0, y0)
&& outcircle(-1.308, 0.0, 0.058, x0, y0)
&& outcircle(0.0, 0.25, 0.35, x0, y0))
{
// code for iteration
}
`

### Cardioid and period-2 checking

One way to improve calculations is to find out beforehand whether the given point lies within the cardioid or in the period-2 bulb. Before passing the complex value through the escape time algorithm, first check if:

$(x+1)^2 + y^2 < \frac{1}{16}$

to determine if the point lies within the period-2 bulb and

$q = \left(x - \frac{1}{4}\right)^2 + y^2$
$q \left(q + \left(x - \frac{1}{4}\right)\right) < \frac{1}{4}y^2.$

to determine if the point lies inside the main cardioid.

### Periodicity checking

Most points inside the Mandelbrot set oscillate within a fixed orbit. There could be anything from ten to tens of thousands of points in between, but if an orbit ever reaches a point where it has been before then it will continually follow this path, will never tend towards infinity and hence is in the Mandelbrot set. This Mandelbrot processor includes periodicity checking (and p-2 bulb/cardioid checking) for a great speed up during deep zoom animations with a high maximum iteration value.

# References

1. parameter plane in wikipedia
2. Mandelbrot set in wikipedia
3. complex quadratic polynomial in wikipedia
4. reenigne blog : mandelbrot-set-taxonomy
5. Mathematics of Divergent Fractals by
6. Java program by Dieter Röß showing result of changing initial point of Mandelbrot iterations
7. Fractal Benchmark by Erik Wrenholt
8. The Computer Language Benchmarks Game
9. LIsp Program by Frank Buss
10. Mandelbrot Set ASCII art at Bill Clementson's blog
11. mandel function from 2.5.11.14 Fractal Patterns at Pov-Ray docs
12. Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
13. http://www.metabit.org/~rfigura/figura-fractal/math.html boundary trace by Robert Figura
14. http://web.archive.org/20010415125044/www.geocities.com/CapeCanaveral/Launchpad/5113/fr27.htm%7C An open letter to Dr. Meech from Joyce Haslam in FRACTAL REPORT 27
15. mandelbrot set n-th-decomposition
16. fractalforums : What range/precision for fractional escape counts for Mandelbrot/Julia sets?
17. Making Mandelbrot Set Movies by Tony Finch
18. Linas Vepstas. Renormalizing the mandelbrot escape.
19. mpmath Python library
20. 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
21. distance rendering for fractals by ińigo quilez
22. fractalforums discussion : How are B&W images from "Beauty of Fractals" created?
23. distance estimation by Claude Heiland-Allen
24. distance rendering for fractals by ińigo quilez
25. fractalforums gallery by Pauldelbrot
26. | The Mandelbrot Function by John J. G. Savard
27. Smooth colouring is the key to the Mandelbrot set by Tony Finch
28. Wolf Jung explanation
29. An algorithm to draw external rays of the Mandelbrot set by Tomoki Kawahira
30. Wolf Jung's test for precision of drawing parameter rays
31. Devaney In Global Analysis of Dynamical Systems, ed.: H. Broer, B. Krauskopf, G. Vegter. IOP Publishing (2001), 329-338 or Homoclinic Points in Complex Dynamical Systems
32. mathoverflow : Parametrization of the boundary of the Mandelbrot set
33. Boundary Scanning by Robert P. Munafo, 1993 Feb 3.
35. Efficient Boundary Tracking Through Sampling by Alex Chen , Todd Wittman , Alexander Tartakovsky , and Andrea Bertozzi
36. http://www.robertnz.net/cx.htm contour integration by Robert Davies
37. Combinatorics in the Mandelbrot Set - Lavaurs Algorithm by
38. Abstract Mandelbrot tree by claudius maximus
39. The Mandelbrot cactus by Robert L. Devaney
40. Burns A M : : Plotting the escape- An animation of parabolic bifurrcation in the mandelbrot set. Mathematics Magazine: Volume 75, Number 2, Pages: 104-116 page 104
41. VISUALIZATION OF THE BIFURCATION LOCUS OF CUBICPOLYNOMIAL FAMILY by HIROYUKI INOU
42. MIsiurewicz point in wikipedia
43. Wolf Jung
44. Ljapunov Exponent and mandelbrot set by janthor
45. Image by Anders Sandberg
46. Fractint : Misc. Options and algorithms
47. Java™ Number Cruncher: The Java Programmer's Guide to Numerical Computing By Ronald Mak
48. Firefly Application Help by Terry W. Gintz
49. Fractint doc by Noel Giffin
50. gnofract4d
51. Fractint doc by Noel Giffin
52. A Series of spiral bays in the Mandelbrot set by Patrick Hahn
53. gnofract4d
54. Mandelbrot set Components : image with descriptions and references
55. Young Hee Geum, Kevin G. Hare : Groebner Basis, Resultants and the generalized Mandelbrot Set - methods for period 2 and higher, using factorization for finding irreducible polynomials
56. Image : parametrisation of boundary with Newton method near centers of componentswith src code
57. Mark McClure "Bifurcation sets and critical curves" - Mathematica in Education and Research, Volume 11, issue 1 (2006).
58. Burns A M : Plotting the Escape: An Animation of Parabolic Bifurcations in the Mandelbrot Set. Mathematics Magazine, Vol. 75, No. 2 (Apr., 2002), pp. 104-116
59. mathoverflow : what-are-good-methods-for-detecting-parabolic-components-and-siegel-disk-componen by Lasse Rempe-Gillen
60. Newsgroups: gmane.comp.mathematics.maxima.general Subject: system of equations Date: 2008-08-11 21:44:39 GMT
61. Thayer Watkins : The Structure of the Mandelbrot Set
62. Brown Method by Robert P. Munafo
63. Mandel: software for real and complex dynamics by Wolf Jung
64. Piers Lawrence and Rob Corless : Mandelbrot Polynomials and Matrices. 1 ; 048 ; 575 roots of period 21
65. [http://www.iec.csic.es/~gerardo/publica/Alvarez98a.pdf Determination of Mandelbrot Set's Hyperbolic Component Centres by G Alvarez, M Romera, G Pastor and F Montoya. Chaos, Solitons & Fractals 1998, Vol 9 No 12 pp 1997-2005]
66. Myrberg, P J, Iteration der reelen polynome zweiten GradesIII, Ann. Acad. Sci. Fennicae, A. I no. 348, 1964.
67. evaldraw script that can compute 1 million root by knighty
68. irreducible divisors at wikipedia
69. V Dolotin , A Morozow : On the shapes of elementary domains or why Mandelbrot set is made from almost ideal circles ?
70. [http://www.cerebralmastication.com/2009/02/mandelbrot-set-in-r/ Mandelbrot Set in R by J.D. Long]
71. Fractal Stream help file
72. youtube video : Mandelbrot Oscillations
73. Multiplatform cpp program Mandel by Wolf Jung
74. Formation of Escape-Time Fractals By Christopher Olah
75. [http://www.dhushara.com/DarkHeart/DarkHeart.htm#_3__The_Dark Exploding the Dark Heart of Chaos by Chris King March-Dec 2009]
76. Really Deep Fractal Zoom Movie – Much Faster by Bruce Dawson
77. Making Mandelbrot Set Movies by Tony Finch
78. MLbrot by Daniel de Rauglaud
79. Discussion : A way to determine the ideal number of maximum iterations for an arbitrary zoom level in a Mandelbrot fractal
80. Gif image by jgabase : a wormhole effect on your fractals by changing the bailout dynamicaly
81. How to use symetry of set
82. FractalNet by Michael Hogg
83. ASCII graphic