QB Fractal Tutorial

Fractal Exploration in Quickbasic

The purpose of this tutorial is to explain the basic theory of and demonstrate how to render an escape-time fractal in Quickbasic. It assumes only some experience at coding in QB.

Although the exact definition of a fractal is still an open question, for the purposes of this tutorial a fractal is just a shape that has some degree of self similarity - that when you zoom into it suffienctly and in the right place, you'll find an identical copy. The most commonly known escape-time fractal is called the Mandelbrot Set.

The mantelbrot set

The Mandelbrot Set
The interesting parts of the fractal are labeled on the image above:
A : This circlular area is the limit of the fractal. In the case of the Mandelbrot set, it's two units in radius.
B : This large white area is called the lake. The lake is the only area on screen that is actually part of the Mandelbrot set - all the other points are outside it. The gradient of colour that you see approaching the lake is actually just a popular graphical trick that I'll explain later on in the tutorial.
C : This is the horizontal axis.
D : This is the vertical axis.
E : The region around the lake which is almost in the set is called the boundry. This region is important for two reasons - one, it's where you find smaller versions of the Mandelbrot set, and two, it contains all the interesting seed points for the julias, which I'll talk about later.

Since the Mandelbrot Set is generated on a complex plane, a knoweldge of complex numbers is essential. The complex numbers are simply the a field of numbers containing all the real numbers and additonally, something called the imaginary number i, where i^2 = -1. This allows you to express the square root of negative numbers as a real number times 'i', for example, sqrt(-4) = 2*i.

Complex numbers are useful because they allow you to express solutions to problems such as this : find the roots of the equation x^2 - 4*x +5 (ans: x=2+i and x=2-i). A complex number is usually expressed in the form z = a + b*i, where a and b are both real numbers.

A Complex Plane
The following properties of complex numbers are of interest when considering fractals:

(1) Complex numbers can be usefully represented on a plane, with the real part of the complex number on the horizontal axis, and the imaginary part on the vertical axis.

(2) In the field of real numbers, the |x| operator would return the absolute value of x, i.e. if x were negative, it would return -x, if x were positive, it would simply return x, so that your result is always positive. Another way to look at this is that the |x| operator simply returns the length of x. In the field of complex numbers, this operator is called the modulus and returns the distance of the complex number from the origin.

The modulus Operator

|z| = sqrt(a^2 + b^2)
So for example, the complex number z = 3 - 3i has modulus |z| = sqrt(18) - in gereral you can work out the modulus by the simple formulue |z| = sqrt(a^2 + b^2).

Now, The Mandelbrot set is generated by an iterative function on a complex plane - (Z_(n+1)) = (Z_n)^2 + c, where z_0 is 0+0i, and c is any complex number. For any value of c, the rule is that if as n goes to infinity, then if |z_n|<2, that point is inside the Mandelbrot set. It's called the bailout test.

Take an example: Is 0.5 + i inside the Mandelbrot set?

Well, let's test it: Remember, Z_0 is 0 + 0i:

So z_1 = (z_0)^2 + c
       = 0 + 0i + 0.5 + i
       = 0.5 + i
   z_2 = (z_1)^2 + c
       = (0.5)^2 + i + i^2 + 0.5 + i
       = (0.25) + 2i - 1 + 0.5
       = -0.25 + 2i
   z_3 = (-0.25)^2 + -i + 4i^2 + 0.5 + i
       = (0.25/4) + 0.5 - 4

..etc

You've probably noticed a problem at this stage - a computer cannot possibly test until infinity. The simple solution that's used is to limit the number of iterations used - in other words, to say that if after say 170 iterations, |z_170| is still less than 2, we'll say that the point c is probably in the Mandelbrot set.

We're now one step away from generating this fractal in QB: Iterations.

Essentially the algorithm is as follows:
Say we're going to test up to 170 iterations:

Set the screen mode as 320x200x8, then setup a suitable palette
(the palette should be a gradient between 0 and 170, preferably with black at 0 and a bright color at 170 for the lake.)
Scan across the screen:
for each point on the screen, convert it's (x,y) position to a suitable complex coordinate.
Test whether or not it's in the Mandelbrot set (counting the iterations)
Plot a pixel on the screen using the color corresponding to the number of iterations used:
eg. so if the pixel failed the bailout test after 15 iterations, use color 15, and if it got all the way to 170, use color 170.

The only realy problem here is converting the x,y coordinates to complex numbers: You can find a good transformation by practice, but the one I use for the Mandelbrot set is:
For x between 0 and 319: ((x% / 319) * 4) - 2.5
For y between 0 and 199: (((y% * .625) / 199) * 4) - 1.2

This centers the fractal on the screen, and scales it just right.

It's high time we wrote some code.

This program will display the Mandelbrot set on screen, with coloured iterations:

'$DYNAMIC
SCREEN 13

'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
'Sets up the first 25 palette entries as a gradient from (0,0,0) 
'to (25,48,63). This fast transition will make the outer borders
'of the fractal more distinct
'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
FOR i% = 0 TO 25
  OUT &H3C8, i%
  OUT &H3C9, i%
  OUT &H3C9, (i% / 25) * 48
  OUT &H3C9, (i% / 25) * 63
NEXT i%

'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
'This completes the transition from (25,48,63) to white with the 
'remaining relevent palette entries.
'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
FOR i% = 26 TO 170
  OUT &H3C8, i%
  OUT &H3C9, 25 + ((i% / 170) * 38)
  OUT &H3C9, 48 + ((i% / 170) * 15)
  OUT &H3C9, 63
NEXT i%

