Mandelbrot Set Calculation

October 29, 2019 • edited April 10, 2020

The Mandelbrot set is one of the most famous fractals and it is really simple to draw. Personally I enjoy a lot seeing how simple rules lead to complex patterns.

Mandelbrot set image

Mandelbrot set representation from wikipedia.

Definition

The Mandelbrot set is defined by a set of complex numbers for which the following function does not diverge when iterated from z = 0. This means that the value remains bounded in absolute value.

$$ f_c(z) = z^2 + c $$

The following image displays the axes to see where the set is located in space. x axis is the real part and y axis is the imaginary part.

Mandelbrot set position

Keep things real

Whit the previous image in mind, we can start playing with some sample numbers to understand how to generate the image. Let’s pick a point outside the set (1,0) and perform the iterations with the previously defined function.

$$ z_0 := 0 $$ $$ c := 1 + 0i $$ $$ z_1 = f_c(z_0) = z^2 + c $$ $$ z_1 = f_c(0) = 0^2 + 1 + 0i = 1$$

After the fist iteration, $z_1 = 1$ and for the next iteration we just need to input the value to the function again.

$$ z_1 := 1 $$ $$ c := 1 + 0i $$ $$ z_2 = f_c(z_1) = z^2 + c $$ $$ z_2 = f_c(0) = 1^2 + 1 + 0i = 2$$

if we continue the iterations, we can see that the value keeps growing forever.

$$z_0 = 0$$ $$z_1 = 1$$ $$z_2 = 2$$ $$z_3 = 5$$ $$z_4 = 26$$

This means that the point doesn’t belong to the set, as we know from the image. Let’s pick a point that we know that belongs to the set. (-1,0). if we are right, the iterations must keep within a bounded value.

$$ z_0 := 0 $$ $$ c := -1 + 0i $$ $$ z_1 = f_c(z_0) = z^2 + c $$ $$ z_1 = f_c(0) = 0^2 + -1 + 0i = -1$$ $$ z_2 = f_c(z_1) = z^2 + c $$ $$ z_2 = f_c(0) = (-1)^2 + -1 + 0i = 0$$

Wops… turns out that $z_2 = z_0$ and this means that we are stuck in a infinite loop. You can iterate as much as you want but values will be either -1 or 0 for ever.

Let your imagination fly

Things get interesting the you add imaginary numbers to the previous equation so lets pick another point, but this time with imaginary values. (-0.5, 0.5)

Remember that $i^2 = -1$

$$ z_0 := 0 $$ $$ c := -0.5 + 0.5i $$ $$ z_1 = f_c(0) = 0^2 + -0.5 + 0.5i = -0.5 + 0.5i $$ $$ z_2 = f_c(z_1) = (-0.5 + 0.5i)^2 + -0.5 + 0.5i = 0.25 - 2*0.25i - 0.25 - 0.5 + 0.5i $$ $$ z_2 = - 0.5 $$

Let’s keep iterating and see what we get. Remember that we expect this point to belong to the set, so it should stay within a bounded value.

$$z_1=(-0.500000+0.500000i)$$

$$z_2=(-0.500000-0.000000i)$$

$$z_3=(-0.250000+0.500000i)$$

$$z_4=(-0.687500+0.250000i)$$

$$z_5=(-0.089844+0.156250i)$$

$$z_6=(-0.516342+0.471924i)$$

$$z_7=(-0.456103+0.012652i)$$

$$z_8=(-0.292130+0.488459i)$$

$$z_9=(-0.653252+0.214613i)$$

$$z_{10}=(-0.119320+0.219608i)$$

$$z_{20}=(-0.170860+0.314224i)$$

$$z_{30}=(-0.238482+0.376027i)$$

$$z_{40}=(-0.305512+0.399379i)$$

$$z_{50}=(-0.360070+0.401208i)$$

$$z_{60}=(-0.403868+0.390242i)$$

After 60 iterations, looks like the value stays within a bounded value but there is no guarantee that the value will remain here and the only way to know it is to keep iterating. That is why we need to define a maximum number of iterations to give up.

Last try with an imaginary number outside the set. (0,1)

$$ z_2 = (1i)^2 + i = -1+i $$ $$ z_3 = (-1+i)^2 + i = 1 -2i -1 + i = -i $$ $$ z_4 = (-i)^2 + i = 1 + i = -1 + i $$ $$ z_5 = (-1+i)^2 + i = 1 -2i -1 + i = -i $$

