Mandelbrot set

The Mandelbrot set

Detail of the Mandelbrot set, each pixel is coloured by the number of iterations it took to escape

The formula for the mandelbrot set using complex numbers is:

Z = Z² + C

If Z escapes to infinity then we are outside the set, otherwise we are inside.

C is a point in the complex plane where i is on the vertical axis and real numbers are on the horizontal axis. A complex number has the form of x+iy where x and y are real numbers and i is the sqare root of -1. so the formula can also be written:

(Zr+Zi) = ((Zr+Zi) * (Zr+Zi)) + (Cr+Ci)

Complex addition has the form:

(a+bi)+(c+di) = (a+c)+i(b+d)

Complex multiplication has the form:

(a+bi)*(c+di) = (ac-bd)+i(ad+bc)

so

Z² = ((Zr+Zi) * (Zr+Zi)) = (Zr*Zr - Zir*Zir) + i(Zr*Zir + Zir*Zr)

and

Z² + C = (Zr*Zr - Zir*Zir) + i(Zr*Zir + Zir*Zr) + (Cr+Ci) = (Zr*Zr - Zir*Zir) + Cr + i((Zr*Zir + Zir*Zr) + Cir)

that gives us

Zr = (Zr*Zr - Zir*Zir) + Cr

Zi = i((Zr*Zir + Zir*Zr) + Cir)

we don't need the imaginary part to determine if Z escapes to infinity so we can use

Zir = (Zr*Zir + Zir*Zr) + Cir

we must remember to use temporary values to calculate Zr so it does not affect the result of Zi.

Z = C when we enter the loop. If Z escapes to infinity at any time during the iterations then C is outside the mandelbrot set. When C is confirmed to be outside or inside we plot a corresponding white or black pixel at coordinate C. After N iterations we give up and assume that Z will never escape to infinity to make sure the program finishes in a reasonable time.

' Floating point version
For Ci = -2 To 2 Step 0.01
    For Cr = -2 To 2 Step 0.008
        Zr = Cr
        Zi = Ci

        For i = 1 To 100  'Give up after 100 iterations to trade accuracy for speed
            Zrt = (Zr * Zr - Zi * Zi) + Cr
            Zi = (Zr * Zi + Zi * Zr) + Ci
            Zr = Zrt
            If (Zr * Zr + Zi * Zi) > 4 Then Exit For
        Next i
       
        If i = 101 Then
            complexPlane.PSet (Cr, Ci), RGB(0, 0, 0)        'Looks like we are inside
        Else
            complexPlane.PSet (Cr, Ci), RGB(255, 255, 255)  'Infinity
        End If
   
    Next Cr
Next Ci

' Integer version with 24 bit fractional resolution
For ci = -2 * (2 ^ 24) To 2 * (2 ^ 24) Step 0.01 * (2 ^ 24)
    For cr = -2 * (2 ^ 24) To 2 * (2 ^ 24) Step 0.008 * (2 ^ 24)
            zr = cr
            zi = ci

            For i = 1 To 100  'Give up after 100 iterations to trade accuracy for speed
                zr2 = zr >> 12
                zr2 *= zr2

                zi2 = zi >> 12
                zi2 *= zi2

                zi = (((zr >> 12) * (zi >> 12)) * 2) + ci
                zr = (zr2 - zi2) + cr

                If (zr2 + zi2) > (4 << 24) Then Exit For
            Next i

            If i = 101 Then
                complexPlane.Pset(cr, ci), RGB(0, 0, 0)
            Else
                complexPlane.Pset(cr, ci), RGB(255, 255, 255)
            End If

    Next cr
Next ci

// F# version using complex datatype
open System.Numerics

let maxIteration = 5000

let mandelbrot c =
  let rec loop z iteration =
    if Complex.Abs(Complex.Pow(z,2.0)) >= 4.0 then iteration
    elif iteration = maxIteration then 0
    else loop ((z * z) + c) (iteration + 1)
  loop c 0

for y in -1.2..0.15..1.2 do
  for x in -2.0..0.07..0.9 do
    printf "%c" " .:-;!+>|)&IH8%#".[mandelbrot(Complex(x, y)) &&& 15]
  printfn ""