## Thursday, April 26, 2012

### How to simply create the Mandelbrot set on the Commodore 64

When I was a kid, I loved my Commodore 64 computer. Back then, I had a subscription to Compute! magazine, which listed out programs that I typed into my Commodore 64. I typed in dozens of these programs, and it usually took four evenings after school for me to type one program into the computer, and another day to debug my mistakes. The games and programs really weren't really much fun, so I only ended up playing the games for about a few hours. I spent way more time typing programs into the computer than I ever spent playing the games on the computer.

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 Zand we iterate it using this formula:

Zn+1-->Zn2+Z0

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 Zn+1-->Zn2+Z0 (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
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.

1. 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.

100 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
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

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

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

4. Nice 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 )
Albeit 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.

1. Wow! EgonOlsen, your project is really cool!