Difference between revisions of "Mandelbrot set"

Revision as of 07:20, 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
                waveformBM.SetPixel((cr + 2) * 50 + 300, (ci + 2) * 50, Color.Black)
            Else
                waveformBM.SetPixel((cr + 2) * 50 + 300, (ci + 2) * 50, Color.White)
            End If

    Next cr
Next ci