GFX Sensations 1/Graphic Sensations - Volume 1.iso/tools/amiga/fractals/rm2.lha/rm2.c

Text File | 1994-01-01 | 33.3 KB | 1,035 lines

/* 80 columns:
123456789 123456789 123456789 123456789 123456789 123456789 123456789 1234567890
*/
/*
 * [e-mail] : Steven C. Den Beste
 *
 * "rm2.c": This is based on the program "retmand" which someone else posted
 * recently. I could not have done this without him, and I hereby announce
 * my gratitude. Nonetheless, I suspect that this is so dramatically modified
 * at this point that the original author wouldn't recognize it.
 *
 * The original program just drew the entire Mandelbrot drawing, using random
 * colors each time. That's all it did - but to do that it had to implement
 * the entire algorithm for CALCULATING it, and all the code
 * for displaying stuff on the Retina. Those are the parts I couldn't have
 * handled. In particular, it used the smart algorithm for calculating
 * the Mandelbrot diagram (the one which optimizes square sections) instead
 * of the stupid algorithm (the one which calculates each point every time).
 *
 * What I DO know how to do is enhance,
 * clean up, and optimize. But to do that I had to
 * perform some modifications. First and foremost: The original code was
 * written for the SAS compiler. What I've got is
 * Aztec C 5.0. Regrettably, they aren't compatible. For a while I tried
 * to maintain compatibility with #ifdef's, but that got so complex that
 * I finally gave up and deleted everything non-aztec, so this code is
 * thoroughly MANXified now. Adapting this back to SAS is probably going to be
 * a pain.
 *
 * The core of any Mandelbrot program is a routine which iterates a
 * rather strange function, waiting until the value of the function
 * exceeds a certain threshold. The value it returns is the integer
 * count of the number of loops it took for this to take place. All the
 * calculations are done in floating point. As a result, Mandelbrot programs
 * have a tendency to take enormous amounts of time unless handled well.
 * One finds implementations using fixed-point integers, but these
 * lead to strange results where the program has a "floor" of detail, for
 * instance. Using floating point prevents this flaw.
 *
 * The innermost loop of this function can iterate literally millions of
 * times depending on the area you are viewing and the dimensions of the
 * screen. Any tiny gains in this function can have significant results in
 * the total compute time.
 *
 * So I dug out my 68882 manual and got to work. It's in ASM now in a separate
 * module (another reason it's probably going to be hard to re-SASify).
 * My code is optimized to
 * take advantage of the '882's capability to perform simultaneous operations.
 * (Don't worry - it works just fine for an '881 - it just takes longer.)
 * The big win is that there are 8 registers out there, and I can
 * store a lot of values permanently in the FPU - effectively,
 * they become "register double" variables.
 *
 * It is possible that  more juice can be squeezed out of this routine, and
 * even a microsecond or two can make a big difference.
 * I invite anyone who is into this sort of thing to try. I'm a programmer
 * by profession, and we tend to avoid ASM like the plague because
 * it is almost always too expensive to debug. (We use it for interrupt
 * routines, task switchers and not a hell of a lot else.)
 *
 * In the mean time I've added a bunch of features. There are now a bunch
 * of parameters you can specify on the command line, as follows:
 *
 * /r number    - specify the "r" (actually horizontal) value of the CENTER
 *          of the screen. Negative is left. Values range from about
 *          minus 2 to plus 2.
 *
 * /i number    - specify the "i" (actually vertical) value of the CENTER
 *          of the screen. Negative is up. Values range from about
 *          minus 1 to plus 1.
 *
 * /w number    - specify the width of the screen in "r" coordinate terms.
 *          The vertical size is computed to maintain an aspect ratio
 *          of 1.2 to 1.
 *
 * I made /r and /i specify the CENTER deliberately. If you can locate the
 * center of a feature, it is convenient to generate the frames
 * of a zoom movie, since all you have to do is change the "w" value. If
 * /r and /i specified the upper left corner, that would be a royal pain.
 * (The names "R" and "I" for these coordinates were established by Mandelbrot
 * himself, and are maintained here for historical reasons.)
 *
 * If you don't specify any of the above three parameters, you get the entire
 * Mandelbrot diagram.
 *
 * /c number    - Choose a color set. Frankly, right now I'm not very happy
 *          with the ones I've chosen. If you don't specify anything,
 *          you get a standard set I think works fairly well. If you
 *          specify "/c 0" you get the random color set the old program
 *          would have given you by default. The others are somewhat
 *          obscure. There is plenty of opportunity here for enhancement.
 *
 * /q          Use a quarter of the screen. Actually a misnomer - this
 *          divides the screen linearly by 4, so it divides the area
 *          by 16. Use this when you're trying to home in on /r and /i
 *          values for interesting features to cut way down on the
 *          compute time.
 *
 * /h          Use a half (that is, 1/4 by area) of the screen.
 *
 * /t number      Raise the top of the drawing. The default is 255. Values
 *          above 255 will cause the color set to recycle. Values can
 *          be specified up to 32767. Please note that depending on the
 *          image you are computing, raising this value can radically
 *          slow down compute time.
 *
 * /d number      Often used in hand with "/t". This takes the output of the
 *          function and divides it by the parameter for every pixel.
 *          In essence it is lowering the color thresholds. For instance,
 *          "/d 3" means that the first color threshold represents the
 *          first THREE step-values, the second represents the
 *          next THREE step-values, and so on. This does not alter
 *          the time to compute the drawing, but it can alter its
 *          image quite a bit.
 *
 * /s filename      If this parameter is NOT present, then when the picture is
 *          completely drawn it will stay on the screen until you hit
 *          the "ENTER" key. If this parameter IS present, then when
 *          the image is complete it will be written out as a 24-bit
 *          PPM file, and then this program will automatically terminate.
 *          This means if you want to create the frames for a zoom
 *          movie, you can set up a batch file and let it run overnight,
 *          or for a week, or for a month, or for ...
 *          WARNING: These files are ENORMOUS. You can chew up disk space
 *          really fast. The file takes 3 bytes per pixel. An
 *          800*600 frame is 1.5 megabytes (sigh). You might want to make
 *          your batch file alternate between creating files with this
 *          program and crunching them with "lha", which can typically
 *          collapse such a file to 100 Kb or less, depending on
 *          image complexity. Alternatively, you can use "ppmtogif"
 *          which won't do as well as "lha" but leaves them in a
 *          displayable form. By definition these files always have
 *          fewer than 256 colors, so "ppmtogif" won't ever tell you
 *          to go use "ppmquant" first.
 *
 * /m threshold      What happens if the Mandelbrot threshold of 4.0 is modified?
 *          Beats heck out of me, but this parameter allows you to
 *          experiment. (My early experimenting shows that values above
 *          about 2.1 all act about the same, and below 1.9 the whole
 *          image turns black.)
 *
 * To get you started, here are some interesting parameter values to try:
 *
 * /r -.93402 /i -.247638 /w .00004 /d 2 /t 512
 *
 * /r -.5554 /i -.5292 /w .006
 *
 * /r -1.2909 /i -.0861 /w .00155
 *
 * /r -.6763 /i -.36742 /w .0003
 *
 * /r -676204 /i -.367276 /w .000017
 *
 * /r -.6762039 /i -.3672763 /w .0000021 /d 2 /t 512    or
 * /r -.6762039 /i -.3672763 /w .0000021 /d 3 /t 768
 *
 * Even with all the optimization, these are not going to be fast drawings,
 * but the results are spectacular (in my own highly biased opinion). These
 * tended to take more than half an hour on my system at 800*600. Your
 * mileage may vary. (My system is a GVP G-force 68040 at 33 MHZ running
 * burst mode, so a 25 MHZ 68030 is probably going to run 4 or 5 times
 * as long.)
 *
 * If you find any interesting coordinates and feel like doing so, please send
 * them to me via email at the address at the top line. Bugs are likewise of
 * interest.
 *
 * Oh, I also added a hell of a lot of comments.
 * There ARE a few comments here which are not mine, but you'll
 * pick up on my style rapidly. Where there's a chance of confusion, I'll
 * include my initials: SCDB.
 *
 * ...and for my next trick: The one obvious need here is a simple way to pick
 * out the coordinate of some object, so as to center and zoom on that object.
 * The obvious way to do that is to use the mouse. Well, I can't. To do that
 * I have to turn on a Sprite, and the part of the Retina includes for that
 * uses a standard SAS include named "utility/tagitem.h", which Aztec doesn't
 * have. I can't find the equivalent structure ("struct TagItem") or define
 * ("TAG_USER") anywhere in the include structure I've got and I haven't got
 * the faintest idea of what they should be. If someone can tell me, I'll try
 * to put it in. In the mean time, the only way to home in on objects is to use
 * "/q" and to tweak the "/r" and "/i" values.
 */

