Colouring of exterior of Mandelbrot set can be :

non-smooth :
- Boolean/binary Escape Time Method ( bETM/M )
- discrete = Level Set Method = LSM/M = iETM/M
Smooth :
- radial measures
  - Real Escape Time Method( rETM/M )
  - Distance Estimation Method( DEM/M )
  - radius of complex potential = Complex Potential Method ( CPM/M )
angular measures
- - argument of complex potential
- SAC = stripe average coloring
- other
  - "triangle inequality"^[1]^[2]
  - Orbit trap ^[3]^[4]

One can also draw curves :

external rays
equipotential lines ( closed curves - quasi circles)

Similar projects:

Escape time or dwellEdit

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 you will get different result ^[5]

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

How to find the number of iterations required to escape the mandelbrot set ?Edit

Boolean escape timeEdit

This algorithm answers the question: “For which values of c will the Julia fractal, J(c), be line-like and for which dust-like?”^[6]

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)Edit

ASCI graphic : Boolean escape time in text mode

// http://mrl.nyu.edu/~perlin/
main(k){float i,j,r,x,y=-16;while(puts(""),y++<15)for(x
=0;x++<84;putchar(" .:-;!/>)|&IH%*#"[k&15]))for(i=k=r=0;
j=r*r-i*i-2+x/25,i=2*r*i+y/10,j*j+i*i<11&&k++<111;r=j);}

-- Haskell code by Ochronus
-- http://blog.mostof.it/mandelbrot-set-in-ruby-and-haskell/
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]^[8]

Graphic file ( PPM )Edit

Here are various programs for creating pbm file ^[9]

GLSL using shadertoy
short c programs by ed Burke^[10]

CEdit

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";
        char *comment="# ";/* comment should start with # */
        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 ;  */
        const double EscapeRadius=2;
        double ER2=EscapeRadius*EscapeRadius;
        /*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 = dwell bandsEdit

Number of details is proportional to maximal number of iterations
Mandelbrot animation based on a static number of iterations per pixel. Here you can see why offset is sometimes used ( because - color gradient changes : for high MaxIteration disapears.

Here color is proportional to last iteration ( of final_n, final iteration).^[11]

This is also called Level Set Method ( LSM )

$L_{n}=\{c:z_{n}\in T~~{\mbox{and}}~~z_{k}\notin T~~{\mbox{where}}~~k<n\}\,$

CEdit

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 C function whithout explicit complex numbers, only doubles:

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;
}

here a short code with complex numbers:

// https://gitlab.com/adammajewski/mandelbrot_wiki_ACh/blob/master/betm.c
int iterate(double complex C , int iMax)
  {
   int i;
   double complex Z= 0.0; // initial value for iteration Z0
   
   for(i=0;i<iMax;i++)
    {
      Z=Z*Z+C; // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
      if(cabs(Z)>EscapeRadius) break;
    }
   return i; 
 }

C++Edit

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++)
    {
      Z=Z*Z+C; // overloading of operators
      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);
}

A touch more optimised :

// optimised from cpp code by Geek3
inline int fractal(double cReal, double cImg, int max_iters)
// gives last iteration 
{
	double zReal = 0, zImg = 0, zReal2 = 0, zImg2 = 0;
	//iteration zero is always 0^2+0^2, it will never escape
	for (int it = 1; it < max_iters; it++)
	{	//because we have zReal^2 and zImg^2 pre-calculated
		//we can calculate zImg first
		//then we don't need to calculate/store the "old" zReal
		zImg = (2 * zReal * zImg ) + cImg;
		zReal = zReal2 - zImg2 + cReal;

		// calculate next iteration: zReal^2 and zImg^2
		// they are used twice so calculate them once
		zReal2 = zReal * zReal;
		zImg2 = zImg * zImg;
		if (zReal2 + zImg2 > 4.0) return(it);
	}
	return(max_iters);
}

See also :

Rendering The Mandelbrot Set In C++: Part 1 - Saving Images , Video by Superdavo0001

GLSLEdit

shadertoy

JavaEdit

//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 ScriptEdit

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.

Khan AcademyEdit

khan academy mandelbrot-set

Lisp programEdit

Whole Lisp program making ASCII graphic based on code by Frank Buss ^[12]^[13]

; 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 GimpEdit

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-RayEdit

Pov-Ray has a built-in function mandel^[14]

Wolfram MathematicaEdit

Here is code by Paul Nylander

Level Curves of escape time Method = LCM/MEdit

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

edge detection of Level sets.
- Algorithm described in paper by M. Romera et al.^[15] This method is fast and allows looking for high iterations.
boundary trace^[16]
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 target set for Mandelbrot set drawingEdit

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/MEdit

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)\leq 0\}\,$

Division of target set induces decomposition of level sets $L_{n}\,$ into $2^{n+1}\,$ parts ( cells, subsets):

$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})\leq 0\}\,$ which is colored black.

External rays of angles (measured in turns):

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

can be seen as borders of subsets.

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 */    
          };

also GLSL code from Fragmentarium :

#include "2D.frag"
#group Simple Mandelbrot

// maximal number of iterations
uniform int iMax; slider[1,100,1000] // increase iMax
// er2= er^2 wher er= escape radius = bailout
uniform float er2; slider[4.0,1000,10000] // increase er2

