Mandelbrot set

This page uses content from Wikipedia. The original article was at Mandelbrot_set. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

Mandelbrot set
You are encouraged to solve this task according to the task description, using any language you may know.

Generate and draw the Mandelbrot set.

BASIC

This is almost exactly the same as the pseudocode from the Wikipedia entry's "For programmers" section (which it's closely based on, of course). The image generated is very blocky ("low-res") due to the selected video mode, but it's fairly accurate. <lang qbasic> SCREEN 13 WINDOW (-2, 1.5)-(2, -1.5) FOR x0 = -2 TO 2 STEP .01

   FOR y0 = -1.5 TO 1.5 STEP .01
       x = 0
       y = 0

       iteration = 0
       maxIteration = 223

       WHILE (x * x + y * y <= (2 * 2) AND iteration < maxIteration)
           xtemp = x * x - y * y + x0
           y = 2 * x * y + y0

           x = xtemp

           iteration = iteration + 1
       WEND

       IF iteration <> maxIteration THEN
           c = iteration
       ELSE
           c = 0
       END IF

       PSET (x0, y0), c + 32
   NEXT

NEXT </lang>

C

This code uses functions from category Raster graphics operations.

<lang c>#include <stdio.h>

1. include "imglib.h"

1. define MAX_ITERATIONS 500

void draw_mandel(image img){

 unsigned int x, y;
 unsigned int width, height;
 color_component colour;
 double x0, y0, xtemp, cr, ci;
 unsigned int i;

 width = img->width;
 height = img->height;

 for (x = 0; x < width; x++) {
   for (y = 0; y < height; y++) {
     x0 = y0 = 0;
     cr = ((double)x / (double)width - 0.5);
     ci = ((double)y / (double)height - 0.5);
     i = 0;

     while (((x0*x0 + y0*y0) <= 4.0) && (i < MAX_ITERATIONS)){
       xtemp = x0*x0 - y0*y0 + cr;
       y0 = 2*x0*y0 + ci;
       x0 = xtemp;
       i++;
     }

     /* points in the set are coloured as white, others are coloured black. */
     colour = (i == MAX_ITERATIONS ? 255 : 0);

     put_pixel_clip(img, x, y, colour, colour, colour);
   }
 }

}

1. define W 640
2. define H 480

int main() {

 image img;
 FILE *o;

 img = alloc_img(W, H);
 draw_mandel(img);

 o = fopen("mandel.ppm", "w");
 output_ppm(o, img);
 fclose(o);
 free_img(img);

}</lang>

C++

This generic function assumes that the image can be accessed like a two-dimensional array of colors. It may be passed a true array (in which case the Mandelbrot set will simply be drawn into that array, which then might be saved as image file), or a class which maps the subscript operator to the pixel drawing routine of some graphics library. In the latter case, there must be functions get_first_dimension and get_second_dimension defined for that type, to be found by argument dependent lookup. The code provides those functions for built-in arrays. <lang cpp>

1. include <cstdlib>
2. include <complex>

// get dimensions for arrays template<typename ElementType, std::size_t dim1, std::size_t dim2>

std::size_t get_first_dimension(ElementType (&a)[dim1][dim2])

{

 return dim1;

}

template<typename ElementType, std::size_t dim1, std::size_t dim2>

std::size_t get_second_dimension(ElementType (&a)[dim1][dim2])

{

 return dim2;

}

template<typename ColorType, typename ImageType>

void draw_Mandelbrot(ImageType& image,                                   //where to draw the image
                     ColorType set_color, ColorType non_set_color,       //which colors to use for set/non-set points
                     double cxmin, double cxmax, double cymin, double cymax,//the rect to draw in the complex plane
                     unsigned int max_iterations)                          //the maximum number of iterations

