Difference between revisions of "Mandelbrot set"

Revision as of 10:15, 1 May 2010

The Mandelbrot set

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.SetPixel((cr + 2) * 50 + 300, (ci + 2) * 50, Color.Black)
            Else
                complexPlane.SetPixel((cr + 2) * 50 + 300, (ci + 2) * 50, Color.White)
            End If

    Next cr
Next ci