In the April 1987 Compute magazine, there was a program called Hyperscan, which rendered a picture called the Mandelbrot set. I was completely fascinated by this picture and design, because the Hyperscan program allowed me to zoom in to the Mandelbrot picture and explore the seemingly limitless detail of this image. I spend weeks with this program, and I printed out tons of different parts of the Mandelbrot set on my dot matrix printer.

I never understood what the Mandelbrot set actually was until recently, when I decided to research it. I knew that it was a very detailed picture that was created by a relatively simple formula, but the programming of Hyperscan was really complicated, so I could never figure out exactly what that supposedly simple formula was. I also knew that the Mandelbrot set uses complex numbers, which are numbers that have the square root of negative 1, or i.

**Explanation of how to calculate the Mandelbrot set:**

Almost all of the complexity of the formula for the Mandelbrot set is in the complex numbers. So, once you understand the math with complex numbers, and how you can plot complex numbers on a graph, the Mandelbrot set becomes pretty easy to understand.

A complex number is a number like .5+1.1i. The thing that makes this number complex is that part of it is an imaginary number, or i. i is the square root of negative 1, so you could write this complex number as .5+1.1*sqrt(-1).

The square root of negative 1 really doesn't exist, so that's why we just write i, instead. Even though i doesn't exist, then i times i, or i squared does exist, and it's just -1.

If you add two complex numbers, it's easy. So, if you add .5+1.1i to 1+1.2i, then you just add the real parts (.5 and 1) and add the imaginary parts (1.1 and 1.2), to get 1.5+2.3i.

Squaring two complex numbers is a little more complex, (a+bi)^2 is the same as (a+bi)*(a+bi), which multiplies out to a*a+a*b*i+b*i*a+b*b*i*i, which simplifies to a^2+2abi+b^2*i^2, and because i squared is negative 1, then it simplifies to a^2-b^2+2abi.

So, using our example number, (.5+1.1i)^2 = (.5)^2-(1.1)^2+2*2*1.1*i. This equals -.96+4.4i

The last concept with complex numbers is something called magnitude, which with complex numbers is also called absolute value. You use the Pythagorean theorem on the real and imaginary parts of the complex number to find this. The magnitude of a+bi is sqrt(a^2+b^2). For our example number of .5+1.1i, the magnitude is sqrt(.5^2+1.1^2), or 1.21

To graph complex numbers, you use the real part of the number as the x axis, or the horizontal axis, and you use the imaginary part of the number as the y axis, or the vertical axis. The number of .5+1.1i would correspond to the point (.5,1.1) on a regular x,y graph. On these Mandelbrot pictures, this point corresponds to the lower right corner of the graph, as the y axis is usually reversed, but since the set is a mirror image, this doesn't matter.

The upper left point of these Mandelbrot pictures is -2-1.1i, which using graphing coordinates is (-2,-1.1).

Each point in this complex plane is either part of the Mandelbrot set, or it isn't. Or another way of saying this is that any complex number that you can think of is either a part of the Mandelbrot set, or isn't. The black part in the pictures are numbers that are part of the Mandelbrot set, and the parts that aren't black aren't part of the set.

Now that we have the background in complex numbers, the Mandelbrot formula is relatively easy. It takes a complex number, which we call Z

_{0 }and we iterate it using this formula:

Z

_{n+1}-->Z

_{n}

^{2}+Z

_{0}

If we iterate a complex number a bunch of times, and it turns into a super big number, then it's not part of the Mandelbrot set. If we iterate it a bunch of times, and it never gets very big, then it is part of the Mandelbrot set.

Technically, we're supposed to iterate it an infinite number of times, and if it diverges to infinity, then it's not part of the Mandelbrot set. But there are two shortcuts, so we don't actually have to do it an infinite number of times.

- Shortcut #1: Once you iterate this more than about 200 times, then you stop getting very much more information. So, iterating this 200 times is about as good as iterating it an infinite number of times.
- Shortcut #2: There's a theorem that figured out that if the magnitude of Z ever gets bigger than 2, while iterating, then it will always eventually get to be infinity.

So a simple way of saying the Mandelbrot equation is that you iterate Z

_{n+1}-->Z_{n}^{2}+Z_{0 }(you square your last result, and add in the original number again) for 200 times, or until Z's magnitude gets to be greater than two. If Z's magnitude gets to be greater than 2, then you plot a light colored dot. If Z's magnitude never got to be greater than 2, then you plot a black colored dot.
If computers had complex data types, then the computer code for the Mandelbrot set would be pretty simple.

**Simple pseudo-code for the Mandelbrot set calculation:**

Function isMandelbrot(Z0 as Complex) As Boolean

Dim Z as Complex, Iterations as Integer

Iterations = 255

isMandelbrot = True

Do While x < Iterations

Z = Z^2 + Z0

If ABS(Z) > 2 Then

isMandelbrot = False

Exit Function

End If

x = x + 1

Loop

End Function