{

 std::size_t const ixsize = get_first_dimension(ImageType);
 std::size_t const iysize = get_first_dimension(ImageType);
 for (std::size_t ix = 0; ix < ixsize; ++ix)
   for (std::size_t iy = 0; iy < iysize; ++iy)
   {
     std::complex c(cxmin + ix/(ixsize-1.0)*(cxmax-cxmin), cymin + iy/(iysize-1.0)*(cymax-cymin));
     std::complex z = 0;
     unsigned int iterations;

     for (iterations = 0; iterations < max_iterations && std::abs(z) < 2.0; ++iterations) 
     {
       z = z*z + c;
     }

     if (iterations == max_iterations)
       image[ix][iy] = set_color;
     else
       image[ix][iy] = non_set_color;
   }

} </lang>

C#

<lang csharp>using System; using System.Drawing; using System.Drawing.Imaging; using System.Threading; using System.Windows.Forms;

/// <summary> /// Generates bitmap of Mandelbrot Set and display it on the form. /// </summary> public class MandelbrotSetForm : Form {

   const double MaxValueExtent = 2.0;
   Thread thread;

   static double CalcMandelbrotSetColor(ComplexNumber c)
   {
       // from http://en.wikipedia.org/w/index.php?title=Mandelbrot_set
       const int MaxIterations = 1000;
       const double MaxNorm = MaxValueExtent * MaxValueExtent;

       int iteration = 0;
       ComplexNumber z = new ComplexNumber();
       do
       {
           z = z * z + c;
           iteration++;
       } while (z.Norm() < MaxNorm && iteration < MaxIterations);
       if (iteration < MaxIterations)
           return (double)iteration / MaxIterations;
       else
           return 0; // black
   }

   static void GenerateBitmap(Bitmap bitmap)
   {
       double scale = 2 * MaxValueExtent / Math.Min(bitmap.Width, bitmap.Height);
       for (int i = 0; i < bitmap.Height; i++)
       {
           double y = (bitmap.Height / 2 - i) * scale;
           for (int j = 0; j < bitmap.Width; j++)
           {
               double x = (j - bitmap.Width / 2) * scale;
               double color = CalcMandelbrotSetColor(new ComplexNumber(x, y));
               bitmap.SetPixel(j, i, GetColor(color));
           }
       }
   }

   static Color GetColor(double value)
   {
       const double MaxColor = 256;
       const double ContrastValue = 0.2;
       return Color.FromArgb(0, 0,
           (int)(MaxColor * Math.Pow(value, ContrastValue)));
   }
   
   public MandelbrotSetForm()
   {
       // form creation
       this.Text = "Mandelbrot Set Drawing";
       this.BackColor = System.Drawing.Color.Black;
       this.BackgroundImageLayout = System.Windows.Forms.ImageLayout.Stretch;
       this.MaximizeBox = false;
       this.StartPosition = FormStartPosition.CenterScreen;
       this.FormBorderStyle = FormBorderStyle.FixedDialog;
       this.ClientSize = new Size(640, 640);
       this.Load += new System.EventHandler(this.MainForm_Load);
   }

   void MainForm_Load(object sender, EventArgs e)
   {
       thread = new Thread(thread_Proc);
       thread.IsBackground = true;
       thread.Start(this.ClientSize);
   }

   void thread_Proc(object args)
   {
       // start from small image to provide instant display for user
       Size size = (Size)args;
       int width = 16;
       while (width * 2 < size.Width)
       {
           int height = width * size.Height / size.Width;
           Bitmap bitmap = new Bitmap(width, height, PixelFormat.Format24bppRgb);
           GenerateBitmap(bitmap);
           this.BeginInvoke(new SetNewBitmapDelegate(SetNewBitmap), bitmap);
           width *= 2;
           Thread.Sleep(200);
       }
       // then generate final image
       Bitmap finalBitmap = new Bitmap(size.Width, size.Height, PixelFormat.Format24bppRgb);
       GenerateBitmap(finalBitmap);
       this.BeginInvoke(new SetNewBitmapDelegate(SetNewBitmap), finalBitmap);
   }

   void SetNewBitmap(Bitmap image)
   {
       if (this.BackgroundImage != null)
           this.BackgroundImage.Dispose();
       this.BackgroundImage = image;
   }

   delegate void SetNewBitmapDelegate(Bitmap image);

   static void Main()
   {
       Application.Run(new MandelbrotSetForm());
   }

}

