/*
							HUGEMAND.C

	druckt die Mandelbrotmenge im Format 1,55 x 1,25 m streifenweise aus.

						(c) PsychoterrorSoft
*/

/**************************************************************************/
/*								Includes								  */
/**************************************************************************/

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

/**************************************************************************/
/*								Makros									  */
/**************************************************************************/

#define WIDTH		1080
#define HEIGHT		1080
#define NX			9
#define NY			4
#define BPL			(((NX*(WIDTH-1)+1+15)/16)*2)
#define SIZE		(BPL*((HEIGHT-1)*NY+1))
#define OUTFILE		"hugemand.pcx"

#define CELLSIZE 	(1.2/(double)(NY*(HEIGHT-1)))
#define MAXITER		100
#define BOUND		2.0
#define OVERFLOW	1E2400
#define DELTA_C2	3e-5


/**************************************************************************/
/*							globale Variablen							  */
/**************************************************************************/

char *pmap, *ptop;

/**************************************************************************/
/*							Prototypen									  */
/**************************************************************************/

void save(void);
void do_rect(int x0, int y0, int w, int h, int xbase, int stage);
int lbl_demM(double rec, double imc);
int check_frame(int x0, int y0, int x1, int y1);
void mirror_all(void);
void setpoint(int x, int y);
void fillrect(int x0, int y0, int x1, int y1);

/**************************************************************************/
/*							Hauptprogramm								  */
/**************************************************************************/

void main(void)
{
	char choice;
	int stripenr, ycount, firststage;
	short *blank;
	
	pmap= malloc( SIZE );
	if( !pmap )	{
		printf("Not enough memory\n");
		return;
	}
	ptop= pmap + SIZE;
	for( blank= (short*)pmap; blank< (short*)ptop; ++blank)
		*blank= 0L;
	for( stripenr= 0; stripenr< NX; ++stripenr )
	{
		if( !stripenr )
			firststage= 5;
		else
			firststage= 3;
		printf("Calculating stripe %d, page %d...\n", stripenr, 0 );
		do_rect( stripenr*(WIDTH-1), (HEIGHT-1)*(NY-1), WIDTH, HEIGHT, 0, firststage );
		--firststage;
		for( ycount= 1; ycount< NY; ++ycount )
		{
			printf("Calculating stripe %d, page %d...\n", stripenr, ycount );
			do_rect( stripenr*(WIDTH-1), (HEIGHT-1)*(NY-1-ycount), WIDTH, HEIGHT, 0, firststage );
		}
	}
	save();
	free(pmap);
}


/**************************************************************************/
/*								Unterprogramme							  */
/**************************************************************************/

/*------------------------------------------------------------------------*/
/*		void save()	speichert die Grafik als PCX-Datei					  */
/*																		  */
/*		->	nichts														  */
/*		<-	nichts														  */
/*------------------------------------------------------------------------*/

void save(void)
{
	char scan= 0;
	int count;
	FILE *handle;
	struct pcx {
	char	id, vers, compr, bits;
	short	xmin, ymin, xmax, ymax, xdpi, ydpi;
	char	coltable[48];
	char	reserved, nplanes;
	short	bpl, palette;
	char	filler[58];
	} header= { 0x0a, 3,0,1, 0,0,NX*(WIDTH-1),2*NY*(HEIGHT-1)+1, 180,180, 
				{0,0,0}, 0, 1, BPL, 1, {0} };
	
	if( !(handle= fopen(OUTFILE, "w")) )
		printf("File error\n");
	else if( fwrite(&header, 1, sizeof(struct pcx), handle)!=sizeof(struct pcx) )	{
		fclose(handle);
		printf("File error\n");
	}
	else if( fwrite(pmap, 1, SIZE, handle)!=SIZE )	{
		fclose(handle);
		printf("File error\n");
	}
	else {
	  mirror_all();
	  if( fwrite(pmap, 1, SIZE, handle)!=SIZE )
		printf("File error\n");
	  fclose(handle);
	}
	
}