What happen here? $z_2 = z_5$ and this means that we are stuck in a infinite loop therefore the point belongs to the set. This is just a special point that belongs to the set but if you move a minimal distance from it, the iteration quickly scapes to infinity. You can visualize this as balancing a ball in a pencil. It is possible but the slightest perturbation will make the ball fall. For this reason, you don’t see anything if you zoom in this zone of the Mandelbrot set.

Let the computer do its work

I’m tired of manually calculating this values so lets create some code to do it for us.

First, we need to have a formal definition of the Mandelbrot set. This can be obtained from wikipedia.

$$ z_{n+1} = z^2_n+c $$

$$ c ∈ M ↔ limsup_{z->∞}|z_{n+1}| < 2 $$

The important information here is that we already know the condition of “diverge”. When de module of the value z is bigger than 2, we know that the point is not in the set.

Mandelbrot point

I’m going to create a structure mandelbrot.Point to hold the necessary information to calculate any given point of the set. This will contain the methods Calculate, Diverges and Iterations.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

package mandelbrot

import ()

type Point struct {
  Point      complex128
  iterations int
  z          complex128
}

// NewPoint returns a new point at a given coordenates
func NewPoint(r, i float64) *Point {
  return &Point{
    Point: complex(r, i),
  }
}

// Calculate performs as many calculations as MaxIterations to determine if the point belongs to the set or not
func (m *Point) Calculate(MaxIterations int) {
  panic("To be implemented")
}

// Diverges returns whether the points diverges from the set.
func (m *Point) Diverges() bool {
  panic("To be implemented")
}

// Iterations returns the number of performed iterations.
func (m *Point) Iterations() int {
  return m.Iterations()
}

it is time to write tests to develop this functionality. TDD fits really well in cases like this were we already have a clear definition of what are all the outputs.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

package mandelbrot_test

import (
  "testing"

  . "github.com/metalblueberry/mandelbrot/mandelbrot"
)

type testMandelbrotPointCases struct {
  point      *Point
  iterations int
  diverges   bool
}

func TestMandelbrotPoint(t *testing.T) {
  tests := []testMandelbrotPointCases{
    testMandelbrotPointCases{
      point:      NewPoint(1, 0),
      iterations: 100,
      diverges:   true,
    },
    testMandelbrotPointCases{
      point:      NewPoint(-1, 0),
      iterations: 100,
      diverges:   false,
    },
    testMandelbrotPointCases{
      point:      NewPoint(-0.5, 0.5),
      iterations: 100,
      diverges:   false,
    },
    testMandelbrotPointCases{
      point:      NewPoint(0, 1),
      iterations: 100,
      diverges:   true,
    },
  }

  for i, test := range tests {
    t.Log(test)
    test.point.Calculate(test.iterations)
    if test.point.Diverges() != test.diverges {
      t.Errorf("Test %d failed, Point %f diverges %t expected %t", i, test.point.Point, test.point.Diverges(), test.diverges)
    }
  }
}

You may notice that the point (0,1) is expected to diverge but we already know that this is a special case that does not diverge. If you debug the program, you will notice that it diverges after a few iterations due to numerical errors. This is not a problem though, the representation of the set is still perfectly fine.

The code is in fact really simple with go because there is support for complex numbers that handles all the special cases with $i^2$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

func (m *Point) Calculate(MaxIterations int) {
  for !m.Diverges() && m.iterations < MaxIterations {
    m.iterations++
    m.z = cmplx.Pow(m.z, 2) + m.Point
    log.Printf("z_%d=%f\n", m.iterations, m.z)
  }
}

func (m *Point) Diverges() bool {
  return cmplx.Abs(m.z) > complex(2, 0)
}
Click here to see all code 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

package mandelbrot

import (
  "log"
  "math/cmplx"
)

type Point struct {
  Point      complex128
  iterations int
  z          complex128
}

// NewPoint returns a new point at a given coordenates
func NewPoint(r, i float64) *Point {
  return &Point{
    Point: complex(r, i),
  }
}

// Calculate performs as many calculations as MaxIterations to determine if the point belongs to the set or not
func (m *Point) Calculate(MaxIterations int) {
  for !m.Diverges() && m.iterations < MaxIterations {
    m.iterations++
    m.z = cmplx.Pow(m.z, 2) + m.Point
    log.Printf("z_%d=%f\n", m.iterations, m.z)
  }
}