struct ComplexNumber {

   public double Re;
   public double Im;

   public ComplexNumber(double re, double im)
   {
       this.Re = re;
       this.Im = im;
   }

   public static ComplexNumber operator +(ComplexNumber x, ComplexNumber y)
   {
       return new ComplexNumber(x.Re + y.Re, x.Im + y.Im);
   }

   public static ComplexNumber operator *(ComplexNumber x, ComplexNumber y)
   {
       return new ComplexNumber(x.Re * y.Re - x.Im * y.Im,
           x.Re * y.Im + x.Im * y.Re);
   }

   public double Norm()
   {
       return Re * Re + Im * Im;
   }

} </lang>

D

Library: qd

using

Library: SDL

Library: Phobos

<lang D> module mandelbrot;

import qd;

double lensqr(cdouble c) { return c.re * c.re + c.im * c.im; }

const Limit = 150;

void main() {

 screen(640, 480);
 for (int y = 0; y < screen.h; ++y) {
   flip; events;
   for (int x = 0; x < screen.w; ++x) {
     auto
       c_x = x * 1.0 / screen.w - 0.5,
       c_y = y * 1.0 / screen.h - 0.5,
       c = c_y * 2.0i + c_x * 3.0 - 1.0,
       z = 0.0i + 0.0,
       i = 0;
     for (; i < Limit; ++i) {
       z = z * z + c;
       if (lensqr(z) > 4) break;
     }
     auto value = cast(ubyte) (i * 255.0 / Limit);
     pset(x, y, rgb(value, value, value));
   }
 }
 while (true) { flip; events; }

} </lang>

Forth

This uses grayscale image utilities. <lang Forth>500 value max-iter

mandel ( gmp F: imin imax rmin rmax -- )

 0e 0e { F: imin F: imax F: rmin F: rmax F: Zr F: Zi }
 dup bheight 0 do
   i s>f dup bheight s>f f/ imax imin f- f* imin f+ TO Zi
   dup bwidth 0 do
     i s>f dup bwidth s>f f/ rmax rmin f- f* rmin f+ TO Zr
     Zr Zi max-iter
     begin  1- dup
     while  fover fdup f* fover fdup f*
            fover fover f+ 4e f<
     while  f- Zr f+
            frot frot f* 2e f* Zi f+
     repeat fdrop fdrop
            drop 0        \ for a pretty grayscale image, replace with: 255 max-iter */
     else   drop 255
     then   fdrop fdrop
     over i j rot g!
   loop
 loop    drop ;

80 24 graymap dup -1e 1e -2e 1e mandel dup gshow free bye</lang>

Haskell

Translation of: Ruby

<lang haskell>import Data.Complex import Control.Monad (forM_)

mandelbrot a = iterate (\z -> z^2 + a) a !! 50

main = mapM_ putStrLn [[if magnitude (mandelbrot (x :+ y)) < 2 then '*' else ' '

                          | x <- [-2, -1.9685 .. 0.5]]
                      | y <- [1, 0.95 .. -1]]</lang>

haXe

This version compiles for flash version 9 or grater. The compilation command is <lang>haxe -swf9 mandelbrot.swf -main Mandelbrot</lang>