/*------------------------------------------------------------------------*/
/*		void do_rect()	berechnet rekursiv die Punkte eines Rechtecks	  */
/*																		  */
/*			->	int x0,y0	linke untere Ecke							  */
/*				int w,h		Groesse des Rechtecks						  */
/*				int xbase	Abstand der Position x==0 vom Anfang des	  */
/*							ganzen Bildes								  */
/*				int stage	Highbyte: 1==naechste Unterteilung in obere/  */
/*							untere Haelfte, 0==linke/rechte				  */
/*							Lowbyte: Welche Seiten schon berechnet?		  */
/*							5: nichts, 4: oben, 3: links, 2: links+oben	  */
/*							1: alle ausser der Abgrenzung zur jeweils	  */
/*							zweiten Haelfte (s. Highbyte), 0: alle		  */
/*			<-	nichts													  */
/*------------------------------------------------------------------------*/

void do_rect(int x0, int y0, int w, int h, int xbase, int stage)
{
	static double re, im;
	static int count;
	int check;
	
	if( (char)stage )
		if( (char)stage==1 )
			if( stage==0x101 )
			{
				im= y0*CELLSIZE + CELLSIZE/2;
				re= -2.1 + CELLSIZE*(xbase+x0+1);
				for( count= x0+1; count< x0+w-1; ++count, re += CELLSIZE )
					if( lbl_demM(re, im)==2 )
						setpoint(count, y0);
			}
			else
			{
				im= (long)(y0+1)*CELLSIZE + CELLSIZE/2;
				re= -2.1 + CELLSIZE*(xbase+x0+w-1);
				for( count= y0+1; count< y0+h-1; ++count, im += CELLSIZE )
					if( lbl_demM(re, im)==2 )
						setpoint(x0+w-1, count);
			}
		else
		{
			im= (long)y0*CELLSIZE + CELLSIZE/2;
			re= -2.1 + CELLSIZE*(xbase+x0+1);
			for( count= x0+1; count< x0+w; ++count, re += CELLSIZE )
				if( lbl_demM(re, im)==2 )
					setpoint(count, y0);
			im += CELLSIZE;	re -= CELLSIZE;
			for( count= y0+1; count< y0+h-1; ++count, im += CELLSIZE )
				if( lbl_demM(re, im)==2 )
					setpoint(x0+w-1, count);
			if( stage & 1 )
			{
				for( count= x0+w-1; count> x0; --count, re -= CELLSIZE )
					if( lbl_demM(re, im)==2 )
						setpoint(count, y0+h-1);
				if( (char)stage>3 && lbl_demM(re, im)==2 )
						setpoint(count, y0+h-1);
			}
			else if( (char)stage>3 )
				re -= CELLSIZE*(long)(w-1);
			if( (char)stage>3 )
			{
				im -= CELLSIZE;
				for( count= y0+h-2; count>= y0; --count, im -= CELLSIZE )
					if( lbl_demM(re, im)==2 )
						setpoint(x0, count);
			}
		}	/*	(char)stage > 1  */ 
	
	if( w<3 || h<3 )		return;
	
	if( (check= check_frame(x0, y0, x0+w-1, y0+h-1))!=0 )
		if( check==2 )
			fillrect(x0,y0, x0+w-1, y0+h-1);
		else;
	else
		if( stage & 0x100 )
		{
			do_rect( x0, y0, w/2+1, h, xbase, 1);
			do_rect( x0+w/2, y0, w-w/2, h, xbase, 0);
		}
		else
		{
			do_rect( x0, y0+h/2, w, h-h/2, xbase, 0x101);
			do_rect( x0, y0, w, h/2+1, xbase, 0x100);
		}
	
}	/*	do_rect()  */


/*------------------------------------------------------------------------*/
/*		int lbl_demM()	fuehrt den distance estimation Algorithmus fuer	  */
/*						die Mandelbrotmenge aus							  */
/*																		  */
/*		->	double rec, imc		Parameter c								  */
/*		<-	0 c Element von M, +/-1 nahe dran, 2 weder noch				  */
/*------------------------------------------------------------------------*/