// compute color of pixel
 vec3 color(vec2 c) {
   vec2 z = vec2(0.0);  // initial value

 // iteration
  for (int i = 0; i < iMax; i++) {
    z = vec2(z.x*z.x-z.y*z.y,2*z.x*z.y) +  c; // z= z^2+c
    if (dot(z,z)> er2)   // escape test 
      // exterior
      if (z.x>0){ return vec3( 1.0);} // upper part of the target set 
      else return vec3(0.0); //lower part of the target set 
  }
  return vec3(0.0); //interior
}

// zoomasm -- zoom video assembler
// (c) 2019,2020,2021,2022 Claude Heiland-Allen
// SPDX-License-Identifier: AGPL-3.0-only

// recommended KF bailout settings: linear smoothing, custom radius 25
vec3 colour(void)
{
  if (getInterior())
  {
    return vec3(1.0, 0.0, 0.0);
  }
  bool decomp = getT() < 0.5;
  return vec3(decomp ? 0.0 : 1.0);
}

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

nth-decompositionEdit

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).^[18]

Real Escape TimeEdit

Other names of this method/algorithm are :

the fully-renormalized fractional iteration count ( by Linas Vepstas in 1997)^[19]
smooth iteration count for generalized Mandelbrot sets ( by Inigo Quilez in 2016)^[20]
continuous iteration count for the Mandelbrot set
Normalized Iteration Count Algorithm
Continuous coloring
smooth colour gradient
fractional iterations
fractional escape time

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 :^[21]

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

where :

log(log(bailout) can be precalculated

CEdit

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^[22]

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 ^[23]

\nu (z)=n-\log _{2}\log(z_{n})\,

C++Edit

// 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));

javaEdit

 /**
Smooth coloring algorithm
https://gitlab.com/shreyas.siravara/mandelbrot-with-smooth-coloring/blob/master/Mandelbrot.java
Mandelbrot with Smooth Coloring by Shreyas Siravara

*/
                double nsmooth = (iterations + 1 - Math.log(Math.log(Zn.getMagnitude())) / Math.log(ESCAPE_RADIUS));

                double smoothcolor = nsmooth / MAX_ITERATIONS;

                if (iterations < MAX_ITERATIONS) {
                    int rgb = Color.HSBtoRGB((float) (0.99f + 1.9 * smoothcolor), 0.9f, 0.9f);
                    g2d.setColor(new Color(rgb));
                } else {
                    g2d.setColor(Color.black.darker());
                }

MatemathicaEdit

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

$cet=n+log_{2}ln(R)-log_{2}ln|z|$

PythonEdit

Python code using mpmath library^[24]

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

Distance estimation DEM/MEdit

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

Variants :

exterior DEM/M
interior DEM/M

Description

Complex potentialEdit

Description

ReferencesEdit

[1] Mathematics of Divergent Fractals by

[2] ussi harkonen : coloring-techniques

[3] wikipedia : Orbit trap

[4] Mandelbrot Orbit Trap Rendering! Programming How-To Video by DKM101

[5] Java program by Dieter Röß showing result of changing initial point of Mandelbrot iterations

[6] Julia fractals in PostScript by Kees van der Laan, EUROTEX 2012 & 6CM PROCEEDINGS 47

[7] Fractal Benchmark by Erik Wrenholt

[8] 12-minute Mandelbrot: fractals on a 50 year old IBM 1401 mainframe

[9] Computer Language Benchmarks Game

[10] xample-code-from-presentation-ways-of-seeing-julia-sets by ed Burke

[11] Computing the Mandelbrot set by Andrew Williams

[12] LIsp Program by Frank Buss

[13] Mandelbrot Set ASCII art at Bill Clementson's blog

[14] unction from 2.5.11.14 Fractal Patterns at Pov-Ray docs

[15] Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.

[16] ttp://www.metabit.org/~rfigura/figura-fractal/math.html boundary trace by Robert Figura

[17] ttp://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

[18] rot set n-th-decomposition

[19] s.org : Renormalizing the Mandelbrot Escape

[20] I Quilez : mset_smooth

[21] ractalforums : What range/precision for fractional escape counts for Mandelbrot/Julia sets?

[22] Making Mandelbrot Set Movies by Tony Finch

[23] Linas Vepstas. Renormalizing the mandelbrot escape.

[24] th Python library

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

Fractals/Iterations in the complex plane/MandelbrotSetExterior

Contents

Escape time or dwellEdit

How to find the number of iterations required to escape the mandelbrot set ?Edit

Boolean escape timeEdit

ASCI graphic ( on screen)Edit

Graphic file ( PPM )Edit

CEdit

Integer escape time = LSM/M = dwell bandsEdit

CEdit

C++Edit

GLSLEdit

JavaEdit

Java ScriptEdit

Khan AcademyEdit

Lisp programEdit

MathMap plugin for GimpEdit

Pov-RayEdit

Wolfram MathematicaEdit

Level Curves of escape time Method = LCM/MEdit

Decomposition of target set for Mandelbrot set drawingEdit

Binary decomposition of LSM/MEdit

nth-decompositionEdit

Real Escape TimeEdit

CEdit

C++Edit

javaEdit

MatemathicaEdit

PythonEdit

Distance estimation DEM/MEdit

Complex potentialEdit

See alsoEdit

ReferencesEdit