⟦12d3a9b0c⟧

TextFile


/*
 * Mandelzoom (from Scientific American, August 1985).
 * Recommended places for exploration (according to the article):
 *	 real  imaginary	side
 *	-2.00	 -1.25		2.50	Entire Mandelbrot set
 *	  .26	  0.00		 .01
 *	 -.76	  0.01		 .02
 *	-1.26	  0.01		 .02
 *
 * This program computes the Mandelbrot set, writing the pixels
 * to xxx.pix and the parameters (origin, size, and a pixel density
 * histogram to xxx.his).
 *
 * Option command line (for spawning):
 * 	mandel origin_real origin_imag size npixel niter filename
 * or	mandel xxx.cmd			read commands from a file.
 * or	mandel				interactive
 *
 * Values are:
 *  orig_real	(double) lower-left hand corner of the picture (real part)
 *  orig_imag	(double) lower-left hand corner of the picture (imag part)
 *  size	(double) size of a side (defines the view into the set
 *  npixel	(int)    number of pixels on each side (max == 512)
 *				default = 100
 *  niter	(int)    number of iterations per point.
 *				default = 100
 *  file	(string) name of output file (no type)
 *				default = mandel
 *
 * In "mandel file" format, each line generates a separate picture.
 *
 * Each picture causes two files to be written:
 *   file.his	(readable text format) contains picture definition parameters
 *		and a histogram of pixel values:
 *	orig_real orig_imag side npixel niter \n
 *	bottom top sum	(first non-zero histogram, last non-zero,
 *			  sum of histogram values).
 *	hist[bottom]	(number of pixels with count == bottom)
 *	...		(etc.)
 *	hist[top]	(number of pixels with specified count)
 *  file.pix	(binary format) contains npixel rows, each containing npixel
 *		integers, defining the count at each pixel location.
 *
 * Note: this program is very cpu intensive.  There are a maximum of
 *	4 * niter * (npixel**2) floating-point multiplies
 *     10 * niter * (npixel**2) other floating-point (add, compare, store)
 * per picture.
 *
 * Decus C bug note:
 *   all printf's that format floats must have only one floating-point
 *   parameter, and it must be the last parameter in the argument list.
 *   Also, Decus C doesn't support (double) x++.
 */

#include	<stdio.h>
#include	<ctype.h>
#ifdef decus
#define W_BIN_MODE	"wn"		/* Write binary file in Decus C	*/
#else
#define W_BIN_MODE	"w"
#endif
#ifdef vms
#include		errno
#define	IO_ERROR	errno
#define	IO_SUCCESS	1		/* SS$_NORMAL			*/
#endif
#define FALSE	0
#define TRUE	1
#define	EOS	'\0'
#ifndef	IO_ERROR
#define	IO_SUCCESS	0		/* Unix definitions		*/
#define	IO_ERROR	1
#endif

#ifdef decus
int		$$narg = 1;		/* Don't prompt for commands	*/
#endif
#define	MAX_NPIXEL	400		/* Forced by display format	*/
#define	MAX_NITER	1000		/* Max. number of iterations	*/
#define DEF_FILENAME	"mandel"	/* Default output filename	*/
#define	DEF_NITER	100		/* Default number of iterations	*/
#define	DEF_NPIXEL	100		/* Default number of pixels	*/
#define	MAX_ARGS	8		/* Command line arguments	*/

/*
 * These values may be reset by the user.
 */
int		niter = DEF_NITER;	/* Number of iterations		*/
int		npixel = DEF_NPIXEL;	/* Number of pixels/side	*/
double		orig_real = -2.0;	/* South-West (i.e. lower-left)	*/
double		orig_imag = -1.25;	/* corner of the picture	*/
double		side = 2.5;

int		pixel[MAX_NPIXEL];	/* Stores one pixel row		*/
double		gap[MAX_NPIXEL];	/* Real axis position		*/
double		hist[MAX_NITER];	/* Pixel density histogram	*/
char		line[256];		/* General text work area	*/
char		filename[81];		/* Output file name work area	*/
char		*myargv[MAX_ARGS];	/* To build command arguments	*/
int		myargc;			/* Index into myargv[]		*/
FILE		*pixfd;			/* Pixel file (binary)		*/
FILE		*hisfd;			/* Pixel histogram (text)	*/