// Diverges returns whether the points diverges from the set.
func (m *Point) Diverges() bool {
  return cmplx.Abs(m.z) > complex(2, 0)
}

// Iterations returns the number of performed iterations.
func (m *Point) Iterations() int {
  return m.Iterations()
}

Mandelbrot set

The mandelbrot.Set will hold all the mandelbrot.Points. To ease the usage, I’ve added a constructor of rectangular Areas centered at a given point. This will come handy when zooming in the future.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

type Area struct {
  HorizontalResolution int
  VerticalResolution   int
  TopLeft              complex128
  BottomRight          complex128
  MaxIterations        int
  Points               []Point
}

//NewAreaCentered creates a mandelbrot.Area with squared shape centered area at x,y of width = 2*area
func NewAreaCentered(Resolution, MaxIterations int, x, y, area float64) *Area {
  return &Area{
    HorizontalResolution: Resolution,
    VerticalResolution:   Resolution,
    MaxIterations:        MaxIterations,
    TopLeft:              complex(x-area, y+area),
    BottomRight:          complex(x+area, y-area),
    Points:               make([]Point, Resolution*Resolution),
  }
}

You may notice that the Points variable is a single dimension array. The idea behind this is to have a single memory allocation for all the data. This will increase the difficulty to index a given point, but we will easily solve this by two helper methods ForIndex and IndexFor that will server to obtain the valid index for a given x,y coordinates and the x,y coordinates of a given index.

1
2
3
4
5
6
7
8
9

func (a *Area) IndexFor(x, y int) int {
  return x + y*a.HorizontalResolution
}

func (a *Area) ForIndex(i int) (x, y int) {
  y = i / a.VerticalResolution
  x = i % a.HorizontalResolution
  return x, y
}

Before using the Area, we need to initialize the Points array depending on the resolution. For this purpose let’s add a method Init to the struct that for each point. To get the complex number for a x,y position, I’m creating the method getNumber. So in the Init method I’m iterating over all the pixels defined by the resolution getting a NewPoint for each one.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

func (a *Area) getNumber(x, y int) (r, i float64) {
  TopLeftReal := real(a.TopLeft)
  TopLeftImag := imag(a.TopLeft)
  BottomRightReal := real(a.BottomRight)
  BottomRightImag := imag(a.BottomRight)

  r = TopLeftReal + (float64(x)/float64(a.HorizontalResolution))*(BottomRightReal-TopLeftReal)
  i = TopLeftImag + (float64(y)/float64(a.VerticalResolution))*(BottomRightImag-TopLeftImag)
  return r, i
}

func (a *Area) Init() {
  for x := 0; x < a.HorizontalResolution; x++ {
    for y := 0; y < a.VerticalResolution; y++ {
      point := NewPoint(a.getNumber(x, y))
      a.SetPoint(x, y, point)
    }
  }
}

Now I’m going to create a Calculate method to start the computational task. As it will be a long process, I want to do it in a separate go routine and have some kind of progress give feedback to the main thread. The method returns a channel that we can read to get the point that has been just analyzed by the function and we can divide this by the total count to know the overall progress. I don’t want the function to get blocked if the user is not reading the progress channel and that is why you see the select block at the end of the for loop. Another think is that the Points variable is a slice of struct, this means that we need to reassign the value in order to update the content of the slice.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

func (a *Area) Calculate() (progress chan int) {
  progress = make(chan int)
  go func() {
    defer close(progress)
    for i, pixel := range a.Points {
      pixel.Calculate(a.MaxIterations)
      a.Points[i] = pixel

      select {
      case progress <- i:
      default:
      }
    }
  }()
  return
}
Click here to see all code 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

package mandelbrot

type Area struct {
  HorizontalResolution int
  VerticalResolution   int
  TopLeft              complex128
  BottomRight          complex128
  MaxIterations        int
  Points               []Point
}

func NewAreaCentered(Resolution, MaxIterations int, x, y, area float64) *Area {
  return &Area{
    HorizontalResolution: Resolution,
    VerticalResolution:   Resolution,
    MaxIterations:        MaxIterations,
    TopLeft:              complex(x-area, y+area),
    BottomRight:          complex(x+area, y-area),
    Points:               make([]Point, Resolution*Resolution),
  }
}

func (a *Area) Init() {
  for x := 0; x < a.HorizontalResolution; x++ {
    for y := 0; y < a.VerticalResolution; y++ {
      point := NewPoint(a.getNumber(x, y))
      a.SetPoint(x, y, point)
    }
  }
}