#define UTILITY_TAGITEM_H           /* This disables some nested includes */
                     /* which Aztec doesn/t support. See above. */
#include <intuition/intuitionbase.h> /* Needed to define a major structure. */
                     /* SAS got it from one of the includes */
                     /* files that aztec didn't have, which */
                     /* I deleted. */


#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <exec/libraries.h>
#include <devices/timer.h>
#include <time.h>


#include "retina/retina.h"
#include "clib/retina_protos.h"
#include "pragmas/retina_pragmas.h"

#define MAX_ORBITS 256



/* uncomment the following line to avoid repeating calculations */
/* WARNING: For some reason this creates code that will cause
   stack overflow errors on any image larger than 320x400 or so */

/*
 * SCDB: As far as I can determine, this was a problem the original author
 * encountered early and never pursued. Uncommenting this causes compile
 * errors, as parts of the code which are NOT controlled by ifdefs try
 * to use structure members which ARE controlled by ifdefs on this symbol.
 * At this point, I consider the code controlled by this symbol to be dead,
 * but I invite anyone to try to fix it if they wish. The benefit is further
 * speed increases, at the expense of code complexity.
 */
/* #define NOREPEATCALC */

struct BrotPoint {
  unsigned int x, y;
  double r, i;
#ifndef NOREPEATCALC
  int orbits; /* -2 if untested, -1 if never escapes,
         >= 0 otherwise */
#endif
};

