Software Notes for Graphing the Mandelbrot Set

by George Taylor, (Version 1.01 Dec/91)


From: arun@ac.dal.ca
Newsgroups: alt.graphics
Subject: How to plot Mandelbrot set
Date: 18 Dec 91 05:29:41 GMT

The Mset is defined by the recursive equation Z = -Z^2+C, where Z and C are complex values. The starting value for Z is 0. If for some C, Z approaches infinity, C is considered to not be a member of the Mset. That is, with Z starting at zero, and the operation Z = -Z^2+C performed over and over forever, if Z at some point approaches infinity then C is not a member of the Mset.

Fortuneately it can be shown that if, after an iteration, ABS(Z) has become greater than 2, the value of Z will increase towards infinity from then on. Also, it can be shown that after one iteration, Z=C. Therefore, you need only test values for C where ABS(C)<=2, that is where the real and imaginary parts of C vary from -2 to 2 inclusive.

It is impractical to carry out the iteration forever to determine membership, so an iteration limit can be established. Even with powerful computers, an iteration limit of 1000 is considered satisfactory to correctly determine a points' membership. However, if after this iteration limit ABS(Z)<=2, it can only be assumed with high probability, but not proven, that the associated point is a member. If the loop exited prematurely due to ABS(Z)>2, that point is not a member.

The number of iterations performed before the loop is exited due to ABS(Z)>2 is called the height or potential associated with a point C. Those heights equal to the iteration limit are assumed to be members. For an interesting display, the varying heights may be plotted as a z-axis value, leading to a 3-d image, or colored according to some arbitrary scheme at the point C. Traditionally those points assumed to be members are always shown as black pixels.

One scheme for showing non-member heights, called monochrome banding, works by plotting ((iteration limit-height) mod 2)=0 heights as black and all other heights as white, e.g. if iteration limit=1000, height 4 pixels are black and height 5 pixels are white. If iteration limit=25, height 15 pixels are black and height 10 pixels are white. In other words, the member pixels are always black and the pixels of height (iteration limit-1) are white, with the black-white pattern repe

You can also try using different modulos and more colors. For example the first 16 heights could be colored from white in darkening shades of grey to black. Another variation which can show the heights more clearly is taking the logarithm of the height before applying the coloring scheme. This is because the heights tend to rise sharply as you approach a member point. Areas near member points are referred to as boundary areas, or the edges of the Mset.

Taking the information given so far, here is a program to plot the Mset with monochrome banding on a 320 by 200 pixel screen:

INPUT "Initial Real value of C",CRL
INPUT "Final Real value of C",CRH
INPUT "Initial Imaginary value of C",CIL
INPUT "Final Imaginary value of C",CIH
INPUT "Iteration limit",MAXITER
WIDTH=320
HEIGHT=200
XSIZE=(CRH-CRL)/(WIDTH-1)
YSIZE=(CIH-CIL)/(HEIGHT-1)
FOR Y=0 TO 199
  CI=CIH-Y*YSIZE
  FOR X=0 TO 319
    CR=X*XSIZE+CRL
    REM Now CR, CI are scaled
    REM Do Z<-Z^2+C iteration
    MEMBER=TRUE
    ZR=0
    ZI=0
    FOR I=1 TO MAXITER
      NEWZR=ZR*ZR-ZI*ZI+CR
      NEWZI=2*ZR*ZI+CI
      ZR=NEWZR
      ZI=NEWZI
      Z=SQR(ZR*ZR+ZI*ZI)
      IF Z>2 THEN
        MEMBER=FALSE
        EXIT  (exits the FOR loop)
      ENDIF
    NEXT I
    IF MEMBER=TRUE THEN
      PLOT X,Y
    ELSEIF ((MAXITER-I) AND 1)=0 THEN
      PLOT X,Y
    ENDIF
  NEXT X
NEXT Y
END

Where the function PLOT X,Y will color a pixel black, with the screen initialized to white pixels, and point 0,0 being the upper-left corner.

Speed considerations

It should be noted that the test of the absolute value of Z ABS(Z)>2 is eqv. to SQR(ZR*ZR+ZI*ZI)>2 is eqv. to (ZR*ZR+ZI*ZI)>4, i.e. squaring both sides eliminates the need to do a square root. Also, the values ZR*ZR and ZI*ZI can be saved for use in the next iteration, where they are used in the calculation of NEWZR. Finally, the temporary values NEWZR and ZEWZI can now be eliminated by reversing the order of assignment of ZR and ZI:

ZR=0
ZI=0
ZRSQUARED=0
ZISQUARED=0
FOR I=1 TO MAXITER
  ZI=2*ZR*ZI+CI
  ZR=ZRSQUARED-ZISQUARED+CR
  ZRSQUARED=ZR*ZR
  ZISQUARED=ZI*ZI
  ZSQUARED=ZRSQUARED+ZISQUARED
  IF ZSQUARED>4 THEN
  ...

Some points can be prematurely eliminated by testing ABS(ZR)>2 or ABS(ZI)>2; that is, if either of the above tests are true, then you know that the ZSQUARED>4 test will also be true. On a small percentage of points this will save doing the ZR*ZR and ZI*ZI calculations by prematurely exiting the loop. Even with the extra instructions there may be a savings in time.

Let the magnification factor of a graph of the Mset be defined as 1.0 if it is graphed with CR and CI from -2 to +2. Hence the width and height of the resulting graph is 4 units. At higher magnifications, smaller areas of this initial graph are shown on the screen. The higher the magnification, the higher MAXITER must be to show accurate results. For example at a magnification of 1.0 on a 320x200 screen, a MAXITER of only 50 can give satisfactory results. At higher magnifications accordingly higher iteration limits must be used.