int lbl_demM(double rec, double imc)
{
	char flag;
	int count, end;
	double rez, imz, buf1, buf2, modcsqr;
	static struct {
		double re, im;
	}	orbit[MAXITER];
	
	orbit->re= rez= rec;
	orbit->im= imz= imc;
	modcsqr= rec*rec+imc*imc;
	if( modcsqr< 0.25*0.25 )
		return 0;
	if( (buf1= rez*rez)+(buf2= imz*imz) > BOUND*BOUND )
		end= 0;
	else
	{
		for( count= 1; count< MAXITER; ++count )
		{
			buf1= buf1-buf2 + rec;
			buf2= rez*imz;
			orbit[count].re= rez= buf1;
			orbit[count].im= imz= buf2+buf2 + imc;
			if( 0!=(flag= ((buf1= rez*rez)+(buf2= imz*imz) > BOUND*BOUND)) )
				break;
		}
		if( !flag )		return(0);
		else	end= count;
	}
	
	rez= rec+rec+1.0;
	imz= imc+imc;
	if( rez*rez+imz*imz>OVERFLOW*OVERFLOW )	return(-1);
	for( count= 1; count<end; ++count )
	{
		buf1= rez*orbit[count].re-imz*orbit[count].im;
		buf2= rez*orbit[count].im+imz*orbit[count].re;
		rez= buf1+buf1+1.0;
		imz= buf2+buf2;
		if( rez*rez+imz*imz>OVERFLOW*OVERFLOW )	return(-1);
	}
	buf1= orbit[end].re*orbit[end].re+orbit[end].im*orbit[end].im;
	buf2= (buf1+buf1)/(rez*rez+imz*imz);
	buf1= buf2 * log(buf1);
	if( buf1<DELTA_C2*modcsqr )	return(1);
	else	return(2);
	
}	/*	lbl_demM()  */


/*------------------------------------------------------------------------*/
/*		void mirror_all()	spiegelt die Druckerbitmap					  */
/*																		  */
/*			->	nichts													  */
/*			<-	nichts													  */
/*------------------------------------------------------------------------*/

void mirror_all(void)
{
	int count;
	char *upper= pmap, *lower= ptop-BPL, tmp;
	
	while( lower>upper ) 
	{
		for( count= BPL; count> 0; --count )
		{
			tmp= *upper;
			*upper++= *lower;
			*lower++ = tmp;
		}
		lower -= 2*BPL;
	}
	
}	/*	mirror_all()  */


/*------------------------------------------------------------------------*/
/*		int check_frame()	Prueft, ob ein Rechteckrahmen einfarbig ist.  */
/*																		  */
/*		->	int x0,y0	linke untere Ecke								  */
/*			int x1,y1	rechte obere Ecke								  */
/*		<-	0 mehrfarbig, 1 schwarz, 2 weiss							  */
/*------------------------------------------------------------------------*/

int check_frame(int x0, int y0, int x1, int y1)
{
	char *leflow, *rilow, *lefup, *riup, *check;
	char andval, byteflag;
	
	leflow= ptop - BPL - BPL*y0 + (x0>>3);
	lefup= leflow - BPL*(y1-y0);
	/*	falls nur schmaler Streifen, der innerhalb eines Bytes liegt:  */
	if( (x0>>3) == (x1>>3) )
		if( !(*leflow & (andval= (char)(0xFF>>(x0&7)) & (0xFF<<(7-(x1&7))))) )	/* schwarz */
		{
			/*	oberer Rand:  */
			if( *lefup & andval )	return(0);
			/*	Rest:  */
			andval= (0x80>>(x0&7)) | (0x80>>(x1&7));
			for( check= leflow-BPL; check> lefup; check -= BPL)
				if( *check & andval )	return(0);
			return(1);
		}
		else if( !(~*leflow & andval) )				/* weiss */
		{
			/*	oberer Rand:  */
			if( ~*lefup & andval )	return(0);
			/*	Rest:  */
			andval= (0x80>>(x0&7)) | (0x80>>(x1&7));
			for( check= leflow-BPL; check> lefup; check -= BPL)
				if( ~*check & andval )	return(0);
			return(2);
		}
		else return(0);
	else
	{
		rilow= ptop - BPL - BPL*y0 + (x1>>3);
		riup= rilow - BPL*(y1-y0);
		if( !(*leflow & (andval= 0xFF>>(x0&7))) )		/*	weiss  */
		{
			/*	linke obere Ecke:  */
			if( (*lefup & andval)!=0 )	return(0);
			/*	linker Rand:  */
			andval= 0x80>>(x0&7);
			for( check= leflow - BPL; check> lefup; check -= BPL )
				if( *check & andval )	return(0);
			/*	rechter Rand:  */
			andval= 0x80>>(x1&7);
			for( check= rilow - BPL; check> riup; check -= BPL )
				if( *check & andval )	return(0);
			/*	rechte Ecken:  */
			andval= 0xFF<<(7-(x1&7));
			if( (*rilow & andval)!=0 )	return(0);
			if( (*riup & andval)!=0 )	return(0);
			/*	unterer Rand:  */
			for( check= leflow + 1; check< rilow; ++check )
				if( *check )	return(0);
			/*	oberer Rand:  */
			for( check= lefup + 1; check< riup; ++check )
				if( *check )	return(0);
			return(1);
		}
		else if( !(~*leflow & andval) )		/*	schwarz  */
		{
			/*	linke obere Ecke:  */
			if( (~*lefup & andval)!=0 )	return(0);
			/*	linker Rand:  */
			andval= 0x80>>(x0&7);
			for( check= leflow - BPL; check> lefup; check -= BPL )
				if( ~*check & andval )	return(0);
			/*	rechter Rand:  */
			andval= 0x80>>(x1&7);
			for( check= rilow - BPL; check> riup; check -= BPL )
				if( ~*check & andval )	return(0);
			/*	rechte Ecken:  */
			andval= 0xFF<<(7-(x1&7));
			if( (~*rilow & andval)!=0 )	return(0);
			if( (~*riup & andval)!=0 )	return(0);
			/*	unterer Rand:  */
			for( check= leflow + 1; check< rilow; ++check )
				if( ~*check )	return(0);
			/*	oberer Rand:  */
			for( check= lefup + 1; check< riup; ++check )
				if( ~*check )	return(0);
			return(2);
		}
		else return(0);
	}
	
}	/*	check_frame()  */