/*
 * SCDB: "threshold" is traditionally "4.0" as defined by Mandelbrot.
 * However, I thought it might be fun to allow it to be modified, so I've
 * made it a variable here whose value can be changed on the command line.
 * It gets loaded by the ASM routine into fp1.
 */
double threshold = 4.0;

int max_orbits = MAX_ORBITS;
int step_size = 1;

int screen_width, screen_height;
int quarter_flag = FALSE;
int half_flag = FALSE;
int h_offset = 0, v_offset = 0;

void Compute_Block(struct BrotPoint, struct BrotPoint);
extern short int Test_Point(double, double);
void my_RectFill(unsigned int, unsigned int,
                    unsigned int, unsigned int);
/* version string */
UBYTE vers[] = "$VER: RetinaMand as modified by SCDB, Dec. 28, 1993";

struct IntuitionBase *IntuitionBase;

struct _xy_RetinaBase *RetinaBase;
struct DefaultScreenInfo screen_info;

/*
 * SCDB: I haven't the faintest idea what this does. It probably doesn't do
 * do anything under Aztec. Anyway, it seems to cause no harm, so I left it in.
 */
int __oslibversion = 37;

/*
 * SCDB: The code depends on the fact that the following number is divisible by
 * 3. If it ceases to be so, this program will crash horribly.
 * "out_buffer" is used to accumulate large blocks of data for writing out to
 * the PPM file if "/s" is used. It is a fact that programs work far more
 * efficiently when they write large blocks rarely than if they write small
 * blocks frequently.
 */
#define BUFF_SIZE 30000
unsigned char    out_buffer[BUFF_SIZE];

/*
 * SCDB: Later on, if we decide to save a file, we're going to have to know
 * what the palette was. The following array will store those values.
 */
unsigned long palette[256];


double dr, di;
struct RetinaScreen *rScreen;