FOR y% = 0 TO 199
  FOR x% = 0 TO 319
    '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    'zr! and zi! refer to the real and imaginary parts of z. 
    'They start at 0 + 0i
    '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    zr! = 0: zi! = 0
    '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    'cr! and ci! refer to the real and imaginary parts of c.
    'They depend on the screen location, and are positioned
    'so that the origin is about center, and that the axis 
    'extend about 2 units in all directions
    '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    cr! = ((x% / 319) * 4) - 2.5
    ci! = (((y% * .625) / 199) * 4) - 1.2
    '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    'The iteration loop starts here
    '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    FOR i% = 0 TO 170
      '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=      
      'A reminder of the formula:
      'z_(n+1)=(z_n)^2 + c
      'When considered as distinct real and imaginary parts,
      'we have = (zr+zi*i)^2 + cr + ci*i
      '        = zr^2 + 2*zr*zi*i + zi*zi*i + cr + ci*i
      '        = (zr^2 - zi^2 +cr) + (2*zr*zi + ci)i (as i^2=-1)
      '        = real part         + imaginary part
      'i.e...
      '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=      
      zrnew! = (zr! * zr!) - (zi! * zi!) + cr!
      zinew! = (2 * zr! * zi!) + ci!
      zr! = zrnew!: zi! = zinew!
      '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
      'This is where we do the bailout test. Instead of checking 
      'whether or not sqrt((zr! * zr!) + (zi! * zi!)) > 2, it 
      'makes a lot of sense to simply square both sides and save
      'on a potentially very time consuming operation...
      '=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
      IF ((zr! * zr!) + (zi! * zi!) > 4) THEN
        GOTO exitloop
      END IF
    NEXT i%
    i% = 170
exitloop:
    PSET (x%, y%), i%
  NEXT x%
NEXT y%
'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
'That's all there is to it.
'=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
DO: LOOP UNTIL INKEY$ <> ""

You can see this code without comments here.

The process is exactly the same for other escape time fractals generated from different formulas. Try, for example, using (z_(n+1))=(z_n)^k+c instead - the results are very interesting for non-integer and negative values of k!

Once you've got a working fractal viewer, you'll probably want to use it to explore the fractal. I suggest allowing modification to the scale and positioning of the screen to complex plane transformation in the program. You may find it useful to use a lookup table for this.

Something quite easy to do without a lot of work is a Julia explorer. Julia's are related fractals that were actually discovered before the Mandelbrot set, even though they are generated in almost the same way.

Julias

Four example Julias based on the Mandelbrot set
There are an infinite number of Julia sets, one for every conceivable complex coordinate. They are generated by providing what's called a seed value. How the Julia looks will depend entirely on the seed you pick. There's an interesting relation between the seed values and the same coordinate on the Mandelbrot set - if you zoom the Mandelbrot at the seed point, you'll see features similar to the Julia generated by that seed. In this way, the Mandelbrot set can be considered a map of Julia fractals.

You can generate a Julia set with the code used earlier for the Mandelbrot, just by making one slight change: In the Mandelbrot, z_0 is set to 0 + 0i, and c is the complex coordinate. However, for the Julia, z_0 is set as the complex coordinate, and c is the seed value. Then all you have to do is generate as normal.

You may find some of these seed values interesting:
0.7 - 0.25i
-1 + 0.25i
0.28 - 0.52i
-0.0672403 - 0.647148i

An animated Julia viewer can be seen here.

There are a few simple optimisations applied to the code, but it should be easy enough to follow. One of the main reasons why somebody would want to program something like this is as a programming practice - in fact, there are very few more interesting coding problems than optimising a fractal viewer. As a finishing touch to this tutorial, allow me to suggest some simple optimisations you could proceed with:

(1) The Mandelbrot set is symmetrical around the x axis. If you set the x axis at exactly y=100, you can simply mirror the top image. Of course, this is no good if you plan to build a zoomer, however, it is useful if you want to make a Julia explorer, as Julias are also symmetrical around the x axis.

(2) The lake is the most time consuming part to draw - but you can reduce the amount of time it takes by creating a lookup table of 170 complex coordinates and saving each value of z_n as you go. Then for each z_k such that k>n, check to see if z_k=z_n. If it is, you know that you are definitely in the lake (as it means that the iterative function will continue in a loop), and you can move onto the next pixel. This will slow down drawing of the boundary by quite a bit, but the speed gain on drawing the lake is definitely worth it.

(3) You could substitute the floating point math for fixed point - you'll lose a little bit of accuracy, and you'll have to be careful with keeping track of your scaling, however, it's worth it as integer math is usually much faster.

(4) One classic trick is called the 'bounding box method'. Rewrite your code so that it skips every second pixel and offsets every second line by one, so that there is a mesh of black pixels on the screen. This means that you only draw half as many pixels. You fill in the rest of the screen as follows: Scan over the screen, and when you come to a black pixel, check every surrounding pixel. If they are all the same color, then fill in that pixel with that color. If they aren't, then you'll have to calculate it's color in the normal way. This optimisation may not sound like it gives much of a speed boost, however, when faced with a screen that's mostly lake, it's one of the most powerful.

(5) QB can plot pixels much faster if you poke them to the screen buffer instead of using the PSET command. That's done as following:

DEF SEG = &HA000
  POKE x% + (y% * 320), colour%
DEF SEG

(6) Finally, and perhaps most simply, study your formula and try to use as few operations as you can. One simple suggestion for the Mandelbrot set is to put zr! * zr! and zi! * zi! in temporary variables, and use them both times instead of squaring them twice. This means that'll you'll be an iteration off, but that's easily fixed, and it means half as much math for QB.

That's it. I hope you've enjoyed this tutorial! If I've failed to make myself clear, send me an email at blinkthedarkness@hotmail.com and I'll try and clarify the best I can. Have fun exploring! Hopefully I'll see you next month with another tutorial!

- Terry Cavanagh