Random Generation of Fractals John Sharp
Enhancing Mandelbrot Fractals Larry Cobb
Mandelbrot Generators on the Macintosh II Mark Datko,
Fractals Maths and Graphics Simon Goodwin
This article explores programming aspects of John de Rivaz's short listing in issue -1 of Fractal Report. Perhaps the first thing I should say is that I don't know how that code works, or - more precisely - why it produces interesting results when apparently similar formulae do not.
I hope that other readers will write in and illuminate us all, and trust you won't feel I'm wasting your time by trying to encourage Fractal Report with an early article, before the intellectual flavour of the publication is established.
I'm not mathematically minded, although I spent a while in so- called 'scientific' programming - almost exclusively working with integer (whole-number) arithmetic. Since then I have written compilers and a display device driver for Z80 and 68000 micros, so I hope I can explain something about speeding up fractal processing even if I can't tell you how to spot a pretty formula off the page.
What I have done - doubtless like many other readers - is entered the program, tweaked it, and tried to deduce things about its behaviour by experimentation. Like John, I ran the program on a Sinclair QL, and a CST Thor XVI QL-compatible. Some of my comments refer to specific QL features, but I'll try to set them in a wider context.
Simple fractals revisited
To start with, a re-tread of John's original program from issue -1. This version has been cut down a little to remove redundant code and hence speed it up.
You can usually get some benefit by trimming code, but it rarely makes a big difference. If you want programs to run quickly on limited hardware it's important to understand the relative efficiency of internal operations performed while the program runs.
Later in this article I'll present and explain further modifications, but for the time being I should explain some of the quirks of SuperBASIC so that you can get the program running on other machines.
10 REMark SuperBASIC fractals by John de Rivaz
20 REMark optimised by Simon N Goodwin 6/12/88
30 REMark Revision 1.1
100 WINDOW 512,256,0,0:INK 7:PAPER 0
110 SCALE 14*10,-10*10,-7*10
120 CLS
130 a=1.01 : b=1 : c=.9
140 x=0 : y=0
150 REPeat loop
160 POINT x,y
170 xx=y-SIGN(x)*(ABS(b*x-c))^.5 : y=a-x
180 x=xx
190 END REPeat loop
200 :
210 DEFine FuNction SIGN(x)
215 IF x<0 : RETurn -1
220 RETurn x>0
230 END DEFine SIGN
Graphic details
The WINDOW, INK and PAPER commands simply tell the program to draw white dots on a black background, using the whole display area of 512 by 256 pixels. The SCALE command indicates that a vertical line of length 140 will fill the Y axis from top to bottom, while the next two co-ordinates put the origin near the middle of the screen, so that the bottom left hand corner of the display is co-ordinate -100, -70.
It's quite rare for micros to work this way, although it can be convenient. Most micros assume that graphics co-ordinates correspond directly to pixel co-ordinates. If this was the case X and Y co-ordinates would be integers between 0 and 511 and 0 and 255 respectively. Such a scheme keeps things simple for the system software, but it can introduce distortion and extra work for the programmer.
Most micros (but not the Apple Mac or Spectrum) push out rectangular pixels, so a displacement of N units gives a different spacing on the screen depending on whether it is horizontal or vertical. Worse still, if you're using a TV display the exact shape of the dots varies depending on the TV standard in use - Amiga users may have noticed that 'circles' drawn by American programs appear foreshortened vertically. This is because the US TV display has 525 lines whereas most European TV systems use 625 lines, foreshortening each point vertically compared with the US system.
The QL uses oblong pixels but translates graphics co-ordinates so that circles appear round regardless of the display mode or TV standard. The machine knows the aspect ratio (shape) of the window, and scales co-ordinates accordingly, adding the origin offset for good measure. If a point or part of a line turns out to be outside the window border, the QL just clips it off without an error message. Floating point co-ordinates allow an enormous range of scales.
This is good news in some ways - it makes it easy to re-scale or move the centre of a drawing, simply by adjusting the SCALE statement, and it appears that the '*10' factors were added later so that the 'outside' of the pattern could be seen. Try removing them to see the centre in more detail.
The snag is that independent X and Y scaling and origin offsets must be applied to each co-ordinate pair, using floating-point arithmetic, before the point can be plotted. Unless you've got a good maths co-processor chip, this arithmetic takes much more time than finding and plotting the point from its pixel co-ordinates.
The cramped QL system ROM makes things slower by treating a point as a very short line, and doing the maths for each 'end' of the line separately. Worse still, a line is treated as a special case of a curve. The result is that the QL draws circular arcs at spectacular speed, but POINT works rather slowly. Later I'll explain how to speed up pixel plotting.
SuperBASIC Quirks
The REPeat ... END REPeat statements form an infinite loop. In pidgin BASIC you'd forget line 150 and write GO TO 160 at line 190.
Line 170 uses the SIGN function, defined later. Most BASICs have a built-in function, SGN, which returns -1, 0 or +1 depending on the sign of the argument. In SuperBASIC we define our own; lines 210 to 230 are not needed in versions of BASIC where you can call SGN directly in line 170. ABS finds the absolute value of its argument.
I have simplified John's original code by removing the temporary variable YY, which was not needed, and re-coding SIGN to take advantage of the fact that a comparison on the QL returns +1 if true and 0 if untrue.
Added colour
This code gives an interesting pattern, but it takes a long time to draw it. The evaluation time is dominated by the square root in line 170, so I decided to try to re-code it to speed things up. John's program takes a square root by raising the number to the power of 0.5. The 'power' operator, '^', in a micro maths package usually works by using logarithms and anti-logs, so it's not the fastest technique. Square roots are a 'special case' that can be found more quickly and directly with a technique like Newton-Raphson, and the QL has a function to work them out that way.
I replaced the '^.5' with a call to the intrinsic function SQRT (SQR in other BASICs) to see if things went any faster. Line 170 became:
170 xx=y-SIGN(x)*SQRT(ABS(b*x-c))) : y=a-x
This seemed to speed things up when I ran the program, but it was hard to be sure... At first the points plotted were just the same as when using '^.5', but later they diverged, slightly at first and then massively, so it was hard to spot distinct phases of drawing because the sequence was different.
This was just the kind of fascinating behaviour I expected of fractals, where tiny rounding errors make an enormous difference to the plotted result. It also proved that the QL maths package is not smart enough to recognise '^.5' as the same case as SQRT, because the results were different.
To explore this difference, I wrote a double version of the program which plotted points in red and green, using red (2) for the '^.5' expression and green (4) for SQRT. The SIGN function is as above:
10 REMark SuperBASIC fractals by John de Rivaz
20 REMark optimised by Simon N Goodwin 6/12/88
30 REMark Revision 2.1
100 WINDOW 512,256,0,0:INK 7:PAPER 0
110 SCALE 14,-10,-7
120 CLS
130 a=1.01 : b=1 : c=.9
140 x=0 : y=0 : x1=0 : y1=0
150 REPeat loop
160 INK 2 : POINT x,y : INK 4 : POINT x1,y1
170 xx=y-SIGN(x)*(ABS(b*x-c))^.5 : yy=a-x
175 x2=y1-SIGN(x1)*SQRT(ABS(b*x1-c)) : y2=a-x1
180 x=xx : y=yy : x1=x2 : y1=y2
190 END REPeat loop
At first points appear in red and then are immediately over-plotted in green. Later the colours diverge, and seem to leave space for one another on the picture. This reminds me of the way that pseudo- random number sequences mesh together if you start them with two different seed values and map the outputs onto the screen. Perhaps a mathematician can see the similarity, or point out where I've made an invalid connection - preferably without drowning us all in Gauss' Law of Quadratic Reciprocity.
This experiment seemed a sucess, so I decided to try other crude techniques to get the code working in colour. The highest QL graphics mode allows four colours; allowing a black background, that left three colours: red, green and white. There seemed to be some symmetry of order 3 about the picture so I naively told the program to plot each subsequent point in a different colour, using a three colour cycle. To re-create this experiment, add these lines to the first listing:
145 colour%=2
155 INK colour% : colour%=colour%+2 : IF colour%>6 : colour%=2
This gave immediate results, with three coloured circles appearing near the middle of the screen. For a while the colours appeared to be scrambled together, but then a clear ring of triangles appeared around the picture, with each triangle in a clear, distinct colour.
Going faster
The next step is to speed things up. The easiest way to speed up an interpreted program is to compile it. A compiler converts a program once and for all into a sequence of machine-code instructions. Normally interpreted BASIc spends most of its execution time finding lines and variables, and relatively little working out results, so acompiler can make a big difference to execution speed. Eight bit micro compilers tend to be disappointing when you compile fractal programs, because the time taken to work out floating point values swamps the interpretation overhead.
Bigger processors - including the QL's 68008 - can multiply and divide integers and work with 16 or 32 bits at a time internally, so their floating-point takes a smaller proportion of the total interpretation time. I used TURBO on the QL, mainly because it's quick and I wrote it.
The final compiled version of the program was 5.25 times faster than the same code when interpreted. This is quite a high factor in view of the device output and floating-point overhead in the program. What's more, compiled code can multi-task, whereas you can only interpret one QL program at a time. It's a pity I don't get royalties!
To speed things up a lot we must remove the implicit maths involved in plotting pixels with POINT. The QL has a BLOCK command for filling a rectangular area. Parameters are all expressed in pixels, so the command does not have to do much translation to find screen addresses, albeit distorted in aspect ratio. We must scale the co-ordinates ourselves. That means more explicit floating-point arithmetic, but TURBO optimises enough to make this worth while, especially if we use factors that are powers of 2, like 4, 8 or 16.
160 BLOCK 1,1,x*8+256,y*8+128,7
This is over 3 times faster than POINT, despite the extra sums and parameters, and even though BLOCK is optimised for wide areas and plots a single point by drawing a short line and trimming the edges!
The last parameter of BLOCK is the colour, so I decided to switch to the QL's lower resolution 8 colour mode, to see if any cycles of order 7 became apparent. These are the extra lines to be added to listing 1 to try this:
110 MODE 8 : REMark SCALE is irrelevant to BLOCK
145 colour%=0
155 colour%=colour%+1 : if colour%>7 : colour%=1
160 BLOCK 2,1,x*8+256,y*8+128,colour%
At first this doesn't seem to give much away, but if you leave the program to run you should see clear, distinct circles of colour appear around the edge of the pattern, after a while.
Last word?
This newsletter must be a co-operative effort if it is to continue beyond the first couple of issues. If this, or other articles stimulate you to think further, or offer criticism, I hope you will find time to write in and share your thoughts. If you think we should be discussing other things, please get us started. Don't be shy, or none of us will learn anything!
c Simon N Goodwin 14/2/89