mcgrew@dartagnan.rutgers.edu (Charles Mcgrew) (09/13/89)
Submitted-by: P{r Emanuelsson <pell@isy.liu.se>
Posting-number: Volume 1, Issue 60
Archive-name: FastMand
I have written a program to implement the FastM algorithm as described
in the book "The Science of Fractal Images" edited by Peitgen and Saupe.
The FastM algorithm is good for fast calculation of a (rough) outline of
the Mandelbrot set and thus suitable for interactive manipulations. This
makes it possible to iterate e.g. only every 10:th row and column of the image,
reducing the time by a factor of 100 compared to the standard Mandelbrot
algorithm. It's also fun to watch, since it's recursive...
The shar file to "FastMtool" follows.
Have fun! But beware - it's addictive!
/Pell
--
#!/bin/sh
# shar: Shell Archiver (v1.22)
#
# Run the following text with /bin/sh to create:
# README
# Makefile
# FastMtool.1
# FastMtool.c
#
if test -f README; then echo "File README exists"; else
sed 's/^X//' << 'SHAR_EOF' > README &&
X This is a program for interactive Mandelbrot exploration for Sun
Xworkstations using SunView.
X
X What's new and exciting with this program (IMHO :-) ) is that it implements
Xa new algorithm. This is good for fast calculation of a (rough) outline of
Xthe Mandelbrot set and thus suitable for interactive manipulations. This
Xmakes it possible to iterate e.g. only every 10:th row and column of the image,
Xreducing the time by a factor of 100 compared to the standard Mandelbrot
Xalgorithm. It's also fun to watch, since it's recursive...
X
X It has been tested on the following configurations:
X
XSun-2/120 under SunOS 3.5, Sun-3/50, 3/75, 3/110 with FPA, 3/280, 4/110
Xand 4/280, all using SunOS 4.0(.1). Hopefully it should work on others too.
X
X So how does it work? Well, using some intricate mathematics you can
Xarrive at an equation that estimates the distance from a point outside the
XMandelbrot set to the nearest point IN the set. Using this knowledge, it's
Xpossible to exclude all points lying inside a disk with the radius equal to
Xthis distance. You know they can't be part of the set and mark them as
Xuninteresting. Then you repeat the calculation with a couple of points
Xpicked from the edge of this disk. In this fashion, you will get an outline
Xof the set in a very short time.
X
X Time for the bad news: if you want an image with the same resolution
Xas one generated with the standard algorithm, this algorithm will be just
Xas slow. This is because this algorithm quickly excludes points that are
Xnot part of the Mandelbrot set. Unfortunately, the calculations for these
Xpoints will usually be fast since they rapidly diverges to infinity. So the
Xpoints left that need to be iterated are all very time-consuming. Furthermore,
Xthis algorithm cannot provide those good-looking color images you have all
Xseen. It only provides information whether a point is in the set or not.
XOf course, on a B/W screen this will usually suffice. Using this program, you
Xcan quickly find interesting areas and then run off to your cray and feed the
Xcoordinates to a standard Mandelbrot algorithm.
X
X I certainly didn't invent this myself. I urge you to get the book
X"The Science of Fractal Images" edited by Peitgen and Saupe. Springer Verlag,
XISBN 0-387-96608-0. It's filled with good stuff like this.
X
X A possible enhancement would be for the program to draw disks INSIDE
Xthe Mandelbrot set too. Unfortunately, this problem is quite non-deterministic
Xand does not have a simple solution. Besides, the set, at larger magnifications,
Xhas a very intricate structure and it would probably not be possible to
Xdraw any large disks inside it.
X
X I refer you to the man-page for more instructions about running the stuff.
XTo get going, just hit the DRAW button as soon as it starts.
X
X You will probably see a lot of inefficient coding, e.g. using array
Xsubscripting instead of pointers etc. Let me just say that profiling shows
Xthat without the Sun-specific stuff, it can spend as much as 98% of the
Xtime in the MSetDist routine (most of which are floating-point
Xcalculations)... That is why I have kept the inefficient, but more
Xreadable, code. Using GCC, you cannot achive much better speed hand-coding
XMSetDist in assembler, but you are welcome to try...
X
X Finally, don't flame me for the coding. I just wanted to try this
Xalgorithm, and whipped together a demonstration program using parts from
Xother programs in a very short time. Specifically: don't run it through lint!
X(I haven't :-). I don't have the time to work on this anymore but I thought
Xthis was interesting enough to warrant posting.
X
X One thing I would really love is an X version of this (hint, hint!).
XPoor me only knows the workings of SunView... :-(
X
XHave fun!
X
X /Pell
X--
X"Don't think; let the machine do it for you!"
X -- E. C. Berkeley
XDept. of Electrical Engineering ====>>> pell@isy.liu.se
XUniversity of Linkoping, Sweden ...!uunet!enea!isy.liu.se!pell
SHAR_EOF
chmod 0644 README || echo "restore of README fails"
fi
if test -f Makefile; then echo "File Makefile exists"; else
sed 's/^X//' << 'SHAR_EOF' > Makefile &&
XLIBS = -lm -lsuntool -lsunwindow -lpixrect
X
X# Choose the compile line that suits you best.
X# Note that cc -O4 gives 25% faster code that gcc 1.30. I'm not sure why.
X# If you have a later version you may want to try gcc anyway, perhaps
X# with fancier optimizations...
X
XFastMtool: FastMtool.c
X# GCC, Sun-3 with 68881 math chip.
X# gcc -O -traditional -m68881 -o FastMtool FastMtool.c $(LIBS)
X# Sun-3 with SunOS 4, 68881 math chip:
X cc -O4 -f68881 -o FastMtool FastMtool.c $(LIBS)
X# Sun-4 with SunOS 4:
X# cc -O4 -o FastMtool FastMtool.c $(LIBS)
X# Other suns, configure for maximum speed yourself:
X# cc -O -o FastMtool FastMtool.c $(LIBS)
SHAR_EOF
chmod 0644 Makefile || echo "restore of Makefile fails"
fi
if test -f FastMtool.1; then echo "File FastMtool.1 exists"; else
sed 's/^X//' << 'SHAR_EOF' > FastMtool.1 &&
X.TH FastMtool 1 "1989 February 25"
X.UC 4
X.SH NAME
XFastMtool \- Explore the Mandelbrot set under SunView(1)
X.SH SYNOPSIS
X.B FastMtool
X.SH DESCRIPTION
X.I FastMtool
Xuses a new algorithm, presented in the book
X.B
X"The Science of Fractal Images"
Xpublished by Springer Verlag. This algorithm is well suited to interactive
Xexplorations, providing a rough outline of the Mandelbrot set very quickly.
XUsing the mouse, you can then select and enlarge interesting regions.
X.PP
X.B Controls:
X.PP
XThere are three sliders. Their functions are:
X.IP Recur
XAffects the recursive part of the algorithm.
X.I Recur
Xis the minimum distance to the Mandelbrot set at which the recursion
Xshould stop. I.e. if
X.I Recur
Xis equal to 10 then the recursive part of the algorithm will never try
Xto get closer than 10 pixels to the set. 5 is probably a good value
Xto start with.
X.IP Increment
XAffects the iterative part of the algorithm. The algorithm starts by scanning
Xthe image line after line, pixel for pixel. When it hits a pixel not part
Xof the set, it invokes the recursive algorithm. Then it continues scanning.
X
XIf
X.I Increment
Xis set e.g. to 5, then only every 5:th line and pixel will be scanned,
Xreducing the time 25-fold.
X.IP Maxiter
XThis determines how many iterations are required before a pixel is considered
Xpart of the set. One can view this parameter as the "focusing". This parameter
Xneeds to be increased as the magnification increases. If it is set too low
Xyou will not get any details in the image; too many pixels will erroneously
Xbe considered part of the set. Increasing this also increases the time.
XExperiment!
X.PP
X.B Buttons:
X.PP
XThere are four buttons:
X.IP DRAW
XDraws a new image. If you haven't selected new coordinates the old image
Xis redrawn, perhaps with new parameters.
X.IP UP
XReturns you to your previous image, before the last
X.I DRAW
Xwith new coordinates.
X.IP STOP
XEmergency stop. It's only possible to stop the iterative refinement. The
Xrecursive refinement is usually fast anyway.
X.IP Quit
XQuits the program (without confirm on SunOS older than 4.0).
X.PP
X.B
XMouse usage:
X.PP
XTo mark a new area for exploration, press the
X.I left
Xmouse button and move the mouse, keeping the button down. When you release
Xthe button a new area has been marked. Now you can select the
X.I DRAW
Xbutton to magnify this area, perhaps changing some of the sliders first.
X
XYou will see the coordinates in the panel change when you move the mouse.
X
XIf you hit the
X.I middle
Xbutton, your selected area will be cancelled and the coordinates restored.
X.SH RESTRICTIONS
XYou cannot save the images to file.
X.SH SEE ALSO
Xsunview(1)
X.SH AUTHORS
X.PP
XP. Emanuelsson, pell@isy.liu.se.
X.br
X.I
X"The Science of Fractal Images"
X\(em Peitgen & Saupe (eds.)
SHAR_EOF
chmod 0644 FastMtool.1 || echo "restore of FastMtool.1 fails"
fi
if test -f FastMtool.c; then echo "File FastMtool.c exists"; else
sed 's/^X//' << 'SHAR_EOF' > FastMtool.c &&
X/*
X * FastMtool - Perform interactive exploring of the Mandelbrot set on
X * Sun workstations.
X * Copyright (c) 1989 P{r Emanuelsson, pell@isy.liu.se. All Rights Reserved.
X * This program is free software; you can redistribute it and/or modify
X * it under the terms of the GNU General Public License as published by
X * the Free Software Foundation; either version 1, or any later version.
X * This program is distributed in the hope that it will be useful,
X * but WITHOUT ANY WARRANTY; without even the implied warranty of
X * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
X * GNU General Public License for more details.
X * You can receive a copy of the GNU General Public License from the
X * Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
X * See the README for more information.
X */
X
X#include <suntool/sunview.h>
X#include <suntool/panel.h>
X#include <suntool/canvas.h>
X
X/* Change event_action to event_id(event) for pre 4.0 */
X#ifndef event_action
X#define event_action event_id
X#define oldSunOS
X#endif
X
X#ifndef oldSunOS
X#include <suntool/alert.h> /* Only for SunOS 4 */
X#include <alloca.h> /* Cannot find this on 3.5 */
X#endif
X
X#include <math.h>
X
Xtypedef unsigned char Bool;
X
X#define Stacksize 100 /* Dynamic allocation would be cleaner */
X
Xstruct stack {
X Pixrect *pr;
X double xmin, xmax, ymin, ymax;
X} Stack [Stacksize];
X
Xstatic int sp = 0; /* Stackpointer */
X
X#define MAX(a,b) ((a) > (b) ? (a) : (b))
X#define MIN(a,b) ((a) < (b) ? (a) : (b))
X#define sgn(x) ((x) >= 0 ? 1 : -1)
X
X#define SQRT_2 1.41421356237
X#define Nmax 800 /* Maximum image size */
X
Xstatic Bool MSet[Nmax][Nmax];
Xstatic double xmin = -2.0, xmax=0.5, ymin = -1.25, ymax=1.25,
X side=2.5, delta, recur;
Xstatic int maxiter=20, incr=5, rec=5;
Xstatic int start_x, start_y, new_x, new_y; /* For region select */
Xstatic Bool selected = FALSE, draw_new_image=FALSE, abort=FALSE;
X
X#define BUTTONFONT "/usr/lib/fonts/fixedwidthfonts/gallant.r.19"
Xstatic char *default_font =
X "DEFAULT_FONT=/usr/lib/fonts/fixedwidthfonts/screen.r.13";
X
Xstatic Frame fr;
Xstatic Canvas cv;
Xstatic Pixwin *pw;
Xstatic Pixfont *bf, *pf;
Xstatic Panel pn;
Xstatic Panel_item recur_p, incr_p, maxiter_p, stop_p;
X
X#ifdef notdef /* Future enhancement */
Xstatic int N=512; /* Current image size */
X#else
X#define N 512
X#endif
X
Xstatic char Nmax_str[10];
X
Xextern int panel_pw_text(), panel_fonthome(); /* Undocumented, but nice */
X
Xdouble MSetDist();
Xvoid quit(), stop(), up(), draw(), sel_region();
X
Xvoid main()
X{
X (void) putenv(default_font);
X bf = pf_open(BUTTONFONT);
X pf = pf_default();
X
X/* This code is full of strange magic numbers. I apologize for that!!! */
X
X fr = window_create(0, FRAME,
X WIN_X, 400,
X WIN_Y, 150,
X WIN_SHOW, TRUE,
X FRAME_LABEL, "FastMtool - Mandelbrot exploring",
X 0);
X pn = window_create(fr, PANEL, 0);
X
X/* The brain-damaged SunView uses the PANEL_MAX_VALUE to determine the
X position of the slider, instead of e.g. always allocate 10 positions
X for this and right-justify or something. This means that I have to
X use delicate (ugly) pixel-positioning to get the sliders to line up. */
X
X recur_p = panel_create_item(pn, PANEL_SLIDER,
X PANEL_LABEL_X, 5,
X PANEL_LABEL_Y, 5,
X PANEL_VALUE_X, 110,
X PANEL_VALUE_Y, 5,
X PANEL_LABEL_STRING, "Recur:",
X PANEL_VALUE, rec,
X PANEL_MIN_VALUE, 1,
X PANEL_MAX_VALUE, 50,
X 0);
X incr_p = panel_create_item(pn, PANEL_SLIDER,
X PANEL_LABEL_X, 5,
X PANEL_LABEL_Y, 25,
X PANEL_VALUE_X, 110,
X PANEL_VALUE_Y, 25,
X PANEL_LABEL_STRING, "Increment:",
X PANEL_VALUE, incr,
X PANEL_MIN_VALUE, 1,
X PANEL_MAX_VALUE, 10,
X 0);
X maxiter_p = panel_create_item(pn, PANEL_SLIDER,
X PANEL_LABEL_X, 5,
X PANEL_LABEL_Y, 45,
X PANEL_VALUE_X, 78,
X PANEL_VALUE_Y, 45,
X PANEL_LABEL_STRING, "Maxiter:",
X PANEL_VALUE, maxiter,
X PANEL_MIN_VALUE, 10,
X PANEL_MAX_VALUE, 1000,
X 0);
X panel_create_item(pn, PANEL_MESSAGE,
X PANEL_ITEM_X, 5,
X PANEL_ITEM_Y, 100,
X PANEL_LABEL_STRING, "Xmin:",
X PANEL_LABEL_BOLD, TRUE,
X 0);
X panel_create_item(pn, PANEL_MESSAGE,
X PANEL_ITEM_X, 200,
X PANEL_ITEM_Y, 100,
X PANEL_LABEL_STRING, "Xmax:",
X PANEL_LABEL_BOLD, TRUE,
X 0);
X panel_create_item(pn, PANEL_MESSAGE,
X PANEL_ITEM_X, 5,
X PANEL_ITEM_Y, 120,
X PANEL_LABEL_STRING, "Ymin:",
X PANEL_LABEL_BOLD, TRUE,
X 0);
X panel_create_item(pn, PANEL_MESSAGE,
X PANEL_ITEM_X, 200,
X PANEL_ITEM_Y, 120,
X PANEL_LABEL_STRING, "Ymax:",
X PANEL_LABEL_BOLD, TRUE,
X 0);
X
X#ifdef notdef /* Possible future enhancement... */
X sprintf(Nmax_str, "%d", Nmax);
X panel_create_item(pn, PANEL_CYCLE,
X PANEL_ITEM_X, 350,
X PANEL_ITEM_Y, 5,
X PANEL_LABEL_STRING, "Image size:",
X PANEL_LABEL_BOLD, TRUE,
X PANEL_CHOICE_STRINGS, Nmax_str, "512", "256", 0,
X PANEL_VALUE, 1,
X 0);
X#endif
X
X panel_create_item(pn, PANEL_BUTTON,
X PANEL_ITEM_X, 360,
X PANEL_ITEM_Y, 30,
X PANEL_LABEL_IMAGE, panel_button_image(pn, "DRAW", 0, bf),
X PANEL_NOTIFY_PROC, draw,
X 0);
X panel_create_item(pn, PANEL_BUTTON,
X PANEL_ITEM_X, 360,
X PANEL_ITEM_Y, 70,
X PANEL_LABEL_IMAGE, panel_button_image(pn, " UP ", 0, bf),
X PANEL_NOTIFY_PROC, up,
X 0);
X panel_create_item(pn, PANEL_BUTTON,
X PANEL_ITEM_X, 360,
X PANEL_ITEM_Y, 110,
X PANEL_LABEL_IMAGE, panel_button_image(pn, "STOP", 0, bf),
X PANEL_NOTIFY_PROC, stop,
X 0);
X panel_create_item(pn, PANEL_BUTTON,
X PANEL_ITEM_X, 450,
X PANEL_ITEM_Y, 72,
X PANEL_LABEL_IMAGE, panel_button_image(pn, "Quit", 0, pf),
X PANEL_NOTIFY_PROC, quit,
X 0);
X window_fit_height(pn);
X cv = window_create(fr, CANVAS,
X WIN_WIDTH, N,
X WIN_HEIGHT, N,
X WIN_CONSUME_PICK_EVENT, LOC_DRAG,
X WIN_EVENT_PROC, sel_region,
X 0);
X
X window_fit(fr);
X notify_interpose_destroy_func(fr, quit);
X
X pw = canvas_pixwin(cv);
X notify_dispatch(); /* To make the next put_coords work */
X put_coords(0, N-1, N-1, 0);
X
X bzero((char *) MSet, sizeof(MSet));
X pw_writebackground(pw, 0, 0, N, N, PIX_SRC | PIX_COLOR(1));
X
X /* Cannot use window_main_loop() because the notifier can't dispatch
X events inside a PANEL_NOTIFY_PROC. */
X
X while (1) {
X notify_dispatch();
X if (draw_new_image) {
X panel_pw_text(pn, 347, 5 + panel_fonthome(bf), PIX_SRC, bf, "Drawing!");
X compute();
X panel_pw_text(pn, 347, 5 + panel_fonthome(bf), PIX_SRC, bf, " ");
X draw_new_image = FALSE;
X }
X usleep(50000);
X }
X}
X
Xput_coords(ixmin, ixmax, iymin, iymax)
X int ixmin, ixmax, iymin, iymax;
X{
X char str[20];
X
X sprintf(str, "%10.7lf", xmin + side * ixmin/(N-1));
X panel_pw_text(pn, 50, 100 + panel_fonthome(pf), PIX_SRC, pf, str);
X sprintf(str, "%10.7lf", xmin + side * ixmax/(N-1));
X panel_pw_text(pn, 245, 100 + panel_fonthome(pf), PIX_SRC, pf, str);
X sprintf(str, "%10.7lf", ymin + side * (N-1-iymin)/(N-1));
X panel_pw_text(pn, 50, 120 + panel_fonthome(pf), PIX_SRC, pf, str);
X sprintf(str, "%10.7lf", ymin + side * (N-1-iymax)/(N-1));
X panel_pw_text(pn, 245, 120 + panel_fonthome(pf), PIX_SRC, pf, str);
X}
X
Xvoid sel_region(canvas, event, arg)
X Canvas canvas;
X Event *event;
X char *arg;
X{
X static Bool mouseing = FALSE;
X register int maxdist, tmpx, tmpy;
X
X#define mkbox(a,b,c,d) \
X pw_vector(pw, a, b, c, b, PIX_SRC^PIX_DST, 1);\
X pw_vector(pw, a, b, a, d, PIX_SRC^PIX_DST, 1);\
X pw_vector(pw, c, b, c, d, PIX_SRC^PIX_DST, 1);\
X pw_vector(pw, a, d, c, d, PIX_SRC^PIX_DST, 1)
X
X switch (event_action(event)) {
X case MS_LEFT:
X if (event_is_down(event)) {
X if (selected) { /* Remove old box first */
X mkbox(start_x, start_y, new_x, new_y);
X }
X start_x = new_x = event_x(event);
X start_y = new_y = event_y(event);
X put_coords(new_x, new_x, new_y, new_y);
X mouseing = TRUE;
X } else {
X mouseing = FALSE;
X selected = TRUE;
X }
X break;
X case LOC_DRAG:
X if (mouseing) {
X mkbox(start_x, start_y, new_x, new_y); /* Remove old box */
X
X /* We want to restrict the size to be square */
X tmpx = event_x(event) - start_x;
X tmpy = start_y - event_y(event);
X maxdist = MIN(tmpx * sgn(tmpx), tmpy * sgn(tmpy));
X new_x = start_x + maxdist * sgn(tmpx);
X new_y = start_y - maxdist * sgn(tmpy);
X
X mkbox(start_x, start_y, new_x, new_y); /* Draw new box */
X put_coords(MIN(start_x, new_x), MAX(start_x, new_x),
X MAX(start_y, new_y), MIN(start_y, new_y));
X }
X break;
X case MS_MIDDLE:
X if (selected) {
X mkbox(start_x, start_y, new_x, new_y);
X selected = FALSE;
X }
X }
X}
X
Xvoid stop()
X{
X abort = TRUE;
X}
X
Xvoid up()
X{
X if (sp > 0) {
X Pixrect *old;
X register int i, j, k;
X Bool *mem;
X
X selected = FALSE;
X sp--;
X xmin = Stack[sp].xmin;
X xmax = Stack[sp].xmax;
X ymin = Stack[sp].ymin;
X ymax = Stack[sp].ymax;
X side = xmax - xmin;
X put_coords(0, N-1, N-1, 0);
X old = Stack[sp].pr;
X pw_write(pw, 0, 0, N, N, PIX_SRC, old, 0, 0);
X
X /* Restore MSet */
X /* This is ugly - I'm assuming that sizeof(Bool) == 1. Shame on me! */
X mem = (Bool *) mpr_d(old)->md_image;
X for (i=0; i<N; i++) {
X for (j=0; j<N; j+=8) {
X for (k=0; k<8; k++)
X MSet[j+k][N-1-i] = ((*mem & (1<<(7-k))) >> (7-k)) == 0 ? 1 : 0;
X mem++;
X }
X }
X pr_destroy(old);
X }
X}
X
Xvoid draw()
X{
X draw_new_image = TRUE;
X}
X
Xvoid quit()
X{
X#ifdef oldSunOS
X exit(0); /* You won't miss the following fancy stuff.. */
X#else
X if (alert_prompt(fr, (Event *)0,
X ALERT_MESSAGE_STRINGS,
X "Please confirm:",
X "Do you know what you're doing??",
X 0,
X ALERT_BUTTON_YES, "Of course, quit bugging me!",
X ALERT_BUTTON_NO, "Sorry, I hit the wrong button...",
X ALERT_MESSAGE_FONT, bf,
X ALERT_BUTTON_FONT, pf,
X 0)
X == ALERT_YES) {
X exit(0);
X } else
X notify_veto_destroy(fr);
X#endif
X}
X
Xcompute()
X{
X register int ix, iy;
X
X if (selected && start_x != new_x) { /* New region selected */
X Pixrect *save = mem_create(N, N, 1);
X mkbox(start_x, start_y, new_x, new_y); /* Remove the box first */
X pw_read(save, 0, 0, N, N, PIX_SRC, pw, 0, 0);
X Stack[sp].pr = save;
X Stack[sp].xmin = xmin;
X Stack[sp].xmax = xmax;
X Stack[sp].ymin = ymin;
X Stack[sp].ymax = ymax;
X if (sp < Stacksize) sp++; /* Hard to imagine this happen, but... */
X bzero((char *) MSet, sizeof(MSet));
X pw_writebackground(pw, 0, 0, N, N, PIX_SRC | PIX_COLOR(1));
X
X xmax = xmin + side * MAX(start_x, new_x) /(N-1);
X xmin += side * MIN(start_x, new_x) /(N-1);
X ymax = ymin + side * (N-1-MIN(start_y, new_y)) /(N-1);
X ymin += side * (N-1-MAX(start_y, new_y)) /(N-1);
X selected = FALSE;
X } else {
X /* No region selected, just redraw. Perhaps using new parameters. */
X put_coords(0, N-1, N-1, 0);
X }
X
X rec = (int) panel_get_value(recur_p);
X incr = (int) panel_get_value(incr_p);
X maxiter = (int) panel_get_value(maxiter_p);
X
X side = xmax - xmin;
X delta = 0.25 * side / (N-1); /* 0.25 seems OK */
X recur = rec * delta;
X
X abort = FALSE;
X
X/*************************************************************************/
X/*************************************************************************/
X
X/* From now on, you will find the new Mandelbrot algorithm. No Sun specific
X stuff, except the notify_dispatch() and some pw_put() and pw_line(). */
X
X for (iy = 0; iy < N; iy += incr) {
X notify_dispatch(); /* Allow user to hit the STOP button */
X if (abort) break;
X for (ix = 0; ix < N; ix += incr)
X if (!MSet[ix][iy])
X MDisk(ix,iy);
X }
X}
X
XMDisk(ix, iy)
X register int ix, iy;
X{
X register double cx, cy, dist;
X register int irad;
X
X if (ix<0 || ix>=N || iy<0 || iy>=N || MSet[ix][iy]) return;
X
X cx = xmin + (side * ix) / (N-1);
X cy = ymin + (side * iy) / (N-1);
X dist = 0.25 * MSetDist(cx, cy, maxiter);
X irad = dist / side * (N-1); /* Bug in the original algorithm */
X
X if (irad == 1) {
X MSet[ix][iy] = 1;
X pw_put(pw, ix, N-1-iy, 0); /* Sun specific */
X } else if (irad > 1)
X FILLDISK(ix, iy, irad);
X else if (dist > delta) {
X MSet[ix][iy] = 1;
X pw_put(pw, ix, N-1-iy, 0); /* Sun specific */
X }
X
X if (dist > recur) {
X if (irad > 1) irad++;
X MDisk(ix, iy + irad);
X MDisk(ix, iy - irad);
X MDisk(ix + irad, iy);
X MDisk(ix - irad, iy);
X
X/* It will be slightly faster if I leave out the following "45-degree"
X recursions. The reason is that most of these points will probably
X be filled already and MDisk will return immediately. But since
X they are in the original algorithm and the improvement is only marginal
X I will leave them here. */
X
X irad = 0.5 + irad / SQRT_2;
X MDisk(ix + irad, iy + irad);
X MDisk(ix - irad, iy - irad);
X MDisk(ix - irad, iy + irad);
X MDisk(ix + irad, iy - irad);
X }
X}
X
Xdouble MSetDist(cx, cy, maxiter)
X register double cx, cy;
X register int maxiter;
X{
X# define overflow 10e10 /* Don't know if this is foolproof */
X
X register int iter=0;
X register double zx, zy, zx2, zy2;
X register double *xorbit, *yorbit;
X
X /* Could use a static array for this, if you don't have alloca */
X xorbit = (double *) alloca(maxiter * sizeof(double));
X yorbit = (double *) alloca(maxiter * sizeof(double));
X
X /* This is the standard Mandelbrot iteration */
X zx = zy = zx2 = zy2 = xorbit[0] = yorbit[0] = 0.0;
X do {
X zy = (zx * zy) + (zx * zy) + cy; /* gcc generates only one mult for this */
X zx = zx2 - zy2 + cx;
X iter++;
X xorbit[iter] = zx; /* Save the iteration orbits for later */
X yorbit[iter] = zy;
X zx2 = zx * zx;
X zy2 = zy * zy;
X } while ((zx2 + zy2) < 1000.0 && iter<maxiter);
X
X if (iter < maxiter) { /* Generate derivatives */
X register double zpx, zpy, tmp;
X register int i;
X
X zpx = zpy = 0.0;
X
X for (i=0; i<iter; i++) {
X tmp = 2 * (xorbit[i] * zpx - yorbit[i] * zpy) + 1.0;
X zpy = 2 * (xorbit[i] * zpy + yorbit[i] * zpx);
X zpx = tmp;
X if (fabs(zpx) > overflow || fabs(zpy) > overflow)
X return 0.0;
X }
X /* This is the distance estimation */
X return log(zx2 + zy2) * sqrt(zx2 + zy2) / sqrt(zpx*zpx + zpy*zpy);
X }
X return 0.0;
X}
X
XFILLDISK(ix, iy, irad)
X register int ix, iy;
X int irad;
X{
X register int x, y, e;
X
X /* The "Mini"-algorithm. Perhaps I should use display locking around the
X plotline's, but after all, the fun is watching it work... */
X
X x = 0;
X y = irad;
X e = irad / 2;
X while (x <= y) {
X plotline(ix - x, ix + x, iy + y);
X plotline(ix - y, ix + y, iy + x);
X plotline(ix - x, ix + x, iy - y);
X plotline(ix - y, ix + y, iy - x);
X e -= x;
X if (e < 0) {
X e += y;
X y--;
X }
X x++;
X }
X}
X
Xplotline(x1, x2, y)
X register int x1, x2, y;
X{
X register int i;
X if (y<0 || y>N-1 || (x1<0 && x2<0) || (x1>=N-1 && x2 >=N-1)) return;
X
X if (x1 < 0) x1 = 0;
X if (x1 > N-1) x1 = N-1;
X if (x2 < 0) x2 = 0;
X if (x2 > N-1) x2 = N-1;
X
X pw_vector(pw, x1, N-1-y, x2, N-1-y, PIX_SRC, 0); /* Sun specific */
X
X for (i=x1; i<=x2; i++)
X MSet[i][y] = 1;
X}
SHAR_EOF
chmod 0644 FastMtool.c || echo "restore of FastMtool.c fails"
fi
exit 0