<lang haxe> class Mandelbrot extends flash.display.Sprite {

   inline static var MAX_ITER = 255;

   public static function main() {
       var w = flash.Lib.current.stage.stageWidth;
       var h = flash.Lib.current.stage.stageHeight;
       var mandelbrot = new Mandelbrot(w, h);
       flash.Lib.current.stage.addChild(mandelbrot);
       mandelbrot.drawMandelbrot();
   }

   var image:flash.display.BitmapData;
   public function new(width, height) {
       super();
       var bitmap:flash.display.Bitmap;
       image = new flash.display.BitmapData(width, height, false);
       bitmap = new flash.display.Bitmap(image);
       this.addChild(bitmap);
   }

   public function drawMandelbrot() {
       image.lock();
       var step_x = 3.0 / (image.width-1);
       var step_y = 2.0 / (image.height-1);
       for (i in 0...image.height) {
           var ci = i * step_y - 1.0;
           for (j in 0...image.width) {
               var k = 0;
               var zr = 0.0;
               var zi = 0.0;
               var cr = j * step_x - 2.0;
               while (k <= MAX_ITER && (zr*zr + zi*zi) <= 4) {
                   var temp = zr*zr - zi*zi + cr;
                   zi = 2*zr*zi + ci;
                   zr = temp;
                   k ++;
               }
               paint(j, i, k);
           }
       }
       image.unlock();
   }

   inline function paint(x, y, iter) {
       var color = iter > MAX_ITER? 0 : iter * 0x100;
       image.setPixel(x, y, color);
   }

} </lang>

J

