/************************************************************************/
/* Float.c								*/
/*									*/
/* Here are found all of the routines doing the math and control logic	*/
/* for "Apfelkiste" using regular float calculations. Even these should	*/
/* be faster than any other "normal" Mandelbrot set generating program	*/
/* that uses the standard line-by-line algorithm. There is a little	*/
/* more overhead but profiling the program shows that more than 90% of	*/
/* calculation time is done in the function that does the math, so	*/
/* every point that is saved from beeing evaluated means saved time.	*/
/************************************************************************/

#include <intuition/intuition.h>
#include <graphics/gfx.h>
#include <graphics/gfxmacros.h>
#include <proto/intuition.h>
#include <proto/graphics.h>
#include <mffp.h>
#include <math.h>

extern struct RastPort	*rp;
extern int		 maxiter;
extern double		 dx, dy, divergenz, rmin, imin;
extern short		 maxcol;

/* -------------------------------------------------------------------- */
/* Evaluate color of a given point X0 on the complex number plane (hope	*/
/* this is the correct term) by the formula:				*/
/*									*/
/*                                    2					*/
/*                            X    = X  + X				*/
/*                             n+1    n    0				*/
/*									*/
/* and counting the iterations until Xn exceedes an upper boundary	*/
/* "divergenz" or a maximum count "maxiter" is reached. The number of	*/
/* iterations performed determines the color of the point to set.	*/
/* -------------------------------------------------------------------- */

long Iter_FLOAT(double r, double i)
{
	int	count = 0;
	double	x, p, q;

	p = r; q = i;
	while ( r * r + i * i < divergenz && count++ < maxiter ) {
		x = r;
		r = (x + i) * (x - i) + p;
		i = 2.0 * x * i + q;
	}
	if ( count >= maxiter ) {
		return 0;
	}
	else {
		return (long) count % maxcol;
	}
}

/* -------------------------------------------------------------------- */
/* Determine whether a rectangular area is surrounded by all the same	*/
/* color. In case of this no point inside the area will get a different */
/* color and thus there's no need to do any calculations for points	*/
/* inside.								*/
/* This function is also used by the fixpoint version of the algorithm.	*/
/* -------------------------------------------------------------------- */

short __regargs Test(long xx1, long yy1, long xx2, long yy2)
{
	long i;
	int  c;

	c = ReadPixel(rp, xx1, yy1);

	for ( i = xx2; i >= xx1; i-- ) {
		if ( ReadPixel(rp, i, yy1) != c ) return FALSE;
		if ( ReadPixel(rp, i, yy2) != c ) return FALSE;
	}

	for ( i = yy2 - 1; i > yy1; i-- ) {
		if ( ReadPixel(rp, xx1, i) != c ) return FALSE;
		if ( ReadPixel(rp, xx2, i) != c ) return FALSE;
	}
	return TRUE;
}

/* -------------------------------------------------------------------- */
/* Recursive evaluation scheme for "Apfelkiste". Not very much glasnost	*/
/* because of diverse optimizations (?), but briefly it does this:	*/
/* Called with the corners of an already calculated rectangular frame	*/
/* it tests for identity of the color of all points on that frame.	*/
/* If yes, fine! Just fill the points inside the area with this color	*/
/* and return without further work.					*/
/* If not, the area is halved and the function calls itself once for	*/
/* each half until the size of the area degenerates to zero or the	*/
/* first termination criterium jumps in.				*/
/* Watch the execution of the program and you'll see.			*/
/* -------------------------------------------------------------------- */

void __regargs Apfel_FLOAT(int x1, int y1, int x2, int y2)
{
	int xneu, yneu;
	long i;
	double help1, help2;

	if ( x2 - x1 > 1 && y2 - y1 > 1 ) {
		if ( Test(x1, y1, x2, y2) ) {
			SetAPen(rp, ReadPixel(rp, x1, y1));
			RectFill(rp, x1, y1, x2, y2);
		}
		else {
			if ( x2 - x1 > y2 - y1 ) {
				xneu = ( x1 + x2 ) >> 1;
				help1 = rmin + xneu * dx;
				help2 = ( y1 + 1 ) * dy + imin;
				for ( i = y1 + 1; i < y2; i++ ) {
					SetAPen(rp, Iter_FLOAT(help1, help2));
					WritePixel(rp, xneu, i);
					help2 += dy;
				}
				Apfel_FLOAT(x1, y1, xneu, y2);
				Apfel_FLOAT(xneu, y1, x2, y2);
			}
			else {
				yneu = ( y1 + y2 ) >> 1;
				help1 = imin + yneu * dy;
				help2 = ( x1 + 1 ) * dx + rmin;
				for ( i = x1+1; i < x2; i++ ) {
					SetAPen(rp, Iter_FLOAT(help2, help1));
					WritePixel(rp, i, yneu);
					help2 += dx;
				}
				Apfel_FLOAT(x1, y1, x2, yneu);
				Apfel_FLOAT(x1, yneu, x2, y2);
			}
		}
	}
}