But, the problem with this psuedo-code is that most programming languages don't have a data type of Complex. So, we have to do the calculations that I wrote above, in a spelled out way. The good news is that this only adds a few lines of code, so it's still a pretty simple formula.

While recently researching the Mandelbrot set, I became fascinated with the lack of simple instructions that there are out there, on how to calculate this relatively simple formula. I also got really nostalgic for my old Commodore 64, so I decided to combine both. I found a really cool Commodore 64 emulator, VICE, and I relearned Commodore 64 BASIC in the quest to write some simple code that would explain the Mandelbrot set.

I present below some of the simplest Mandelbrot programs ever. The low resolution program is very low res, and produces an image that is only 40x25. The hires, or high resolution image is only 320x200. By today's standards, it would never be considered high resolution, but this does use the maximum pixels on the Commodore 64's display. I made both programs as similar as possible, so you can easily dissect and change them if you'd like.

**Programming notes:**

- {Clear} should be entered as Shift+Home, and will display as a reverse heart.
- ^ does not exist on the Commodore 64 keyboard, and should be the up arrow key (not the movement key). It's located in between * and RESTORE. When entered correctly, it will display as an up arrow.
- The two letters before the line numbers are Compute! magazines checksum, and are not entered. You can find instructions on how to display this checksum in any back issue of Compute! or in any of their fine books.
- Instead of writing (RZ*RZ + IZ*IZ)^(0.5)>2, for the magnitude calculation, I wrote (RZ*RZ + IZ*IZ)>4, which is the same thing.

**Mandelbrot function in Commodore 64 BASIC (low-resolution):**

RB 90 REM **TO CHANGE THE SHAPE, CHANGE**

JR 91 REM **XL,XU,YL,YU. REPS CHANGES **

RA 92 REM **THE RESOLUTION**

DD 100 XL = -2.000:XU = 0.500

HR 110 YL = -1.100:YU = 1.100

JX 115 REPS = 20

HC 120 WIDTH = 40:HEIGHT = 25

KD 130 XINC = (XU-XL)/WIDTH

RD 140 YINC = (YU-YL)/HEIGHT

FR 200 REM MAIN ROUTINE

FM 205 PRINT "{CLEAR}"

PJ 210 FOR J = 0 TO HEIGHT - 1

EA 220 : FOR I = 0 TO WIDTH - 1

DR 230 : GOSUB 300

RR 231 : GOSUB 400

CP 240 : NEXT I

JG 250 NEXT J

SE 290 GET A$:IF A$ = "" THEN 290

MM 299 END

FM 300 REM CALCULATE MANDELBROT

BC 310 ISMND = -1

PF 312 NREAL = XL + I * XINC

RK 313 NIMG = YL + J * YINC

XS 315 RZ = 0:IZ = 0

RS 316 R2Z = 0:I2Z = 0

EQ 320 FOR K = 1 TO REPS

CC 330 : R2Z = RZ*RZ - IZ*IZ

JQ 340 : I2Z = 2*RZ*IZ

QJ 350 : RZ = R2Z+NREAL

HX 360 : IZ = I2Z +NIMG

EE 370 : IF (RZ*RZ + IZ*IZ)>4 THEN ISMND=0:K=REPS

QA 390 NEXT K

GB 399 RETURN

GB 400 REM PLOT MANDELBROT

GE 410 COUNT = I+J*WIDTH

PJ 420 POKE 1024+COUNT,160

GR 430 POKE 55296+COUNT,2 AND NOT ISMND

BH 499 RETURN

**Mandelbrot function in Commodore 64 BASIC (high-resolution):**

RB 90 REM **TO CHANGE THE SHAPE, CHANGE**

JR 91 REM **XL,XU,YL,YU. REPS CHANGES **

RA 92 REM **THE RESOLUTION**

DD 100 XL = -2.000:XU = 0.500

HR 110 YL = -1.100:YU = 1.100

JX 115 REPS = 20

CM 120 WIDTH = 320:HEIGHT = 200

KD 130 XINC = (XU-XL)/WIDTH

RD 140 YINC = (YU-YL)/HEIGHT

FR 200 REM MAIN ROUTINE

FM 205 PRINT "{CLEAR}"

XR 207 GOSUB 500

PJ 210 FOR J = 0 TO HEIGHT - 1

EA 220 : FOR I = 0 TO WIDTH - 1

DR 230 : GOSUB 300

BX 231 : GOSUB 600

CP 240 : NEXT I

JG 250 NEXT J

SE 290 GET A$:IF A$ = "" THEN 290

EJ 299 GOTO 700

FM 300 REM CALCULATE MANDELBROT

BC 310 ISMND = -1

PF 312 NREAL = XL + I * XINC

RK 313 NIMG = YL + J * YINC

XS 315 RZ = 0:IZ = 0

RS 316 R2Z = 0:I2Z = 0

EQ 320 FOR K = 1 TO REPS

CC 330 : R2Z = RZ*RZ - IZ*IZ

JQ 340 : I2Z = 2*RZ*IZ

QJ 350 : RZ = R2Z+NREAL