The characteristic function of the Mandelbrot can be defined as follows: <lang j> mcf=. (<: 2:)@|@(] ((*:@] + [)^:((<: 2:)@|@])^:1000) 0:) NB. 1000 iterations test </lang> The Mandelbrot set can be drawn as follows: <lang j> domain=. |.@|:@({.@[ + ] *~ j./&i.&>/@+.@(1j1 + ] %~ -~/@[))&>/

load'graph' viewmat mcf "0 @ domain (_2j_1 1j1) ; 0.01 NB. Complex interval and resolution </lang>

Modula-3

<lang modula3>MODULE Mandelbrot EXPORTS Main;

IMPORT Wr, Stdio, Fmt, Word;

CONST m = 50;

     limit2 = 4.0;

TYPE UByte = BITS 8 FOR [0..16_FF];

VAR width := 200;

   height := 200;
   bitnum: CARDINAL := 0;
   byteacc: UByte := 0;
   isOverLimit: BOOLEAN;
   Zr, Zi, Cr, Ci, Tr, Ti: REAL;

BEGIN

 Wr.PutText(Stdio.stdout, "P4\n" & Fmt.Int(width) & " " & Fmt.Int(height) & "\n");

 FOR y := 0 TO height - 1 DO
   FOR x := 0 TO width - 1 DO
     Zr := 0.0; Zi := 0.0;
     Cr := 2.0 * FLOAT(x) / FLOAT(width) - 1.5;
     Ci := 2.0 * FLOAT(y) / FLOAT(height) - 1.0;
     
     FOR i := 1 TO m + 1 DO
       Tr := Zr*Zr - Zi*Zi + Cr;
       Ti := 2.0*Zr*Zi + Ci;
       Zr := Tr; Zi := Ti;
       isOverLimit := Zr*Zr + Zi*Zi > limit2;
       IF isOverLimit THEN EXIT; END;
     END;
     
     IF isOverLimit THEN
       byteacc := Word.Xor(Word.LeftShift(byteacc, 1), 16_00);
     ELSE
       byteacc := Word.Xor(Word.LeftShift(byteacc, 1), 16_01);
     END;

     INC(bitnum);
     
     IF bitnum = 8 THEN
       Wr.PutChar(Stdio.stdout, VAL(byteacc, CHAR));
       byteacc := 0;
       bitnum := 0;
     ELSIF x = width - 1 THEN
       byteacc := Word.LeftShift(byteacc, 8 - (width MOD 8));
       Wr.PutChar(Stdio.stdout, VAL(byteacc, CHAR));
       byteacc := 0;
       bitnum := 0
     END;
     Wr.Flush(Stdio.stdout);
   END;
 END;

END Mandelbrot.</lang>

OCaml

<lang OCaml>#load "graphics.cma"

open Graphics;;

let mandelbrot xMin xMax yMin yMax xPixels yPixels maxIter =

 let rec mandelbrotIterator z c n =
   if (Complex.norm z) > 2.0 then false else
     match n with
     | 0 -> true
     | n -> let z' = Complex.add (Complex.mul z z) c in
            mandelbrotIterator z' c (n-1) in
 Graphics.open_graph
   (" "^(string_of_int xPixels)^"x"^(string_of_int yPixels)^"+0-0");
 let dx = (xMax -. xMin) /. (float_of_int xPixels) 
 and dy = (yMax -. yMin) /. (float_of_int yPixels) in
 for xi = 0 to xPixels - 1 do
   for yi = 0 to yPixels - 1 do
     let c = {Complex.re = xMin +. (dx *. float_of_int xi);
              Complex.im = yMin +. (dy *. float_of_int yi)} in
     if (mandelbrotIterator Complex.zero c maxIter) then
       (Graphics.set_color Graphics.white;
        Graphics.plot xi yi )
     else
       (Graphics.set_color Graphics.black;
        Graphics.plot xi yi )
   done
 done;;

mandelbrot (-1.5) 0.5 (-1.0) 1.0 500 500 200; read_line();; </lang>

Octave

This code runs rather slowly and produces coloured Mandelbrot set by accident (output image here).

<lang octave>#! /usr/bin/octave -qf global width = 200; global height = 200; maxiter = 100;

z0 = 0; global cmax = 1 + i; global cmin = -2 - i;

function cs = pscale(c)

 global cmax;
 global cmin;
 global width;
 global height;
 persistent px = (real(cmax-cmin))/width;
 persistent py = (imag(cmax-cmin))/height;
 cs = real(cmin) + px*real(c) + i*(imag(cmin) + py*imag(c));

endfunction

ms = zeros(width, height); for x = 0:width-1

 for y = 0:height-1
   z0 = 0;
   c = pscale(x+y*i);
   for ic = 1:maxiter
     z1 = z0^2 + c;
     if ( abs(z1) > 2 ) break; endif
     z0 = z1;
   endfor
   ms(x+1, y+1) = ic/maxiter;
 endfor

endfor

saveimage("mandel.ppm", round(ms .* 255).', "ppm");</lang>

Perl

translation / optimization of the ruby solution <lang perl> use Math::Complex;

sub mandelbrot {

   my ($z, $c) = @_[0,0];
   for (1 .. 20) {
       $z = $z * $z + $c;
       return $_ if abs $z > 2;
   }

}

for (my $y = 1; $y >= -1; $y -= 0.05) {

   for (my $x = -2; $x <= 0.5; $x += 0.0315)
       {print $_ ? ' ' : '#' for mandelbrot $x + $y * i}
   print "\n"

} </lang>

PostScript

<lang postscript>%!PS-Adobe-2.0 %%BoundingBox: 0 0 300 200 %%EndComments /origstate save def /ld {load def} bind def /m /moveto ld /g /setgray ld /dot { currentpoint 1 0 360 arc fill } bind def %%EndProlog % param /maxiter 200 def % complex manipulation /complex { 2 array astore } def /real { 0 get } def /imag { 1 get } def /cmul { /a exch def /b exch def

   a real b real mul
   a imag b imag mul sub
   a real b imag mul
   a imag b real mul add
   2 array astore

} def /cadd { aload pop 3 -1 roll aload pop

   3 -1 roll add
   3 1 roll add exch 2 array astore

} def /cconj { aload pop neg 2 array astore } def /cabs2 { dup cconj cmul 0 get} def % mandel 200 100 translate -200 1 100 { /x exch def

 -100 1 100 { /y exch def
   /z0 0.0 0.0 complex def
   0 1 maxiter { /iter exch def

x 100 div y 100 div complex z0 z0 cmul cadd dup /z0 exch def cabs2 4 gt {exit} if

   } for
   iter maxiter div g
   x y m dot
 } for

} for % showpage origstate restore %%EOF</lang>

R

<lang R>createmandel <- function(zfrom, zto, inc=0.01, maxiter=100) {

 x <- c()
 y <- c()
 cost <- zfrom
 while( Re(cost) <= Re(zto) ) {
   cost <- Re(cost) + Im(zfrom)*1i
   while ( Im(cost) <= Im(zto) ) {
     j <- 0
     z <- 0
     while( (abs(z) < 2) && (j < maxiter)) {
       z <- z^2 + cost
       j <- j + 1
     }
     if ( j == maxiter ) {
       x <- c(x, Re(cost))
       y <- c(y, Im(cost))
     }
     cost <- cost + inc*1i
   }
   cost <- cost + inc
 }
 plot(x,y, pch=".")

}

createmandel(-2-1i, 1+1i)</lang>

Ruby

Text only, prints an 80-char by 41-line depiction. Found here. <lang ruby>require 'complex'

def mandelbrot(a)

 Array.new(50,a).inject(a) { |z,c| z*z + c }

end

(1.0).step(-1,-0.05) do |y|

 (-2.0).step(0.5,0.0315) do |x|
   print mandelbrot(Complex(x,y)).abs < 2 ? '*' : ' '
 end
 puts

end</lang>

Translation of: Tcl

Uses the text progress bar from Median filter#Ruby <lang ruby>class RGBColour

 def self.mandel_colour(i)
   self.new( 16*(i % 15), 32*(i % 7), 8*(i % 31) )
 end

end

class Pixmap

 def self.mandelbrot(width, height)
   mandel = Pixmap.new(width,height)
   pb = ProgressBar.new(width) if $DEBUG
   width.times do |x|
     height.times do |y|
       x_ish = Float(x - width*11/15) / (width/3)
       y_ish = Float(y - height/2) / (height*3/10)
       mandel[x,y] = RGBColour.mandel_colour(mandel_iters(x_ish, y_ish))
     end
     pb.update(x) if $DEBUG
   end
   pb.close if $DEBUG
   mandel
 end

 def self.mandel_iters(cx,cy)
   x = y = 0.0
   count = 0
   while Math.hypot(x,y) < 2 and count < 255
     x, y = (x**2 - y**2 + cx), (2*x*y + cy)
     count += 1
   end
   count
 end 

end

Pixmap.mandelbrot(300,300).save('mandel.ppm')</lang>

Tcl

Library: Tk

This code makes extensive use of Tk's built-in photo image system, which provides a 32-bit RGBA plotting surface that can be then quickly drawn in any number of places in the application. It uses a computational color scheme that was easy to code... <lang tcl>package require Tk

proc mandelIters {cx cy} {

   set x [set y 0.0]
   for {set count 0} {hypot($x,$y) < 2 && $count < 255} {incr count} {
       set x1 [expr {$x*$x - $y*$y + $cx}]
       set y1 [expr {2*$x*$y + $cy}]
       set x $x1; set y $y1
   }
   return $count

} proc mandelColor {iter} {

   set r [expr {16*($iter % 15)}]
   set g [expr {32*($iter % 7)}]
   set b [expr {8*($iter % 31)}]
   format "#%02x%02x%02x" $r $g $b

} image create photo mandel -width 300 -height 300

1. Build picture in strips, updating as we go so we have "progress" monitoring
2. Also set the cursor to tell the user to wait while we work.

pack [label .mandel -image mandel -cursor watch] update for {set x 0} {$x < 300} {incr x} {

   for {set y 0} {$y < 300} {incr y} {
       set i [mandelIters [expr {($x-220)/100.}] [expr {($y-150)/90.}]]
       mandel put [mandelColor $i] -to $x $y
   }
   update

} .mandel configure -cursor {}</lang>

XSLT

The fact that you can create an image of the Mandelbrot Set with XSLT is sometimes under-appreciated. However, it has been discussed extensively on the internet so is best not reproduced here, not when the code can be executed directly in your browser at that site.