extern double	atof();
extern char	*strchr();
\f


main(argc, argv)
int		argc;
char		*argv[];
{
	if (argc <= 1) {
	    while (interactive()) {
		doit();
	    }
	}
	else if (argc == 2 && !isdigit(argv[1][0])) {
	    if (freopen(argv[1], "r", stdin) == NULL) {
		perror(argv[1]);
		exit(IO_ERROR);
	    }
	    while (comfile()) {
		getarguments(myargc, myargv);
		doit();
	    }
	}
	else {
	    getarguments(argc, argv);
	    doit();
	}
	exit(IO_SUCCESS);
}

doit()
{
	if (setup()) {
	    process();
	    finish();
	}
}
\f


process()
/*
 * Compute the Mandelbrot set.
 */
{
	register double	z_real, z_imag;	/* Pixel accumulator		*/
	register double	z2_real;	/* z_real ** 2			*/
	register double z2_imag;	/* z_imag ** 2			*/
	register int	count;		/* Inner-loop counter		*/
	register int	j;		/* Column counter		*/
	register int	i;		/* Row counter			*/
	register double	c_real, c_imag;	/* Pixel position (constant)	*/
	register double	float_pixels;	/* To computer pixel position	*/
	register double	total;		/* For progress log		*/

	/*
	 * Precompute the position of each pixel on the real axis.
	 * This loop should not be "unrolled" to a succession of
	 * additions as that would lose accuracy.  Given what
	 * follows, the cost isn't excessive.
	 */
	total = 0.0;			/* Progress log counter		*/
	float_pixels = npixel;
	for (j = 0; j < npixel; j++)	/* Precompute column (real) gap	*/
	    gap[j] = orig_real + (side * ((double) j) / float_pixels);
	for (i = 0; i < npixel; i++) {
	    c_imag = orig_imag + (side * ((double) i) / float_pixels);
	    for (j = 0; j < npixel; j++) {
		/*
	 	 * Compute one point (pixel) of the Mandelbrot set:
		 * 1. Set z to 0 + 0i and
		 *    set c to the pixel location (real, imag)
		 * 2. Perform step 3 until either
		 *    1. count reaches the selected number of iterations or
		 *    2. the "size" of z exceeds 2.0, where size is
		 *       defined as sqrt(z_real**2 + z_imag**2)
		 *       (we don't bother with the square root.)
		 * 3. z = z**2 + c;
		 *    count = count + 1;
		 *    size = size of z as defined above.
		 * 4. Set pixel[i,j] to the number of iterations (count).
		 */
		c_real = gap[j];
		z_real = c_real;
		z_imag = c_imag;
		for (count = 0; count < niter; count++) {
		    if (((z2_real = (z_real * z_real))
		       + (z2_imag = (z_imag * z_imag))) > 4.0)
			break;
		    z_imag = (z_real * z_imag * 2.0) + c_imag;
		    z_real = z2_real - z2_imag + c_real;
		}
		pixel[j] = count;		/* Store the pixel	*/
		hist[count] += 1.0;		/* Gray-scale histogram	*/
		total += (double) count;	/* For progress log	*/
		if (i == j && isatty(fileno(stdout))) {
		    scr_move(7, 1);		/* Log progress @ diag.	*/
		    printf("%3d %3d %4d %8.0f", i, j, count, total);
		    scr_eol();
		    fflush(stdout);
		}
	    }					/* All columns this row	*/
	    fwrite(pixel, npixel * sizeof (pixel[0]), 1, pixfd);
	    if (ferror(pixfd)) {
		perror("pixel file write error");
		return;
	    }
	}					/* All rows in picture	*/
}
\f


int
setup()
/*
 * Open the .pix and .his files.  Zero the vectors.
 */
{
	register int		i;
	extern FILE		*fcreate();

	sprintf(line, "%s.pix", filename);
	if ((pixfd = fopen(line, W_BIN_MODE)) == NULL) {
	    perror(line);
	    sleep(2);
	    return (FALSE);
	}
	sprintf(line, "%s.his", filename);
	if ((hisfd = fcreate(line)) == NULL) {
	    perror(line);
	    sleep(2);
	    fclose(pixfd);
	    return (FALSE);
	}
	fprintf(hisfd, "%f", orig_real);
	fprintf(hisfd, " %f", orig_imag);
	fprintf(hisfd, " %f", side);
	fprintf(hisfd, " %d", npixel);
	fprintf(hisfd, " %d\n", niter);
	for (i = 0; i < niter; i++)
	    hist[i] = 0.0;
	return (TRUE);
}
\f


