#ifndef lint
static char *sccsid = "%W% (UKC) %G%";
#endif  lint
/*
 * we keep this routine in a separate file so that it may be optimised,
 * since this is where the program spends all its time.
 */

#include "mandel.h"

/*
 * 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.
 *
 * To 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).
 *
 *	Operations on fixed-point quantities:
 *	fixed + fixed => fixed
 *	fixed - fixed => fixed
 *	fixed * int   => fixed
 *	fixed / fixed => int
 *	fixed / int   => fixed
 *
 *	Other operations must be pre- and/or post-scaled.
 *
 * We have a fractional part 24 bits long. The need to munge by 4 reduces this
 * to 22 bits. For a 512-wide picture, maximum side length is 1/(2^22/512)
 * == 1/8192 == 0.000122
 */

typedef long fixed;	/* A signed integral type */
#define NBITS	(((sizeof(fixed)*8)-8)/2*2)
#define ONE	(1<<NBITS)		/* to hold -127..+127 */
#define SQRT1	(1<<(NBITS/2))
#define MULT(a,b) (((a)/SQRT1)*((b)/SQRT1))	/* fixed * fixed => fixed */
#define DTOF(a)	((fixed) (a * ONE))

static	short	pixel[MAX_NPIXEL];	/* Stores one pixel row		*/
static	fixed	gap[MAX_NPIXEL];	/* Real axis position		*/
extern	int	hist[];			/* Pixel density histogram	*/
extern	char	line[];			/* General text work area	*/
extern	FILE	*pixfd;			/* Pixel file (binary)		*/

/*
 * Compute the Mandelbrot set.
 */
process()
{
	register int	count;		/* Inner-loop counter		*/
	int	j;			/* Column counter		*/
	int	i;			/* Row counter			*/
	register fixed	z2_real;	/* z_real ** 2			*/
	register fixed	z2_imag;	/* z_imag ** 2			*/
	register fixed	z_real, z_imag;	/* Pixel accumulator		*/
	register fixed	c_real, c_imag;	/* Pixel position (constant)	*/
	fixed	o_real = DTOF(start_real);	/* private copy */
	fixed	o_imag = DTOF(start_imag);	/* private copy */
	fixed	o_side = DTOF(side);		/* private copy */

	/* Precompute column (real) gap	*/
	/* fixed can hold up to 128, MAXNPIXEL*max(o_real) == 512*2 => predivide by 8 */
	for (count = 0; count < npixel; count++)
	    gap[count] = o_real + (o_side / 8 * count / npixel * 8);

	for (i = 0; i < npixel; i++) {
	    c_imag = o_imag + (o_side / 8 * i / npixel * 8);
	    for (j = 0; j < npixel; j++) {
		c_real = gap[j];
		z_real = c_real;
		z_imag = c_imag;
		for (count = 0; count < niter; count++) {
#ifdef DEBUG
fprintf(stderr, "%4d %10d (%2.5f)  %10d (%2.5f) ", count, z_real, (double)z_real/ONE, z_imag, (double)z_imag/ONE);
#endif
		    z2_real = z_real/SQRT1; z2_imag = z_imag/SQRT1; /* temps */
		    z_imag = z2_real * z2_imag;		/* == z_real * z_imag */
		    z2_real = z2_real * z2_real;	/* == z_real squared */
		    z2_imag = z2_imag * z2_imag;	/* == z_imag squared */
#ifdef DEBUG
fprintf(stderr, "%10d (%2.5f)\n", z2_real+z2_imag, (double)(z2_real+z2_imag)/ONE);
#endif
		    if (z2_real + z2_imag > DTOF(4))
			break;
		    z_imag = z_imag * 2 + c_imag;
		    z_real = z2_real - z2_imag + c_real;
		}
		pixel[j] = count;		/* Store the pixel	*/
		hist[count] += 1;		/* Gray-scale histogram	*/
	    }					/* 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	*/
}