int main(int argc, char *argv[])
{
  int i;
  double rstep, gstep, bstep;
  int c;
  UWORD x_scan, y_scan;
  long    color_value;

  /*
   * SCDB:
   *
   * "red_value" used to be named "foo".
   * "green_value" used to be named "bar".
   * "blue_value" used to be named "baz".
   *
   * Where I work this would get you hung without a trial.
   *
   * These also used to be UBYTEs, but you can get better palette results if
   * they are floating point numbers, since this allows you fractional steps
   * (which round when ultimately converted to RGB values).
   *
   * (I have to admit a bad habit, however: I tend to use "i", "j" and "k" for
   * arbitrary loop counters when they are used very locally. Depending on how
   * old you are, you may have heard of a language named "FORTRAN". Variables
   * whose names began with letters between "I" and "N" were integers, with
   * all other variables being floating point. It was routine to use i/j/k
   * for loop-control variables. However, I never use such names when
   * their scope is large. In this program I use "i", and I'm only using it
   * twice, both times to step through arrays: I use it to scan 'argv' for
   * parameters, and I use it to construct output buffers when writing PPMs.)
   */
  double red_value, green_value, blue_value;
  struct BrotPoint start, end;

  /*
   * The image created by this program used to be forced to the entire
   * diagram. Part of my redesign is to allow a subsection to be passed in
   * as invocation parameters, and to calculate and display just that part.
   * That required replacing a bunch of magic numbers with variables and
   * calculations in the initialization code.
   */
  double left_coord, top_coord, width, height;
  double i_center, r_center;

  FILE *save_file = NULL;

    /*
     * If there are no parameters later for a subsection, use values which
     * create the entire diagram. "height" is calculated to keep the aspect
     * ratio of the screen we opened.
     * These are not the values used in the original program, which
     * discarded the first contour, which is a large circle.
     */
    r_center = 0;
    i_center = 0;
    width = 6;

    /*
     * These values are used if no color parameter is specified.
     */
    red_value = 0;
    rstep = 4;
    green_value = 0;
    gstep = 2;
    blue_value = 0;
    bstep = 12;

    /*
     * Now for something completely new: Invocation parameters!
     * The parameter checking here is not quite as rigorous as it probably
     * should be. Any switch using parameters only checks to see that there
     * are enough entries in "argv" remaining to satisfy its needs. They
     * do not check to see if they make sense. Let's assume some intelligence
     * on the part of the users. There's prudence, and then there's
     * paranoia.
     */
    for (i=1; i<argc; i++) {
        if (strcmp(argv[i],"/c") == 0) {
            /*
             * "/c n" predefined color sets, as follows:
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /c\n");
                exit(1);
            }
            c = atoi(argv[++i]);
            switch (c) {
            case 0:
                /*
`                 * 0: random values, sort of. This is a direct
                 * steal from the original program.
                 */
                srand((long) time(NULL));
                red_value = (rand() >> 3) & 0xFF;
                green_value = (rand() >> 4) & 0xFF;
                blue_value = (rand() >> 3) & 0xFF;
                rstep = ((rand() >> 5) & 0x0F) + 1;
                bstep = ((rand() >> 5) & 0x0F) + 1;
                gstep = ((rand() >> 5) & 0x0F) + 1;
                break;
            case 1:
                /*
                 * 1: Gray scale, 1 slope total over the course
                 * of 256 contours.
                 */
                red_value = 1;
                green_value = 1;
                blue_value = 1;
                rstep = 1;
                gstep = 1;
                bstep = 1;
                break;
            case 2:
                /*
                 * 2: Gray scale, 4 slopes over 256 contours.
                 */
                red_value = 4;
                green_value = 4;
                blue_value = 4;
                rstep = 4;
                gstep = 4;
                bstep = 4;
                break;
            case 3: /*
                 * 3: Gold. Properly used, this can
                 * yield a picture which looks as if it glows.
                 * This is why we use floating point values
                 * for color construction...
                 */
                red_value = 16;
                green_value = 16;
                blue_value = 16;
                rstep = 0.935;
                gstep = 0.935;
                bstep = 0.5;
                break;
            default:
                fprintf(stderr, "Unknown color set %d\n", c);
                exit(1);

            }
        } else if (strcmp(argv[i],"/s") == 0) {
            /*
             * /s filename
             * store the results into that filename as a
             * ppm-format file. This can get BIG.
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /s\n");
                exit(1);
            }
            save_file = fopen(argv[++i], "w");
            if (save_file == NULL) {
                fprintf(stderr,
                    "Couldn't open file <%s> for output\n",
                    argv[i]);
                exit(1);
            }
        } else if (strcmp(argv[i],"/r") == 0) {
            /*
             * /r CENTER of the R axis
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /r\n");
                exit(1);
            }
            r_center = atof(argv[++i]);
        } else if (strcmp(argv[i],"/w") == 0) {
            /*
             * /w width
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /w\n");
                exit(1);
            }
            width = atof(argv[++i]);
        } else if (strcmp(argv[i],"/i") == 0) {
            /*
             * /i CENTER of the I axis
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /i\n");
                exit(1);
            }
            /*
             * In order to have the same coordinate system as
             * "MandelVroom" this value has to be negated.
             * This allows MandelVroom to be used as a sort of
             * spotter-scope for this program.
             *
             * ...Which is a great theory, but for reasons I don't
             * understand, there is some sort of shift between
             * "MandelVroom" coordinates and those used here. It is
             * entirely possible that this is due to a systemic error
             * I've induced in my optimized computation function.
             * Anyway, "MandelVroom" gets you near, and you can
             * refine it from there using /q.
             */
            i_center = -atof(argv[++i]);
        } else if (strcmp(argv[i],"/q") == 0) {
            /*
             * /q - use only a quarter of the screen (per direction)
             * which means a 16th of the area.
             */
            quarter_flag = TRUE;
        } else if (strcmp(argv[i],"/h") == 0) {
            /*
             * /h - use only a half of the screen (per direction)
             * which means a quarter of the area.
             */
            half_flag = TRUE;
        } else if (strcmp(argv[i],"/d") == 0) {
            /*
             * /d - change the step-size
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /d\n");
                exit(1);
            }
            step_size = atol(argv[++i]);
            if (step_size > 126) {
                fprintf(stderr, "/d divider cannot exceed 126\n");
                exit(1);
            }
        } else if (strcmp(argv[i],"/t") == 0) {
            /*
             * /t - Raise the ceiling.
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /t\n");
                exit(1);
            }
            max_orbits = atol(argv[++i]);
            if (max_orbits > 32767) {
                fprintf(stderr, "/t cannot exceed 32767\n");
                exit(1);
            }
        } else if (strcmp(argv[i],"/m") == 0) {
            /*
             * /m threshold
             * Traditionally the threshold is 4.0. But what if
             * it isn't?
             */
            if (argc < i+1) {
                fprintf(stderr, "Not enough parameters for /m\n");
                exit(1);
            }
            threshold = atof(argv[++i]);
        } else {
            fprintf(stderr, "Unknown parameter <%s>\n", argv[i]);
            exit(1);
        }
    }

    /*
     * SCDB: Compute the height to maintain an aspect ratio of 1.2 to 1.
     */
    height = width * 0.83;

    /*
     * SCDB: Calculate the left and top coordinates based on center, width
     * and height values.
     */
    left_coord = r_center-width/2;
    top_coord = i_center+height/2;

  if ((RetinaBase = (struct _xy_RetinaBase *)
       OpenLibrary("retina.library", 2L)) != NULL) {
    if ((IntuitionBase = (struct IntuitionBase *)
     OpenLibrary("intuition.library", 37L)) != NULL) {
      if ((rScreen = Retina_OpenScreen(RSCR_STDWIDTH, RSCR_STDHEIGHT,
                       MID_DEFAULT_08,
                       RSFF_AUTOADJUST, 0L))
      != NULL) {
    printf("Screen Opened successfully: %d by %d\n",
           rScreen->_rs_Width, rScreen->_rs_Height);

    /*
     * SCDB: Well, ONE of the following flags has to have precedence. If the
     * user is silly enough to specify both, we'll use the smaller one.
     */
    if (quarter_flag) {
        screen_width  = rScreen->_rs_Width/4;
        h_offset      = (rScreen->_rs_Width*3)/8; /* to center it */
        screen_height = rScreen->_rs_Height/4;
        v_offset      = (rScreen->_rs_Height*3)/8; /* to center it */
    } else if (half_flag) {
        screen_width  =  rScreen->_rs_Width/2;
        h_offset      = rScreen->_rs_Width/4; /* to center it */
        screen_height =  rScreen->_rs_Height/2;
        v_offset      = rScreen->_rs_Height/4; /* etc. */
    } else {
        screen_width = rScreen->_rs_Width;
        screen_height = rScreen->_rs_Height;
    }


    /* set up screen colors; color 0 must always be black */
    Retina_SetPalette(rScreen, 0, 0, 0, 0);
    palette[0] = 0;
    for (i = 1; i < 256; i++) {
      palette[i] = ((long)red_value<<16)+((long)green_value<<8)+(long)blue_value;
      Retina_SetPalette(rScreen, i, red_value, green_value, blue_value);
      red_value += rstep;
      green_value += gstep;
      blue_value += bstep;
    }

    Retina_SetDrMd(rScreen, RDM_JAM1);

    /*
     * SCDB: I deduce that 'dr' is the the numeric increment in
     * calculation-space represented by one pixel to the right, and
     * 'di' is the numeric increment in calculation-sapce of one pixel down.
     */
    dr = (width)/((double) screen_width);
    di = -(height)/((double) screen_height);

    /* OK, folks, enough initializing.  Let's call the
       Mbrot routine. */

    /*
     * SCDB:
     * The screen is created in four quadrants. Create the upper left
     * quadrant first.
     * The parameters are:
     *    "x" is the pixel position horizontally.
     *    "y" is the pixel position vertically.
     *    "r" is the value in calculation-space horizontally.
     *        (0 to the left)
     *    "i" is the value in calculation-space vertically.
     *        (0 below, but we inverted the parameter so 0 is
     *         effectively above)
     *
     * The original code used incremental changes on the structures
     * to set up these values, but in the interest of clarity
     * I'm going to do them completely for each quadrant. The
     * increase in execution time is miniscule.
     *
     * A long description of a change I made: I don't think that the
     * original code handled this correctly. Both axes were handled the same
     * way, so I'll describe only one of them for simplicity's sake. The
     * original code described the first quadrant as 0-(width>>1) and the
     * second as (width>>1)-width; which bothered me because
     * different screen pixels were used for the last of the first
     * quadrant and the first of the other quadrant, but they were computed
     * using the same value. I modified this so that the first quadrant goes
     * 0-((width>>1)-1) but I think I missed something in the logic.
     * As I've observed it writing the screen, it
     * sometimes redundantly writes over pixels it has previously calculated
     * BUT WITH DIFFERENT COLORS. There's no way that can be correct.
     * The result is still esthetically pleasing and shows no important
     * gaps or asymmetries, so I haven't pursued it. I don't understand the
     * code which recursively divides squares below this point, and the
     * problem lies in there somewhere.
     * This is another invitation for change to anyone interested.
     *
     * First describe the pixel set of the quadrant
     */
    start.x = 0;
    end.x = (screen_width >> 1)-1;
    start.y = 0;
    end.y = (screen_height >> 1)-1;

    /*
     * Then derive the numeric set for the quadrant.
     * This is the same for every quadrant. I've described it this way
     * because it has negligible effect on execution time, but makes the
     * code clearer.
     */
    start.r = left_coord + dr*(double)start.x;
    end.r = left_coord + dr*(double)end.x;
    start.i = top_coord + di*(double)start.y;
    end.i = top_coord + di*(double)end.y;


#ifndef NOREPEATCALC
    start.orbits = -2;
    end.orbits = -2;
#endif
    Compute_Block(start, end);


    /*
     * Create the lower left quadrant.
     */
    start.x = 0;
    end.x = (screen_width >> 1)-1;
    start.y = (screen_height>>1);
    end.y = screen_height-1;

    start.r = left_coord + dr*(double)start.x;
    end.r = left_coord + dr*(double)end.x;
    start.i = top_coord + di*(double)start.y;
    end.i = top_coord + di*(double)end.y;

    Compute_Block(start, end);

    /*
     * Create the lower right quadrant.
     */
    start.x = (screen_width>>1);
    end.x = screen_width-1;
    start.y = (screen_height>>1);
    end.y = screen_height-1;

    start.r = left_coord + dr*(double)start.x;
    end.r = left_coord + dr*(double)end.x;
    start.i = top_coord + di*(double)start.y;
    end.i = top_coord + di*(double)end.y;

    Compute_Block(start, end);

    /*
     * Create the upper right quadrant;
     */
    start.x = (screen_width>>1);
    end.x = screen_width-1;
    start.y = 0;
    end.y = (screen_height >> 1)-1;

    start.r = left_coord + dr*(double)start.x;
    end.r = left_coord + dr*(double)end.x;
    start.i = top_coord + di*(double)start.y;
    end.i = top_coord + di*(double)end.y;

    Compute_Block(start, end);


    if (save_file != NULL) {
        /*
         * Then the user wants the results saved to a PPM file.
         * The file has already been opened, so the only thing we
         * have to do is write it out.
         *
         * The first part of a PPM is "P6" to indicate that
         * it IS a PPM, the x and y sizes and "255" to
         * indicate full brightness level. This is then followed
         * by byte triples (R,G,B) scanning x within y.
         * It took me a while to realize that the placement of spaces
         * and newlines here is critical for some programs which are
         * picky.
         */
        fprintf(save_file, "P6\n%d %d\n255\n",
            screen_width, screen_height);
        i = 0;
        for (y_scan=0; y_scan<screen_height; y_scan++) {
            for (x_scan=0; x_scan<screen_width; x_scan++) {

            /*
             * The documentation says that Retina_ReadPixel returns
             * the RGB value. Actually for 8 bit screens it returns
             * the palette index for that pixel.
             * We therefore have to translate
             * the palette index into RGB based on the values we
             * saved when we set the palette up in the first place.
             */
            color_value = Retina_ReadPixel(rScreen,
                x_scan+h_offset, y_scan+v_offset);
            color_value = palette[color_value];

            out_buffer[i++] = (unsigned char)
                ((color_value >> 16) & 0xFFL);
            out_buffer[i++] = (unsigned char)
                ((color_value >> 8) & 0xFFL);
            out_buffer[i++] = (unsigned char)(color_value & 0xFFL);

            /*
             * The output has been accumulated in a big buffer
             * because programs run faster if they write out a few
             * big blocks than if they write out a lot of small
             * ones, which is what would have happened if we'd
             * been using "fwrite" on the color values directly.
             */
            if (i>=BUFF_SIZE) {
                fwrite(out_buffer, i, 1, save_file);
                i = 0;
            }
            }
        }

        /*
         * We've put all values into the output buffer, and most of
         * it has probably been written out. If there is any residue,
         * write it out, too.
         */
        if (i != 0)
            fwrite(out_buffer, i, 1, save_file);
        fclose(save_file);

    } else {

        /*
         * If no output file was specified, leave the picture on the
         * screen until the user indicates it should go away.
         * For people (like me) who use retina workbench emulation,
         * the following printf is useless - but
         * there might be someone out there who is keeping
         * their workbench on something else. The rest of us are just
         * going to have to remember to hit CR afterwards.
         *
         * In theory, any keypress should work. In practice, "getchar"
         * (at least for Aztec) stacks up all input until CR is pressed
         * and only presents it to the program at that point.
         */
        printf("Press -- ENTER -- to quit.\n");
        getchar();
    }

    Retina_CloseScreen(rScreen);
      }
      else {
    fprintf(stderr, "Unable to open Retina screen.\n");
      }
    }
    else {
      fprintf(stderr, "Can't open intuition.library V37+\n");
    }
  }
  else {
    fprintf(stderr, "Can not open retina.library!\n");
  }
  exit(0);
}