HX 360 : IZ = I2Z +NIMG

EE 370 : IF (RZ*RZ + IZ*IZ)>4 THEN ISMND=0:K=REPS

QA 390 NEXT K

GB 399 RETURN

HM 500 REM CLEAR SCREEN

BB 530 POKE53272,29:POKE53265,59

CC 540 FOR R=8192 TO 16383:POKE R,0:NEXT

QS 599 RETURN

QH 600 REM DRAW HIRES

CA 610 P=8192+INT(J/8)*320+INT(I/8)*8+(JAND7)

JE 620 IF ISMND=0 THEN POKEP,PEEK(P) OR INT(2^(7-((I/8-INT(I/8))*8)))

JF 699 RETURN

CM 700 REM RETURN TO NORMAL

QH 740 POKE 53272,21:POKE 53265,27

**Mandelbrot function in Excel (VBA or Visual Basic, too):**

Function isMandelbrot(NReal As Double, NImg As Double) As Boolean

Dim ZReal As Double, ZImg As Double

Dim ZReal2 As Double, ZImg2 As Double

Dim Iterations As Integer, x As Integer

Iterations = 255

isMandelbrot = True

Do While x < Iterations

ZReal2 = (ZReal) ^ 2 - (ZImg) ^ 2

ZImg2 = 2 * ZReal * ZImg

ZReal = ZReal2 + NReal

ZImg = ZImg2 + NImg

If ((ZReal) ^ 2 + (ZImg) ^ 2) > 4 Then

isMandelbrot = False

Exit Function

End If

x = x + 1

Loop

End Function

**About the pretty colors you see on other Mandelbrot images:**

To keep this simple, I used only two colors, black, and non-black. The most common way of doing colors is to record how many iterations before a number fails the Mandelbrot test, and map various colors to the numbers of iterations before failing the test. A common number of iterations is 255, because that usually maps well to color formulas.

Links:

I added multicolor hi-res support, in 4-color glory! Also, some machine code to erase the hi-res screen, which saves a lot of time before anything interesting happens.

ReplyDelete100 xl =-2.000:xu =0.500

110 yl =-1.100:yu=1.100

115 reps =20

120 width =160:height =200

130 xinc =(xu -xl)/width

140 yinc =(yu-yl)/height

200 rem main routine

205 print "{clr}"

207 gosub 500

210 for j = 0 to height - 1

220 : for i = 0 to width -1

230 gosub 300

231 gosub 600

240 next i

250 next j

290 get a$: if a$ = ""then 290

299 goto 700

300 rem calculate mandelbrot

310 ismnd =-1

312 nreal = xl + i * xinc

313 nimg = yl + j * yinc

315 rz = 0:iz = 0

316 r2z = 0:i2z = 0

320 for k = 1 to reps

330 r2z=rz*rz-iz*iz

340 i2z=2*rz*iz

350 rz=r2z+nreal

360 iz=i2z+nimg

370 if(rz*rz+iz*iz)>4 then ismnd=0:goto395

390 next

395 if k<8 then c=3:return

396 if k<12then c=2:return

397 if k<16then c=1:return

398 c=0

400 return

500 rem clear screen

530 poke 53272, 29:poke 53265, 59

531 poke53270,peek(53270)or16

532 poke 53280, 0

541 gosub 800

549 rem set screen colors to grays

550 for r=1024 to 2023:poke r,207:next

560 for r=55296to56295:poke r,1:next

599 return

600 rem draw hires

610 p =8192 + int(j/8)*320 +int(i/4)*8+(j and7)

625 lb=2^(7-(1+2*(iand3)))

626 hb=2^(7-(2*(iand3)))

630 if c=1 then pokep,peek(p)or lb

640 if c=3 then pokep,peek(p)or lb

650 if c=2 then pokep,peek(p)or hb

699 return

700 rem return to normal

705 poke 53280, 14

740 poke 53272,21:poke 53265,27

750 poke 53270,peek(53270)and239

760 end

799 rem some machine code to erase the hires screen

800 for i= 0to 20:read a:poke 49152+i, a:next

810 data 169,0,162,32,160,0,132,33,134,34,145,33,200,208,251,232,224,64,208,244

815 data 96

820 sys49152

830 return

I improved my color code, and creating a post: http://grahamsneurons.blogspot.com/2013/11/commodore-64-mandelbrot-in-color.html

ReplyDeleteThis is so amazingly cool! Thank you so much for adding to my code, and making this so awesome!!!

ReplyDeleteNice little program. I added as a test case to my BASIC V2-compiler-for-Java-project (which can be found here: https://github.com/EgonOlsen71/basicv2 )

ReplyDeleteAlbeit that thing isn't an emulator, i.e. it doesn't actually display graphics and such, I added the option to save the rendered image to disk.

It runs approx. 60.000 times faster than an actual C64 on a 4Ghz machine and still approx. 4000 times faster than VICE in warp mode on the same machine.

Wow! EgonOlsen, your project is really cool!

Delete