Fractals/Iterations in the complex plane/MandelbrotSetExterior

Colouring of exterior of Mandelbrot set can be : non-smooth = Escape Time = dwell Boolean/binary Escape Time Method ( bETM/M ) discrete = Level Set Method = LSM/M = integer ETM = 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"Mathematics of Divergent Fractals by jussi harkonen : coloring-techniques Orbit trap wikipedia : Orbit trapMandelbrot Orbit Trap Rendering! Programming How-To Video by DKM101 One can also draw curves : external rays equipotential lines ( closed curves - quasi circles) Similar projects: Mandelbrot Notebook by Claude Heiland-Allen mandelbrot-book mandelbrot-book/book : one file c programs , images and description fork of the mandelbrot-book and book in pdf format different drawing techniques and algorithms by Arnaud Cheritat Mandelbrot set and example C programs ( not by A Cheritat) Linas Vepstas' Art Gallery: original site fork at gitlab Escape time or dwell Here for given point c on parameter plane one checks how critical point behaves on dynamical plane under forward iteration. If you change initial point you will get different result Java program by Dieter Röß showing result of changing initial point of Mandelbrot iterations 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 ? math.stackexchange question: is-there-an-equation-for-the-number-of-iterations-required-to-escape-the-mandelb math.stackexchange question: a-way-to-determine-the-ideal-number-of-maximum-iterations-for-an-arbitrary-zoom? Boolean escape time This algorithm answers the question: “For which values of c will the Julia fractal, J(c), be line-like and for which dust-like?” Julia fractals in PostScript by Kees van der Laan, EUROTEX 2012 & 6CM PROCEEDINGS 47 Here complex plane consists of 2 sets : Mandelbrot set and its complement : ASCI graphic ( on screen) right|thumb|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 Fractal Benchmark by Erik Wrenholt12-minute Mandelbrot: fractals on a 50 year old IBM 1401 mainframe Graphic file ( PPM ) Here are various programs for creating pbm file Computer Language Benchmarks Game GLSL using shadertoy short c programs by ed Burkeexample-code-from-presentation-ways-of-seeing-julia-sets by ed Burke 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 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 bands Here color is proportional to last iteration ( of final_n, final iteration).Computing the Mandelbrot set by Andrew Williams This is also called Level Set Method ( LSM ) C thumb|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 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++ 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 C++ High-Precision Perturbation CLI Renderer thumb|Mandelbrot set fragment using perturbation theory. 4.31e-51 thumb|Mandelbrot set fragment using perturbation theory. 2.4e-57 thumb|Mandelbrot set fragment using perturbation theory. 1.01e-55 thumb|Mandelbrot set fragment using perturbation theory. 7.17e-36 thumb|Mandelbrot set fragment using perturbation theory. 1.53e-62 thumb|Mandelbrot set fragment using perturbation theory. 2.84e-82 thumb|Mandelbrot set fragment using perturbation theory. 1.15e-119 Here is a complete, highly optimized C++ implementation for high-precision rendering of the Mandelbrot set using perturbation theory, 8x8 Super-Sampling Anti-Aliasing (SSAA), and OpenMP multi-threading. The calculation can be performed significantly faster by using bilinear approximation. Created by User:Aokoroko. #include <atomic> #include <cmath> #include <cstdint> #include <cstdio> #include <vector> #include <algorithm> #include <mpfr.h> #include <omp.h> using std::vector; using std::min; const char * CENTER_RE = "-1.99999543561201124623198345433951143502785679245726844745821388800402678499411681518036306219179273434395557574279985918047221291197081186140687781560831995"; const char * CENTER_IM = "-0.00000000000000000000000026198152173811047783694060060607013913873144250985383083459221663448338433592617272786772587281530484110756597337683912309313885172"; const char * VIEW_SIZE = "1.15e-119"; const int WIDTH = 2160; const int HEIGHT = 2160; const int AA = 8; const int MAX_ITER = 50000; const double ESCAPE_RADIUS_SQUARED = 50000.0; const int PALETTE_FRAME = 200; const char * OUTPUT_FILE = "Mandelbrot Set Image 112.bmp"; const mpfr_prec_t PRECISION_BITS = 1000; const int REF_SIZE = MAX_ITER + 200; struct Complex { double re; double im; }; #pragma pack(push, 1) struct BMPHeader { uint16_t type{0x4D42}; uint32_t size{0}; uint32_t reserved{0}; uint32_t offBits{54}; uint32_t structSize{40}; int32_t width{0}; int32_t height{0}; uint16_t planes{1}; uint16_t bitCount{24}; uint32_t compression{0}; uint32_t sizeImage{0}; int32_t xPixelsPerMeter{2834}; int32_t yPixelsPerMeter{2834}; uint32_t colorsUsed{0}; uint32_t colorsImportant{0}; }; #pragma pack(pop) int main() { const double startTime = omp_get_wtime(); const long rawWidth = static_cast<long>(WIDTH) * AA; const long rawHeight = static_cast<long>(HEIGHT) * AA; mpfr_t centerRe, centerIm, zReMp, zImMp, tmp1, tmp2, viewSizeMp; mpfr_inits2(PRECISION_BITS, centerRe, centerIm, zReMp, zImMp, tmp1, tmp2, viewSizeMp, static_cast<mpfr_ptr>(nullptr)); mpfr_set_str(centerRe, CENTER_RE, 10, MPFR_RNDN); mpfr_set_str(centerIm, CENTER_IM, 10, MPFR_RNDN); mpfr_set_str(viewSizeMp, VIEW_SIZE, 10, MPFR_RNDN); const double sampleStep = mpfr_get_d(viewSizeMp, MPFR_RNDN) / rawWidth; vector<Complex> referenceOrbit; referenceOrbit.reserve(REF_SIZE); mpfr_set_ui(zReMp, 0, MPFR_RNDN); mpfr_set_ui(zImMp, 0, MPFR_RNDN); for (int iter = 0; iter < REF_SIZE - 1; ++iter) { Complex z{mpfr_get_d(zReMp, MPFR_RNDN), mpfr_get_d(zImMp, MPFR_RNDN)}; referenceOrbit.push_back(z); if (z.re * z.re + z.im * z.im > ESCAPE_RADIUS_SQUARED) { break; } mpfr_mul(tmp1, zReMp, zImMp, MPFR_RNDN); mpfr_sqr(tmp2, zReMp, MPFR_RNDN); mpfr_sqr(zReMp, zImMp, MPFR_RNDN); mpfr_sub(zReMp, tmp2, zReMp, MPFR_RNDN); mpfr_add(zReMp, zReMp, centerRe, MPFR_RNDN); mpfr_mul_2ui(tmp1, tmp1, 1, MPFR_RNDN); mpfr_add(zImMp, tmp1, centerIm, MPFR_RNDN); } const int referenceLength = static_cast<int>(referenceOrbit.size()); mpfr_clears(centerRe, centerIm, zReMp, zImMp, tmp1, tmp2, viewSizeMp, static_cast<mpfr_ptr>(nullptr)); std::fprintf(stderr, "Reference orbit: %d points\n", referenceLength); std::fprintf(stderr, "Precomputing skip100 matrices...\n"); vector<Complex> coeff_A(REF_SIZE, {1.0, 0.0}); vector<Complex> coeff_B(REF_SIZE, {0.0, 0.0}); vector<double> rad_R(REF_SIZE, 2.0); vector<double> aS_squared(REF_SIZE, 0.0); for (int i = 0; i < referenceLength; ++i) { double r2 = referenceOrbit[i].re * referenceOrbit[i].re + referenceOrbit[i].im * referenceOrbit[i].im; aS_squared[i] = (r2 < ESCAPE_RADIUS_SQUARED) ? r2 : 0.0; } const int loop_limit = min(static_cast<int>(MAX_ITER), referenceLength - 105); #pragma omp parallel for for (int i = 0; i < loop_limit; ++i) { double min_r2 = ESCAPE_RADIUS_SQUARED; for (int k = 0; k < 100; ++k) { if (i + k >= referenceLength) break; if (aS_squared[i + k] < min_r2) min_r2 = aS_squared[i + k]; } rad_R[i] = std::sqrt(min_r2); for (int k = 0; k < 100; ++k) { if (i + k >= referenceLength) break; double r_re = referenceOrbit[i + k].re; double r_im = referenceOrbit[i + k].im; double next_A_re = 2.0 * (r_re * coeff_A[i].re - r_im * coeff_A[i].im); double next_A_im = 2.0 * (r_re * coeff_A[i].im + r_im * coeff_A[i].re); double next_B_re = 2.0 * (r_re * coeff_B[i].re - r_im * coeff_B[i].im) + 1.0; double next_B_im = 2.0 * (r_re * coeff_B[i].im + r_im * coeff_B[i].re); coeff_A[i].re = next_A_re; coeff_A[i].im = next_A_im; coeff_B[i].re = next_B_re; coeff_B[i].im = next_B_im; } } const double PI = 3.14159265358979323846; uint8_t palette[256][3]; for (int i = 0; i < 255; ++i) { palette[i][0] = static_cast<uint8_t>(std::lround(127.0 + 127.0 * std::cos(2.0 * PI * i / 255.0))); palette[i][1] = static_cast<uint8_t>(std::lround(127.0 + 127.0 * std::sin(2.0 * PI * i / 255.0))); palette[i][2] = palette[i][1]; } palette[255][0] = 255; palette[255][1] = 255; palette[255][2] = 255; const int rowBytes = (WIDTH * 3 + 3) & ~3; vector<uint8_t> image(static_cast<size_t>(rowBytes) * HEIGHT, 0); std::atomic<int> completedRows{0}; const Complex * reference = referenceOrbit.data(); #pragma omp parallel for schedule(dynamic) for (int y = 0; y < HEIGHT; ++y) { uint8_t * row = image.data() + static_cast<size_t>(y) * rowBytes; for (int x = 0; x < WIDTH; ++x) { unsigned blueSum = 0; unsigned greenSum = 0; unsigned redSum = 0; for (int sampleY = 0; sampleY < AA; ++sampleY) { const double deltaCIm = (static_cast<long>(y) * AA + sampleY - rawHeight / 2) * sampleStep; for (int sampleX = 0; sampleX < AA; ++sampleX) { const double deltaCRe = (static_cast<long>(x) * AA + sampleX - rawWidth / 2) * sampleStep; double deltaRe = 0.0; double deltaIm = 0.0; double zRe = 0.0; double zIm = 0.0; int referenceIndex = 0; int iter = 0; while (iter < MAX_ITER) { if (zRe * zRe + zIm * zIm >= ESCAPE_RADIUS_SQUARED) { break; } double eps_abs2 = deltaRe * deltaRe + deltaIm * deltaIm; double limit_r2 = 1e-60 * rad_R[referenceIndex] * rad_R[referenceIndex]; if (eps_abs2 < limit_r2 && (referenceIndex + 100 < loop_limit) && (iter + 100 < MAX_ITER)) { double backup_deltaRe = deltaRe; double backup_deltaIm = deltaIm; int backup_refIdx = referenceIndex; int backup_iter = iter; double next_eps_re = (coeff_A[referenceIndex].re * deltaRe - coeff_A[referenceIndex].im * deltaIm) + (coeff_B[referenceIndex].re * deltaCRe - coeff_B[referenceIndex].im * deltaCIm); double next_eps_im = (coeff_A[referenceIndex].re * deltaIm + coeff_A[referenceIndex].im * deltaRe) + (coeff_B[referenceIndex].re * deltaCIm + coeff_B[referenceIndex].im * deltaCRe); deltaRe = next_eps_re; deltaIm = next_eps_im; referenceIndex += 100; iter += 100; zRe = reference[referenceIndex].re + deltaRe; zIm = reference[referenceIndex].im + deltaIm; if (zRe * zRe + zIm * zIm >= ESCAPE_RADIUS_SQUARED) { deltaRe = backup_deltaRe; deltaIm = backup_deltaIm; referenceIndex = backup_refIdx; iter = backup_iter; } else { continue; } } const double a = 2.0 * reference[referenceIndex].re + deltaRe; const double b = 2.0 * reference[referenceIndex].im + deltaIm; const double nextDeltaRe = a * deltaRe - b * deltaIm + deltaCRe; deltaIm = a * deltaIm + b * deltaRe + deltaCIm; deltaRe = nextDeltaRe; ++referenceIndex; ++iter; zRe = reference[referenceIndex].re + deltaRe; zIm = reference[referenceIndex].im + deltaIm; if (zRe * zRe + zIm * zIm < deltaRe * deltaRe + deltaIm * deltaIm || referenceIndex >= loop_limit) { deltaRe = zRe; deltaIm = zIm; referenceIndex = 0; } } const int remaining = MAX_ITER - iter; const uint8_t colorIndex = (remaining == 0) ? 255 : static_cast<uint8_t>(remaining % 254); const int paletteIndex = (colorIndex == 255) ? 255 : (colorIndex - PALETTE_FRAME + 255) % 255; blueSum += palette[paletteIndex][0]; greenSum += palette[paletteIndex][1]; redSum += palette[paletteIndex][2]; } } const int samples = AA * AA; row[x * 3 + 0] = static_cast<uint8_t>(blueSum / samples); row[x * 3 + 1] = static_cast<uint8_t>(greenSum / samples); row[x * 3 + 2] = static_cast<uint8_t>(redSum / samples); } const int done = ++completedRows; if (done % 50 == 0 || done == HEIGHT) { std::fprintf(stderr, "\rProgress: %d/%d rows (%.1f%%)", done, HEIGHT, 100.0 * done / HEIGHT); } } BMPHeader header; header.width = WIDTH; header.height = HEIGHT; header.sizeImage = static_cast<uint32_t>(image.size()); header.size = header.sizeImage + 54; FILE * file = std::fopen(OUTPUT_FILE, "wb"); if (!file) { std::perror(OUTPUT_FILE); return 1; } std::fwrite(&header, sizeof(header), 1, file); std::fwrite(image.data(), 1, image.size(), file); std::fclose(file); std::fprintf(stderr, "\nDone! Saved to %s in %.2f seconds.\n", OUTPUT_FILE, omp_get_wtime() - startTime); return 0; } See the full project and compiled releases on the official Mandelbrot CLI GitHub repository. C++ Perturbation Theory Optimization showcased in the global algorithm repository on Rosetta Code. GLSL shadertoy 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. Khan Academy khan academy mandelbrot-set Lisp program Whole Lisp program making ASCII graphic based on code by Frank Buss LIsp Program by Frank Buss Mandelbrot Set ASCII art at Bill Clementson's blog ; 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 mandelmandel function from 2.5.11.14 Fractal Patterns at Pov-Ray docs Wolfram Mathematica Here is code by Paul Nylander Level Curves of escape time Method = LCM/M right|thumb|edge detection of Level sets thumb|right|Lemniscates of Mandelbrot set 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.Drawing the Mandelbrot set by the method of escape lines. M. Romera et al. This method is fast and allows looking for high iterations. boundary tracehttp://www.metabit.org/~rfigura/figura-fractal/math.html boundary trace by Robert Figura drawing curves , see explanation and source code. This method is very complicated for iterations > 5. Decomposition of target set for Mandelbrot set drawing 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 thumb|binary decomposition: image with full code in C Here target set on dynamic plane is divided into 2 parts (binary decomposition = 2-decomposition ): upper half ( white) lower half (black) Division of target set induces decomposition of level sets into parts ( cells, subsets): which is colored white, which is colored black. "The Level Sets and Field Lines are superimposed, creating a sort of grid, and the "squares" of the grid are filled with N-digit binary numbers giving the first N binary digits of the external angles of field lines passing through the square. (Alternately, only the Nth binary digit is used.) Each level set is divided into 2n squares. It is easy to "read" the external arguments of points in the boundary of the Mandelbrot Set using a binary decomposition." Robert P. Munafo For binary decomposition use exp(pi) as escape radius, so that the boxes appear square (a tip from mrob). External rays of angles (measured in turns): 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.http://web.archive.org/20010415125044/www.geocities.com/CapeCanaveral/Launchpad/5113/fr27.htm| An open letter to Dr. Meech from Joyce Haslam in FRACTAL REPORT 27 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 than 2 parts ( for example 8).mandelbrot set n-th-decomposition Real Escape Time Other names of this method/algorithm are : the fully-renormalized fractional iteration count ( by Linas Vepstas in 1997)linas.org : Renormalizing the Mandelbrot Escape smooth iteration count for generalized Mandelbrot sets ( by Inigo Quilez in 2016)I Quilez : mset_smooth 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 : Other methods and speedups Colouring formula in Ultrafractal :fractalforums : What range/precision for fractional escape counts for Mandelbrot/Julia sets? smooth iter = iter + 1 + ( log(log(bailout)-log(log(cabs(z))) )/log(2) where : log(log(bailout) can be precalculated theory Description by Claude : First description : If R is large, the first z to escape satisfies (approximately)fractalforums : gradient-pallet-with-two-colors so taking logs so taking logs again so by algebra when at escape is bigger, the smooth iteration count should be smaller, so this value needs to be subtracted from the integer iteration count Alternatively this fraction can be used for interpolation, or used with arg(z) for exterior tiling / binary decomposition. Second descriptionfractalforums.org : can-anyone-help-me-understand-smooth-coloring pick a radius R > 2, then |Z| > R implies that |Z^2 + C| > |Z| and more generally that |Z| -> infinity, this gives R the name escape radius. proof is on math.stackexchange.com somewhere now suppose R is large, and n is the first iteration where |Z_n| > R. consider what happens when |Z_n| increases as you move the point C a bit further from the Mandelbrot set boundary. eventually |Z_n| > R^2, but then |Z_{n-1}| > R, so the iteration count should be n - 1. for smoothing, we want a value to add to n that is 0 when |Z_n| = R and -1 when |Z_n| = R^2. taking logs, get log |Z| is between log(R) and 2 log(R) taking logs again, get log log |Z| is between log log R and log log R + log 2 dividing by log 2, get log_2 log |Z| is between log_2 log R and log_2 log R + 1 subtracting log_2 log R gives (log_2 log |Z| - log_2 log R) is between 0 and 1 negating it gives a value between 0 and -1, as desired so the smooth iteration count is (replace 2 by P if you do Z^P + C) see also http://linas.org/art-gallery/escape/escape.html which makes a value independent of R, but that is not so useful for some colouring algorithms (e.g. smooth part of escape count doesn't align with angle of final iterate) 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 FinchMaking Mandelbrot Set Movies by Tony Finch 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 Linas Vepstas. Renormalizing the mandelbrot escape. 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)); java /** 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()); } Matemathica Here is code by Paul Nylander. It uses different formula : Python Python code using mpmath librarympmath Python library 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/M Variants : exterior DEM/M interior DEM/M Description Complex potential Description See also https://web.archive.org/web/20071008112609/http://rgba.scenesp.org/iq/trastero/fieldlines/ http://fraktal.republika.pl/mset_bottcher.html References