Fractals are commonly referred to as images which are obtained by performing several calculations. One property that sets those images apart from others is the fact that if they are zoomed in, some particular shapes will show off again and again on a regular basis.
A Mandelbrot set is a collection of points that have some particular properties. Those points can be found by using a math law. This law is based on several iterations. That is, given a point of the plane, we must use its coordinates to make a recursive computation, several times. If no matter how many times we do the computation, the result remains a finite number, than the current point belongs to the Mandelbrot set. Instead, if after a certain number of computations the result goes to infinity (that is, it keeps becoming considerably bigger after each computation), then the point doesn’t belong to the Mandelbrot set.
The law of the Mandelbrot fractal is as follows:
|z0 = 0 | |z_n+1 = z_n^2 + c
This law allows to compute each zn number (for n = 0, 1, 2, ….). So we start from z0 = 0 (n = 0). Then, we must evaluate zn+1 (that is, z1), which is: z1 = z0^2+c = 0 + c = c. Then, we have z2 = z1^2 + c = c^2 + c and so on…. The law is actually pretty simple. Each computation is a recursive iteration, as any previous result is used to obtain a new result.
Now, the problem is that those numbers are not conventional numbers (that is, real numbers). Instead, those are complex numbers. And a complex number is the sum of a real number and an imaginary number. So, we just have to take a look at imaginary numbers.
As a computer deals only with real numbers, if you ask your Commodore 64 to compute the square root of -1, all you will get is the ?ILLEGAL QUANTITY ERROR message. But, in mathematics such a quantity does exit, and it is called i.
So, i equals the square root of -1. As a consequence, i^2 equals -1. This is quite important to remember.
An imaginary number is any real number multiplied by i. So, numbers like:
4i, 5i, -1.5i
are all examples of imaginary numbers.
A complex number, being the sum of a real number and an imaginary number, is a number such as:
4.7 + 2.3i
So, a generic complex number is in the form:
Z = a + bi
or
C = a + bi (which is the same).
In the Mandelbrot calculations, complex numbers to test are referred as C, complex numbers being calculated are referred as Z (Z1, Z2, … Zn).
As you can see from the Mandelbrot law, the square of a complex number must be evaluated. How can we do that?
From algebra, we know that the square of the quantity (a+b) is:
(a + b)^2 = a^2 + b^2 + 2ab
Well, we just have to apply this rule to a complex number. So:
Z^2 = (a+bi)^2 = a^2 – b^2 + 2 a b i
Why -b^2? The minus sign is there because the square of i is -1.
As we know, we must compute a set of numbers Z0, Z1, Z2, Z3… Zn, and see if the values obtained go to infinity or not. How many values are we supposed to compute? Just a limited number of values, as we will see in a moment.
But, how can we evaluate a quantity such as 4.5 + 6i for instance? We must calculate its module. First of all, we must take the real part and the imaginary part of the number. As it can be easily guessed, the imaginary part is the one along with the i. So, given the number 4.5 + 6i, we have:
Real part = 4,5; imaginary part = 6
Given the generic form of the square of a complex number:
Z^2 = a^2 – b^2 + 2 a b i
the real part is: a^2 – b^2
and the complex part is: 2 a b
This is quite important to notice, as we will use the above formula in the code.
Given the real part and the imaginary part of a complex number, the module can be evaluated as follows:
module = sqrt(real_part^2 + imaginary_part^2)
It can be shown that if the module of a certain Zn ever gets bigger than 2, then the point doesn’t belong to the Mandelbrot set. That allows for a faster computation of the fractal. Still, if after a fixed number of iterations, the module of Zn keeps being less than 2, then the point belongs to the Mandelbrot set.
Instead of evaluating the square root in the computation of the module, we can just take out the square root and test the remaining quantity for the value 4 (2^2, of course). So, we can test if a point belongs to the Mandelbrot set or not by using the following test:
real_part^2 + imaginary_part^2 < 4 ?
Now, where are we going to take the points to test? Given a part of a plane, we must take each point of it, and test it with the Mandelbrot laws. We just have to take the coordinates of each point, and use that to make the iterations.
As we need to work with complex numbers, we are going to take the test points from the complex plane. The complex plane is a plane made up of two axis.
- one axis, is the real axis, made up of real numbers;
- the other axis, is the imaginary axis, made up of imaginary numbers.
Those axis have a 90 degrees angle between them, just like conventional x and y axis.
So, a complex number in the complex plane can be found by using its real part and its imaginary part as coordinates.
To summarize, before starting the iterations, a C point in the complex plane is taken, and its coordinates are used for the computations. A complex number is taken for each point of the plane. Of course, a plane has infinite points. We just take a number of points in the complex plane that equals the numbers of points available on the screen, and we test each of them. So, we must take a set of “input real parts” and a set of “input imaginary parts”. Those will make up the complex numbers to be tested.
Mandelbrot fractals can be drawn by using the complex numbers from -2-1.1i to 0.5 + 1.1i. So, the first calculation will be performed by using the numbers -2 for the real part and – 1.1 for the imaginary part as input. Then, the next c will be obtained by incrementing the previous real part and imaginary part by a small amount.
A PC implementation using Freebasic
In order to obtain a fast computation of the fractal with simple code, I have developed the following simple program using Freebasic, on the PC.
ScreenRes 800,600,32 Dim as integer w, h dim as double r_squared, i_squared, real_part, im_part, im0, real0, module_squared dim as double real_start, real_end, im_start, im_end dim as double dx, dy dim as integer x, y, k, is_mandelbrot ScreenInfo w, h real_start = -2.0 real_end = 0.5 im_start = -1.1 im_end = 1.1 dx = (real_end-real_start)/w dy = (im_end-im_start)/h for y = 0 to (h-1) im0 = im_start + y * dy for x = 0 to (w-1) is_mandelbrot = 1 real0 = real_start + x * dx real_part = 0 : im_part = 0 : r_squared = 0 : i_squared = 0 for k = 1 to 250 r_squared = real_part^2 - im_part^2 i_squared = 2 * real_part * im_part real_part = r_squared + real0 im_part = i_squared + im0 module_squared = real_part ^ 2 + im_part^2 if module_squared > 4 then is_mandelbrot = 0: goto escape next k escape: rem plot point if is_mandelbrot = 0 then pset (x,y), rgb (255,255,255) if is_mandelbrot = 1 then pset (x,y), 0 next x next y sleep
A fixed number of iterations has been used. 250 iterations are more than enough to obtain a detailed shape of the fractal. The above code sets the points of the Mandelbrot set to black, and the other points to white. The program is able to create the following image.
As there is no provision for complex numbers on most computer languages, the real part and the imaginary part of each complex number must be treated on its own, using the above formulas. The code makes use of the following numbers:
- real_start, real_end: those numbers are the limits of the real axis. Those are the real parts of the starting complex test number and of the final complex test number.
- im_start, im_end: those numbers are the limits of the imaginary axis. Those are the imaginary parts of the starting complex test number and of the final complex test number.
- dx and dy are the increments used to obtain the various numbers to test. Those are computed taking into account the screen resolution in use, so that each point of the screen will be tested.
- im0, real0: the imaginary part and the real part of the complex number that is currently being tested.
- r_squared, i_squared: those are the real part and the imaginary part of the squared complex number (those are evaluated on each iteration). They are NOT the square of the real part and the square of the complex part!
- real_part, im_part: those are the real part and the imaginary part of the Z_n+1 complex number, evaluated according the Mandelbrot law. Again, those are evaluated on each iteration as well.
- module_squared is the square of the module of the complex number Z_n+1. This is used to test if any C point in the complex plane belongs to the Mandelbrot set or not.
The above program basically consists of three nested loops. The first two loops are used to take points from the complex plane. The inner loop performs the computations to draw the Mandelbrot set. If this loop is not escaped, the point belongs to the Mandelbrot set, otherwise it doesn’t. Of course, the loop is escaped whenever the square of the module of Z_n+1 gets bigger than 4. That means, if that is the case there is no need to compute other values: after that point, the Zn values will keep growing and will go to infinite (a theoreme states this).
The instruction ScreenRes is used to set the screen resolution. The rgb function is used to set the color of the point to be plotted. Given the red, green and blue components, it evaluates the matching color value. White is obtained by the values 255, 255, 255, which are the maximum level of each color component. Please notice that if we had used the instruction Screen [number] to set the screen resolution, the rgb function wouldn’t have worked properly. While the rgb function is of little use for this black and white version of the program, it will be of help for the color version.
Mandelbrot set in colours: Freebasic and Commodore 64 BASIC V2 programs
Someone on the web managed to reverse engineer the colors used on the Wikipedia images of the Mandelbrot fractal. I used those colors in the following Freebasic program, which shows off some nice color shapes just outside the Mandelbrot set.
Some points may suddenly “diverge”, while others may take more time to “diverge”. That means, the complex number C matching a given point on the screen, may require a different number of iterations before the square of the module of the number Zn_+1 becomes greater than four. If we select the color of the point according to the number of iterations required to identify it as a non-Mandelbrot set point, we can obtain nice shades. So, the colours that we can admire on fractal images are actually points that do NOT belong to the Mandelbrot set.
The following Freebasic program uses the MOD function to set the color of non-Mandelbrot points. That way, given the colour palette of the integer array “mapping”, points that diverge faster are coloured darker, points that diverge slower have lighter colours instead.
ScreenRes 800,600,32 Dim as integer w, h dim as double r_squared, i_squared, real_part, im_part, im0, real0, module_squared dim as double real_start, real_end, im_start, im_end dim as double dx, dy dim as integer x, y, k, is_mandelbrot, iterations, color_index dim mapping(15) as integer mapping(0) = Rgb(66, 30, 15) rem brown 3 mapping(1) = Rgb(25, 7, 26) rem dark violette mapping(2) = Rgb(9, 1, 47) rem darkest blue mapping(3) = Rgb(4, 4, 73) rem blue5 mapping(4) = Rgb(0, 7, 100) rem blue4 mapping(5) = Rgb(12, 44, 138) rem blue3 mapping(6) = Rgb(24, 82, 177) rem blue2 mapping(7) = Rgb(57, 125, 209) rem blue1 mapping(8) = Rgb(134, 181, 229) rem blue0 mapping(9) = Rgb(211, 236, 248) rem lightest blue mapping(10) = Rgb(241, 233, 191)rem lightest yellow mapping(11) = Rgb(248, 201, 95) rem light yellow mapping(12) = Rgb(255, 170, 0) rem dirty yellow mapping(13) = Rgb(204, 128, 0) rem brown 0 mapping(14) = Rgb(153, 87, 0) rem brown 1 mapping(15) = Rgb(106, 52, 3) rem brown 2 ScreenInfo w, h real_start = -2.0 real_end = 0.5 im_start = -1.1 im_end = 1.1 dx = (real_end-real_start)/w dy = (im_end-im_start)/h for y = 0 to (h-1) im0 = im_start + y * dy for x = 0 to (w-1) is_mandelbrot = 1 real0 = real_start + x * dx real_part = 0 : im_part = 0 : r_squared = 0 : i_squared = 0 for k = 1 to 250 r_squared = real_part^2 - im_part^2 i_squared = 2 * real_part * im_part real_part = r_squared + real0 im_part = i_squared + im0 module_squared = real_part ^ 2 + im_part^2 if module_squared > 4 then is_mandelbrot = 0: iterations = k : goto escape next k escape: rem plot point color_index = iterations mod 16 if is_mandelbrot = 0 then pset (x,y), mapping(color_index) if is_mandelbrot = 1 then pset (x,y), 0 next x next y sleep
The MOD function just returns the integer remainder of a division. The remainder of the division between the number of iterations performed before escaping the loop and the numbers of colours available/chosen can be used to select the color of each non-Mandelbrot point.
The above program creates the following image.
Despite being written for clarity rather than for speed, the above program only takes a little amount of time to be executed, even on an old Pentium 4 PC. The pset instruction doesn’t provide excellent performances, but it can be replaced with other instructions. I have used it as I am quite new to Freebasic. Still, other languages may offer better performances, but I have chosen Freebasic as it is quite similar to the BASIC syntax of old BASIC interpreters.
On the Commodore 64, a computation of such a fractal will not be fast by today’s standard. If we want to use simple code, BASIC is what we may look at. The program will be quite slow, but still, we can take advantage of the warp mode of the Commodore 64 X64 VICE emulator. Machine language implementation is possible, but quite complex.
To make the computation not too slow, we may just use screen characters as points for the fractal image. That means, we only have to compute 1000 points, as we have 40 columns and 25 rows on the C64 text screen. But, that sadly means we will have a low resolution 8×8 image.
Anyway, we can have a color advantage. Each 8×8 character cell can be freely set to any of the 16 colours available. And that helps for shading. Of course, even higher resolution modes allow for quite good shadings (especially if using dithering), but this low-res mode at least offers an easy way to get many colors on screen.
We can use the character with ASCII code 166 to create some sort of “mask”. That will somehow trick our eyes a little bit, and the low-res image will appear more detailed than what it really is. So, we can fill the screen with this character, then POKE the screen memory according to the above calculations.
Another trick to smooth things up is to use some sort of “interlacing”. Since we have an 8×8 image, we can “split” each 8×8 pixel character by moving the screen (using the scrolling registers). If we move the screen by 4 pixels in both X and Y directions, then we reset it to its original position, then we move it again and so on, quickly alternating two different images, we somehow manage to obtain an apparent improved interlaced resolution. But, as the movement is really coarse (4 pixels!), we will see the image moving considerably. But still, we will have a rough idea of how interlacing works.
The following BASIC program uses the above techniques. On a real Commodore 64, the fractal takes some minutes to be drawn. As we have said, the warp mode of the VICE emulator really helps here.
I haven’t crunched the program much, as the speed would not improve noticeably. So, I have kept the program quite readable.
DOWNLOAD: Mandelbrot set (Commodore 64 BASIC V2)
Note: a faster version can be found here.
After the fractal has been drawn, you can hit the spacebar to interlace it. Of course, the image will be very unstable, and it is advisable to watch the screen by a 2 meters distance. But, despite the big movement, you should notice that this rough interlacing somehow manages to improve the smoothness of shadings and to provide a better apparent resolution.
The interlace routine has been coded in BASIC as well. It makes use of the VIC-II scrolling registers, and fast writing to these registers is provided by using the “PRINT on memory technique”. Synchronization with the rasterbeam is achieved by using the WAIT statement.
The image is best viewed on a small screen (a little VICE window is OK). I have tested the interlacing on a real Commodore 64 using a CRT TV set. The image needs proper color adjusting, otherwise it nearly seems to explode. The second version of the program, “mand.less.flick”, uses a reduced color palette with “soft” colors, thus reducing the movement of the image during interlacing. That offers a better result on a CRT video screen. However, this is just a quite strange experiment of mine. But, I have read some interesting posts dealing with fractals on the Commodore 64, and I just wanted to dig those things out a bit. Many thanks to the authors of the posts linked below. As for my rude interlacing technique, I took the idea from a demo I saw many years ago, it showed some 8×8 animations but I can’t remember its name any longer.
References:
- http://semioriginalthought.blogspot.it/2012/04/how-to-simply-create-mandelbrot-set-on.html
- http://grahamsneurons.blogspot.it/2013/11/commodore-64-mandelbrot-in-color.html