When the calculation is performed in machine language, a scaled signed integer format can be used for speed. Using a scaling factor of 4096, ie 1.000=4096 and 7.999=32767, a 16 bit word size will have 3 bits to represent the integer part, 1 bit for the s part, that is a binary number of the form SIII.FFFFFFFFFFFF where S is the sign (S=0 means positive), I is an integer bit, and F is a fractional bit. With scaled arithmetic, a test of ABS(ZSQUARED)>4 is done in assembly with:

    (make d0 positive)
    and.w d0,$4000
    beq loop ;d0<=4

It can be shown that with CR and CI limited to +-2, the ranges of the other variables are:

    NEWZR and NEWZI from almost -6 to almost +6
    ZR^2+ZI^2 from 0 to almost 8
    The values of CR and CI that cause the minimum ZR^2+ZI^2 value are:
    CR and CI are 0
    The values of CR and CI that cause the maximum ZR^2+ZI^2 value are:
    CR and CI are almost 2
    The values of CR and CI that cause the minimum NEWZR and NEWZI values are:
    CR and CI are 0
    The values of CR and CI that cause the maximum NEWZR and NEWZI values are:
    CR is almost +-2, CI is 0

Therefore all numbers involved in calculating the Mandelbrot set can be performed with the scaler of 4096 in a 16 bit word without ever causing an overflow error. Note that you cannot calculate the points C=+-exactly 2.

With some care, a scaler of 8192 (where 1.000=8192) can be used. In this case any overflow in any operation can be taken as a premature indicator of the ZSQUARED>4 test. That is, if any operation overflows, the associated point is not a member. Note that with a scaler of 8192 there are 2 bits in the integer part, the rest in the fractional part, and the range that can be represented is +- almost 4. The binary format is: SII.FFFFFFFFFFFFF.

Plotting can be further sped up by observing that the Mset is symetrical about the x-axis. Any positive CI value can be reflected to it's negative counterpart and plotted if that point is on the screen. Of course when only the upper 1 or 2 or lower 1 or 2 quadrants are being displayed this technique cannot be used, since the reflected symetrical part is not in view.

The stability condition can be used too to exit the loop. If the value of Z remains the same for two iterations in a row, it can be shown that it will remain the same from then on. If this value passes the membership test, further iterations need not be done. Note that the actual values at this point may be changing at some decimal point further on, but within the scope of the precision you're using, it is accurate.

The method of edge detection can also be used. All member points are connected together, that is there are no set of member points which form an island, separate from another group of member points. Hence if you can outline an area whose border contains only member points, all points within that border must be member points. this is true in the real Mandelbrot set, but in the approximated version you see on the computer screen, there may be islands of assumed member points which may be connected by tendrils smaller than your current pixel size. But the principle can still be used - each apparent island on the screen can be outlined.

You can implement this technique by plotting the Mset as per usual, but when encountering a member point, use edge detection to outline the member set area. once outlined, all points within that area can be colored black. Edge detection works like this: starting at a pixel x, y which is a member, test all 8 surrounding points for membership. When one is found make that your new base point, and color it black. Continue until you have traced your way to your original starting point. If you should encounter the edge of the screen, your surrounding points are limited to fewer directions. Finally, do a flood fill of the interior.

There are a number of other optimizing techniques, such as the distance estimator, which will not be covered in this document. However you can feel satisfied in writing the fastest possible loop for the membership test knowing that no further techniques can give the same result. The distance estimator technique for example can speed up graphing of non-member points, but the graph won't look the same as the standard technique, and usually the standard technique is used for the remaining pixels that are very near the boundary and for the members.

The membership test loop is essential in any Mset plotter regardless of what fancy techniques are used.

Testing your program

It can be shown that c=0+0i is a member. For a speed benchmark, I suggest timing the membership test routine with this point for an arbitrary number of iterations to find an iterations/second value. This number will be indicative of the plotting speed of the basic graphing program. You can't precisely predict plotting times even if the average height per pixel of a certain range of C were known, as multiplication of numbers on a microcomputer depends on the number of 1 bits in the multiplier. Yet this can't be taken as a lower boundary on the time if any sort of pre-ZSQUARED>4 testing is implemented, as this could save two multiplies and an addition on some pixels. For accuracy testing, you can try the point C=.25+0i which can be represented exactly in scaled integer form. This point will cause a Z value approaching 5+0i. After 255 iterations you should reach a value of .496188090413988+0i. My own testing has shown that for MAXITER<=256, a carefully written routine will be accurate to the last bit. In other words the maximum magnification possible for a particular scalar can be plotted without error.

For comparison I have compiled a (very) short list of benchmarks from various computers.

Commodore 64 with 256k memory expansion which contains a times table and a squares table, does about 1100 iterations/second. Amiga 2000 does 128700 iterations/second.

Revision 1.01 Notes

Text prettyed up and typos corrected. Changed example program so rightmost pixels represent CRH and top row pixels represent CIH. Flipped the graph to be consistent with a plot routine having 0,0 as upper-left corner. Eliminated swap variables in second example. Fixed explanation of scaled integer format to represent 2's complement signed integers. Problems: a test value should be given for an actual calculation in 16 bit arithmetic. The test value given now is the ideal result. Another test point with an imaginary part<>0 should be given. A table of maximum magnification for various precisions should be given. Table of benchmarks can be updated. Accuracy of 16 bit result compared to ideal result should be checked. Number of iterations, points plotted should be given for a few sets of input parameters. Insert the IBM CGA, Amiga, and C64 BASIC example programs, 68000 assembler, DSP560001 assembler and generic C.