/*
 * SCDB: I do not understand the following code and have made no attempt
 * to debug or optimize it. Heare bee dragones.
 */

void Compute_Block(struct BrotPoint uleft, struct BrotPoint lright)
{
  unsigned int i;
  double j;
  short int q, z;

  if (lright.x - uleft.x < 2 && lright.y - uleft.y < 2) {
#ifdef NOREPEATCALC
    q = Test_Point(uleft.r, uleft.i)/step_size;
    if (q > 0) Retina_SetPixel(rScreen, uleft.x+h_offset, uleft.y+v_offset, q);
#else
    if (uleft.orbits == -2) {
      uleft.orbits = Test_Point(uleft.r, uleft.i)/step_size;

    }
    if (uleft.orbits > 0) Retina_SetPixel(rScreen, uleft.x+h_offset,
                uleft.y+v_offset, uleft.orbits);
#endif
    q = Test_Point(uleft.r, lright.i)/step_size;

    if (q > 0) Retina_SetPixel(rScreen, uleft.x+h_offset, lright.y+v_offset, q);
#ifdef NOREPEATCALC
    q = Test_Point(lright.r, lright.i)/step_size;
    if (q > 0) Retina_SetPixel(rScreen, lright.x+h_offset, lright.y+v_offset, q);
#else
    if (lright.orbits == -2) {
      lright.orbits = Test_Point(lright.r, lright.i)/step_size;
    }

    if (lright.orbits > 0) Retina_SetPixel(rScreen, lright.x+h_offset,
                lright.y+v_offset, lright.orbits);
#endif
    q = Test_Point(lright.r, uleft.i)/step_size;

    if (q > 0) Retina_SetPixel(rScreen, lright.x+h_offset, uleft.y+v_offset, q);
    return;
  }

  /* draw a box -- if the borders are all the same, then fill it that
     color */
  q = Test_Point(uleft.r, uleft.i)/step_size;

  uleft.orbits = q;
  if (q > 0)
    Retina_SetPixel(rScreen, uleft.x+h_offset, uleft.y+v_offset, q);
  z = Test_Point(lright.r, lright.i)/step_size;

  lright.orbits = z;
  if (z > 0)
    Retina_SetPixel(rScreen, lright.x+h_offset, lright.y+v_offset, z);
  if (z != q) {
    goto recurse;
  }

  /* top & bottom */
  for (i = uleft.x, j = uleft.r; i <= lright.x; i++, j += dr) {
    z = Test_Point(j, uleft.i)/step_size;

    if (z != q) {
      goto recurse;
    }
    else {
      if (z > 0)
    Retina_SetPixel(rScreen, i+h_offset, uleft.y+v_offset, z);
    }
    z = Test_Point(j, lright.i)/step_size;

    if (z != q) {
      goto recurse;
    }
    else {
      if (z > 0)
    Retina_SetPixel(rScreen, i+h_offset, lright.y+v_offset, z);
    }
  }
  /* left and right */
  for (i = uleft.y, j = uleft.i; i <= lright.y; i++, j += di) {
    z = Test_Point(uleft.r, j)/step_size;

    if (z != q) {
      goto recurse;
    }
    else {
      if (z > 0)
    Retina_SetPixel(rScreen, uleft.x+h_offset, i+v_offset, z);
    }
    z = Test_Point(lright.r, j)/step_size;

    if (z != q) {
      goto recurse;
    }
    else {
      if (z > 0)
    Retina_SetPixel(rScreen, lright.x+h_offset, i+v_offset, z);
    }
  }

  /* if we passed through all that, we draw a filled rectangle */
  if (q > 0) {
    Retina_SetAPen(rScreen, q);

    my_RectFill(uleft.x, uleft.y, lright.x, lright.y);
  }
  return;

 recurse:
  {
    struct BrotPoint top, bot, left, right, middle;

    top.x = bot.x = middle.x = (uleft.x + lright.x) >> 1;
    left.y = right.y = middle.y = (uleft.y + lright.y) >> 1;
    top.r = bot.r = middle.r = (uleft.r + lright.r)/2.0;
    left.i = right.i = middle.i = (uleft.i + lright.i)/2.0;
    top.y = uleft.y;
    top.i = uleft.i;
    bot.y = lright.y;
    bot.i = lright.i;
    left.x = uleft.x;
    left.r = uleft.r;
    right.x = lright.x;
    right.r = lright.r;
#ifndef NOREPEATCALC
    middle.orbits = top.orbits = bot.orbits = left.orbits =
      right.orbits = -2;
#endif

    Compute_Block(uleft, middle);
#ifdef NOREPEATCALC
    left.y += 1;
    left.i += di;
#endif
    Compute_Block(left, bot);
#ifdef NOREPEATCALC
    middle.x += 1;
    middle.y += 1;
    middle.r += dr;
    middle.i += di;
    top.x += 1;
    top.r += dr;
#endif
    Compute_Block(middle, lright);
    Compute_Block(top, right);
  }
}