finish()
/*
 * Compute a histogram of the counts.
 */
{
	register int		bottom, top, i;
	register double		sum;

	fclose(pixfd);
	for (bottom = 0; bottom < niter && hist[bottom] == 0.0; bottom++)
	    ;
	for (top = niter - 1; top > bottom  && hist[top] == 0.0; top--)
	    ;
	for (sum = 0.0, i = bottom; i <= top; i++)
	    sum = sum + hist[i];
	fprintf(hisfd, "%d %d %f\n", bottom, top, sum);
	for (i = bottom; i <= top; i++)
	    fprintf(hisfd, "%10.2f\n", hist[i]);
	fclose(hisfd);	   
}	
\f


getarguments(argcount, argstring)
int		argcount;
char		*argstring[];
/*
 * Process argv[] from command line or pseudo argv[] from a file.
 */
{
	register int		i;

	if (argcount < 5) {
	    printf("Missing arguments, need at least\n");
	    printf("orig_real, orig_imaginary, side\n");
	}
	else {
	    orig_real	= atof(argstring[1]);
	    orig_imag	= atof(argstring[2]);
	    side	= atof(argstring[3]);
	    strcpy(filename, DEF_FILENAME);
	    niter	= DEF_NITER;
	    npixel	= DEF_NPIXEL;
	    switch (argcount) {
	    default:
		printf("Extra arguments ignored:\n");
		for (i = 7; i < argcount; i++)
		    printf(" arg[%d] = \"%s\"\n", i, argstring[i]);
	    case 7:
		strcpy(filename, argstring[6]);
	    case 6:
		if ((niter = atoi(argstring[5])) > MAX_NITER) {
		    printf("%d iterations max.\n", MAX_NITER);
		    niter = MAX_NITER;
		}
	    case 5:
		npixel = atoi(argstring[4]);
	    }
	    if (npixel > MAX_NPIXEL)
		npixel = MAX_NPIXEL;
	    printf("corner = [%f,",	orig_real);
	    printf("%f], ",		orig_imag);
	    printf("side = %f, ",	side);
	    printf("%d pixels, %d iterations\n", npixel, niter);
	}
}
\f


int
interactive()
/*
 * Read commands from the terminal -- assumed to be a VT100 or similar.
 */
{
	if (isatty(fileno(stdin)))
	    scr_clear();
	if (!getcomplex(1, "Lower-left Corner (real, imaginary)",
		&orig_real, &orig_imag)
	 || !getdouble(2,  "Side Length", &side)
	 || !getint(3,"Number of pixels on a side",
		&npixel, DEF_NPIXEL, MAX_NPIXEL)
	 || !getint(4,	"Iterations",
		&niter, DEF_NITER, MAX_NITER)
	 || !getstring(5, "Output filename", filename, DEF_FILENAME))
	    return (FALSE);
#if 0
/* Debug */
	scr_move(2, 1);
	printf("corner = [%f,",	orig_real);
	printf("%f], ",		orig_imag);
	printf("side = %f, ",	side);
	printf("%d pixels, %d iterations\n", npixel, niter);
	scr_eol();
	sleep(2);
#endif
	return (TRUE);
}
\f


int
comfile()
/*
 * Read commands from an indirect command file.
 */
{
	register char	*lp;

	if (gets(line) == NULL)
	    return (FALSE);
	myargv[0] = "";
	for (myargc = 1, lp = line; *lp != EOS && myargc < MAX_ARGS;) {
	    while (isspace(*lp))
		lp++;
	    myargv[myargc++] = lp;
	    while(!isspace(*lp))
		lp++;
	    if (*lp != EOS)
		*lp++ = EOS;
	}
	return (TRUE);
}
DataMuseum.dk

DKUUG/EUUG Conference tapes

⟦12d3a9b0c⟧ TextFile

Derivation

TextFile