func (a *Area) Calculate() (progress chan int) {
  progress = make(chan int)
  go func() {
    defer close(progress)
    for i, pixel := range a.Points {
      pixel.Calculate(a.MaxIterations)
      a.Points[i] = pixel

      select {
      case progress <- i:
      default:
      }
    }
  }()
  return
}

func (a *Area) IndexFor(x, y int) int {
  return x + y*a.HorizontalResolution
}

func (a *Area) ForIndex(i int) (x, y int) {
  y = i / a.VerticalResolution
  x = i % a.HorizontalResolution
  return x, y
}

// SetPoint changes the address of the point located at x,y
func (a *Area) SetPoint(x, y int, p Point) {
  a.Points[a.IndexFor(x, y)] = p
}

// GetPoint changes the address of the point located at x,y
func (a *Area) GetPoint(x, y int) Point {
  return a.Points[a.IndexFor(x, y)]
}

// getNumber gives the real and imaginary parts for the complex number located at the given x,y coordenates in the given resolution.
func (a *Area) getNumber(x, y int) (r, i float64) {
  TopLeftReal := real(a.TopLeft)
  TopLeftImag := imag(a.TopLeft)
  BottomRightReal := real(a.BottomRight)
  BottomRightImag := imag(a.BottomRight)

  r = TopLeftReal + (float64(x)/float64(a.HorizontalResolution))*(BottomRightReal-TopLeftReal)
  i = TopLeftImag + (float64(y)/float64(a.VerticalResolution))*(BottomRightImag-TopLeftImag)
  return r, i
}

Generate Mandelbrot Data

Now it is time to generate a simple main program to calculate the data for a given section of the complex number space. I’m going to use the image package to generate a png file of the given section and write it to the standard output. I’ve chosen a squared area centered at (-0.5, 0) of size 1.5 with 5000x5000 resolution using 200 iterations at maximum.

After calling set.Calculate I’m using range progress to get the updates and to print feedback to Stderr. To avoid printing a line for every pixel, I’ve added a <-time.After(time.Second * 1) to limit the prints to 1 per second.

The process of generating the image is as simple as initializing it with the Horizontal and Vertical resolutions and then fill each pixel using the method setRGBA. I’m not using any complex algorithm for choosing the color. For simplicity I’m just using red and the intensity depends on the iterations needed to tell if the point belongs to the set or not. the more iterations, the darker the color. That’s why the points that belong to the set are black.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

package main

import (
  "fmt"
  "image"
  "image/color"
  "image/jpeg"
  "os"
  "time"

  "github.com/metalblueberry/mandelbrot/mandelbrot"
)

func main() {
  set := mandelbrot.NewAreaCentered(5000, 200, -0.5, 0, 1.5)
  set.Init()
  progress := set.Calculate()
  for p := range progress {
    fmt.Fprintf(os.Stderr, "%d\r", 100*p/len(set.Points))
    <-time.After(time.Second * 1)
  }

  img := image.NewRGBA(image.Rect(0, 0, set.HorizontalResolution, set.VerticalResolution))

  for x := 0; x < set.HorizontalResolution; x++ {
    for y := 0; y < set.VerticalResolution; y++ {
      point := set.GetPoint(x, y)
      intensity := 255 - (255 * point.Iterations() / set.MaxIterations)
      img.SetRGBA(x, y, color.RGBA{R: uint8(intensity), A: 255})
    }
  }
  err := jpeg.Encode(os.Stdout, img, &jpeg.Options{Quality: 90})
  if err != nil {
    panic(err)
  }
}

Now you can build or run the code directly and pipe the Stdout to a png file. This code generates the following image.

mandelbrot-img

And with just modifying the input parameters, we can zoom in a given section. For example, this image is generated from the following input.

1

mandelbrot.NewAreaCentered(5000, 350, -0.7463, 0.1102, 0.005)

mandelbrot-zoom

Conclusion

I’ve plans to keep improving this code. In one hand, I would like to create an API to send the image data to a web browser. Then we will be able to use JS to build a simple viewer to ease the work of requesting zooms in a given area. In the other hand, I would like to speedup the calculation by splitting the area in chunks and compute them in parallel.

You can explore the repository to see the whole code and to play with it here.

Thank you for reading and feel free to leave a comment bellow, I will be really happy to hear from you.

howtogomandelbrot

←Go CPU profiling

→Create a Blog Like This