/*
 * SCDB: The routine which follows used to be controlled by an ifdef. I've
 * made it unconditional, since it works very well.
 * "my" is the original author, not me. He did good.
 */
void my_RectFill(unsigned int x1, unsigned int y1,
                    unsigned int x2, unsigned int y2)
{
  short int tx1, tx2, v;

  /* a little magic because Retina_RectFill() ignores the low-order
     3 bits of the x-coords for 256-color screens */
  tx1 = (x1 >> 3) + 1; tx2 = (x2 >> 3) - 1;
  tx1 <<= 3;
  tx2 <<= 3;

  if (tx2 > tx1) {
    Retina_RectFill(rScreen, tx1+h_offset, y1+v_offset, tx2+h_offset, y2+v_offset);
    if (x2 - tx2 > 1) {
      for (v = tx2; v <= x2; v++)
    Retina_Line(rScreen, v+h_offset, y1+v_offset, v+h_offset, y2+v_offset);
    }
    if (tx1 - x1 > 1) {
      for (v = x1; v <= tx1; v++)
    Retina_Line(rScreen, v+h_offset, y1+v_offset, v+h_offset, y2+v_offset);
    }
  }
  else {
    /* fake it with line drawing */
    for (v = y1; v <= y2; v++) {
      Retina_Line(rScreen, x1+h_offset, v+v_offset, x2+h_offset, v+v_offset);
    }
  }
}

#ifdef UNUSED
/*
 * SCDB: Sorry folks, a bit of a debugging hook. I never got the hang of the
 * Aztec debugger, powerful though it is, since I've never really had the need.
 * My traditional way of debugging is to sprinkle printf's all through the code
 * to figure what the hell is happening. These are mechanisms which can
 * be called from the FPU routine to see what's happening there.
 * Unless further debugging is required there by someone using as primitive
 * of techniques as I used, they are unnecessary.
 */
int    trace_buff[32] = {-1,-1,-1,-1,-1,-1,-1,-1
              -1,-1,-1,-1,-1,-1,-1,-1,
              -1,-1,-1,-1,-1,-1,-1,-1,
              -1,-1,-1,-1,-1,-1,-1,-1};

int    trace_index=0;

void trace(val)
int val;
{
    trace_buff[trace_index] = val;
    if (++trace_index >= 32)
        trace_index = 0;
}

void trace_dump()
{
    int    i;

    for (i=trace_index; i<32; i++)
        printf("%d: %d\n", i, trace_buff[i]);

    for (i=0; i<trace_index; i++)
        printf("%d: %d\n", i, trace_buff[i]);
}
#endif