Software Notes for Graphing the Mandelbrot Set by George Taylor 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 not to 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. If abs(z)=2 the absolute value of z will not change in subsequent iterations, and is therefore 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 3d 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 heights as black and all other heights as white, EG 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 repeated for decreasing heights. You can also 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 ysize=(cih-cil)/height for y=0 to 199 ci=y*ysize+cil 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 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. 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 ie 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: zr=0 zi=0 zrsquared=0 zisquared=0 for i=1 to maxiter newzr=zrsquared-zisquared+cr newzi=2*zr*zi+ci zr=newzr zi=newzi zrsquared=zr*zr zisquared=zi*zi zsquared=zrsquared+zisquared if zsquared>4 then ... Also many 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. Since with scaled arithmetic, a test of abs(x)>2 can be done with: (make d0 positive) and.w d0,$c000 bnz true: 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 integer format can be used for speed. Using a scaling factor of 8192, IE 1.000=8192 and 7.99=65535, a 16 bit word size will have 3 bits to represent the integer part and the remaining bits representing the fractional part. 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 8192, in a 16 bit word without ever causing an overflow error. With some care, a scaler of 16384 (where 1.000=16384) 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 16384 there are 2 bits in the integer part, the rest in the fractional part, and the range that can be represented is +- almost 4. 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 real 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 nonmember 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 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 abszsquared>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 for 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.