In the first part of this article I have presented a Commodore 64 program to draw the Mandelbrot set coded in assembly language. Some explanations dealing with signed arithmetic have been provided. Now, we are going to talk about fixed point math, that is, using integers to simulate decimal numbers.
Fixed point math – using an integer offset
Suppose we want to evaluate the sinus function of a number, then we want to compute the value y = sin(x) * 34.
As the sinus function can return decimal values from -1 to 1, problems arise. Since we can only store numbers from 0 to 255 in a byte, how can we represent small steps such as 0.005 and so on?
For instance, if X equals 28, then:
sin(28) = 0.27 (two decimals value, other decimals
are discarded, we are using
a two digits offset)
Y = sin(28) * 34 = 0.27 * 34 = 9.18
How can we represent the number 0.27 in memory? We can multiply it by an integer offset, let’s say 100. So, sin(28) becomes 27.
Now, Y can be evaluated as follows:
Y' = sin(28)*34 = 27*34 = 918 Y = Y'/offset = 918 / 100 = 9.18
It’s important to remember to divide by offset after multiplication: it’s the only way to preserve the information of original decimal numbers. We have already dealt with this kind of math for trig functions in this article.
Now, suppose we want to evaluate the expression Y = sin(34) * sin(78)
Again, we must find the matching integer values using the same offset for both decimal values. So:
S1 = sin(34) * off = sin(34) * 100 = 55 S2 = sin(78) * off = sin(78) * 100 = 97
As both numbers have been multiplied by 100, to properly evaluate the product, everything must be kept in balance. So, we must divide the product by 100*100, the square of the offset.
Y' = S1 * S2 = sin(34)*sin(78) = 55 * 97 = 5335 Y = 5335 / (100 * 100) = 0.5335
The exact value is 0.5469, which is quite close to our approximate result. Of course, we can use an offset greater than 100 and get a more accurate result. For instance:
S1 = sin(34) * 10000 = 5591 S2 = sin(78) * 10000 = 9781 Y = (S1 * S2) / (10000 * 10000) = 54685571 / (100000000) = 0.54685
Now, the result is quite close to the exact value. So, the bigger the offset, the bigger the accuracy we obtain. But, larger offsets do require bigger integer numbers to deal with. So, we just have to select a good value for the offset, so that we get a good accuracy without slowing down calculations too much (larger numbers will more than likely result in slower multiplications).
Sometimes, we just don’t need a final result with decimal numbers. For instance, when computing a Mandelbrot fractal, we don’t care about knowing the values. We just want to compute a value and see if it is bigger or not than a constant value. Let’s see this in practice.
Suppose we want to compute the product:
R*I
Both real and imaginary part are signed decimal numbers. So, we must use an offset for both numbers. On the program, an offset value of 64 has been used. This number is a power of two. This is important, as it allows us to divide by offset very fast using shift and rotate instructions.
Let’s put that R = 1.51 and I = 0.78. Then
R' = R * offset = R * 64 = 96 I' = I * offset = I * 64 = 50
Now, since we have an iterative process, and R* I is one of the computed values used to create the new R and I for the next iteration, then this product must be expressed on an integer form. So, we just divide by offset, and not by offset * offset:
R * I = (R' * I') / offset = (96* 50) /64 = 75
This way, the new I = 2 * R * I + i0 can be created without any issue.
On the formula used by the program we also have additions. Those do not require any offset to be applied, as everything is in balance when adding.
The integer math used in this program is as follows (Freebasic code for simplicity):
for y = 0 to (h-1) im0 = (im_start + y * dy)*offset for x = 0 to (w-1) is_mandelbrot = 1 real0 = (real_start + x * dx)*offset real_part = 0 : im_part = 0 : r_squared = 0 : i_squared = 0 for k = 1 to 250 r_squared = (real_part * real_part) / offset i_squared = (im_part * im_part) / offset if (r_squared + i_squared) > (4 * offset) then is_mandelbrot = 0: goto escape im_part = (2 * real_part* im_part)/offset + im0 real_part = r_squared - i_squared + real0 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
Notice that the limit for I^2 + R^2 is no longer 4, but 4*offset!
As for the values i0 and r0, the assembly program uses tables. These are generated using a BASIC program, multiplying decimal values by the offset and converting the result to integers, then into two’s complement form.
Please note that values in the tables are 8 bit signed integers. The program converts them into 16 bit signed integer. As the code shows, the process is quite straight. As the offset is only 64, 8 bit signed integers are enough, but in the program these values must be used along with 16 bit signed values. That’s why such a conversion is needed. Still, using the 8 bit form in tables makes accessing to them easier.
Speeding things up – using tables for squares. Using an 8 bit and a 16 bit multiply routine
To store the squares of all 8 bit numbers, only 512 bytes are required. Since we need to compute the terms R^2 and I^2, and since numbers are quite slow and most of the time they fall in the 8 bit range, we can use tables to evaluate those squares, and avoid using the slow bit-shifting multiply routine.
The only problem is that we have to recognize 8 bit numbers and use the “square with table” routine only with them. A 16 bit squares table would be huge so that’s not an option – we just have to use this approach.
As you can see from the code, the high byte is checked to see if a number is 8 bit or not. On 8 bit signed numbers stored in 16 bit form, the higher byte must be either 0 of $ff. Then, if it’s negative, its module must be evaluated before squaring it. The whole process is repeated for both I and R.
Taking the square it’s just a matter of reading two bytes, the high byte and the low byte, from the squares tables (there are two tables for low bytes and high bytes). Square values in the tables are already divided by offset. That saves a division in the code.
If a number is 16 bit, then it is squared by using the 16 bit multiply routine.
Likewise, numbers are checked if they fall in the 8 bit range before performing multiplications. If they are 8 bit numbers, the 8 bit multiply routine is used. Otherwise, the 16 bit one is used. Again, the high byte of 16 bit numbers is checked to see if they can be treated as 8 bit signed numbers.
In the 8 bit multiplication, the formula R* I/32 is used. That’s just the same as 2 * R * I /64, but with an added advantage: R is smaller than 2*R. That means, the 8 bit multiplication can be used more often. As a disadvantage, there’s a little accuracy loss that results in two fake Mandelbrot points (this has been discussed in the first part of this article).
One note: when R is multiplied by two with ASL, it is directly changed with no temporary value. That is correct, since after that, R is just reassigned using the values I^2 and R^2 computed much earlier, so no need to preserve R this time.
Multicolor points plotting
Plotting on a multicolor bitmap screen is not much different than plotting on a hires bitmap screen. We just have to plot a couple of horizontal hires points at a time, for each multicolor “big” point (bit pair).
We just have to decide whether each bit in the bit pair must be either 0 or 1. When we switch to multicolor mode, each pair will have a different color. Each bitpair is assigned a color according to the following rules.
[0][0] = background color ($d021) [0][1] = higher nibble in video matrix memory [1][0] = lower nibble in video matrix memory [1][1] = nibble in color RAM memory
The Mandelbrot fractal program uses two tables to define bitpairs. One table contains the first bit. The other table contains the second bit. By looking at these tables, the code decides wheter to turn on a bit or not.
For instance, if the value taken from color_table1 is 0, and the value taken from color_table 2 is 1, then the bitpair is 01, and the matching color is the one in the higher nibble in the video matrix memory.
The color of diverging points is selected depending on the number of iterations.
Please note that colors in the video matrix and in the color RAM are all set the same. Despite each 8×8 cell can actually have different colors, this algorithm just requires random colors for each multicolor point. That just restricts our choice to only four colors. Please also note that the background color would be common anyway.
Those mathematical drawings just make severe restrictions for colors. Art drawings enjoy more colors instead, since the artist can adapt the drawing so that different colors can be selected for any 8×8 cell, virtually allowing to use all 16 colors on screen.
An early version of the program
I think that I have covered most of the topics related to my simple Mandelbrot fractal drawing program, so the code should be quite clear. Still, I am going to provide a simpler code that should be more suitable for studying purposes.
The following is the 8×8 version, plotting coloured chars on the screen. It is quite similar to the Commodore 64 BASIC version presented earlier, so you can use that BASIC version to study the assembly code.
Download: Mandelbrot early version (executable)
Download: Mandelbrot early version (source code)
The program is quite slow. Based on that code, I coded a bitmap version. That version required as much as 24 minutes to draw the fractal with only 20 iterations. Newer versions turned out to be much faster – and there is still plenty of room for optimization. Still, even simple optimizations like the ones presented in these series do allow for better performances.