/*------------------------------------------------------------------------*/
/*		void setpoint()	setzt einen weissen Punkt in die Druckerbitmap    */
/*																		  */
/*			->	int x,y		Koordinaten des Punktes (0,0 ist links unten) */
/*			<-	nichts													  */
/*------------------------------------------------------------------------*/

void setpoint(int x, int y)
{
	if( y<0 || x<0 || y>(HEIGHT-1)*NY || x>(WIDTH-1)*NX )	return;
	
	ptop[-BPL - BPL*y + (x>>3)] |= 0x80>>(x&7);
	
}	/*	setpoint()  */


/*------------------------------------------------------------------------*/
/*		void fillrect()	fuellt in der Druckerbitmap ein Rechteck weiss   */
/*							aus; der Rahmen wird nicht veraendert.		  */
/*																		  */
/*		->	int x0,y0	linke untere Ecke								  */
/*			int x1,y1	rechte obere Ecke								  */
/*		<-	nichts														  */
/*------------------------------------------------------------------------*/

void fillrect(int x0, int y0, int x1, int y1)
{
	char *boc, *eoc, *black;
	int count, nfull;
	char orval;
	
	if( y1-y0 < 2 || x1-x0 < 2 )	return;
	++x0;	++y0;	--x1;
	/*	falls nur schmaler Streifen, der innerhalb eines Bytes liegt:  */
	if( (x0>>3) == (x1>>3) )
	{
		orval= (char)(0xFF>>(x0&7)) & (0xFF<<(7-(x1&7)));
		black= ptop - BPL - BPL*y0 + (x0>>3);
		eoc= black - BPL*(y1-y0);
		for( ; black> eoc; black -= BPL)
			*black |= orval;
	}
	else
	{
		nfull= ((x1+1)>>3) - ((x0-1)>>3) - 1;
		boc= ptop - BPL - y0*BPL + (x0>>3);
		eoc= boc - BPL*(y1-y0);
		if( (x0&7)!=0 )	/*	Streifen am linken Rand des Rechtecks:  */
		{
			orval= 0xFF>>(x0&7);
			for( black= boc; black> eoc; black-=BPL )
				*black |= orval;
			++boc;	++eoc;
		}
		for( count= nfull; count> 0; --count )
		{
			for( black= boc; black> eoc; black -= BPL )
				*black= 0xFF;
			++boc;	++eoc;
		}
		if( ((x1+1)&7)!=0 )	/*	Streifen am oberen Rand des Rechtecks:  */
		{
			orval= (char)(0xFF<<(7-(x1&7)));
			for( black= boc; black> eoc; black -= BPL )
				*black |= orval;
		}
	}	/*	nicht nur schmaler Streifen  */
	
}	/*	fillrect()  */