/* Test program for new-ra 
 * From ljapunov calc program
 * (c) Kurt Garloff, 10/2003, GNU GPL
 */ 

#define _GNU_SOURCE
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#if 1 //def __alpha__
#include <setjmp.h>
#include <signal.h>
#include <fenv.h>
#endif


#define LONG_DOUBLE long double
#define MAX(a,b) ((a)>(b)?(a):(b))
#define MIN(a,b) ((a)>(b)?(b):(a))
#if defined(__GNUC__) && __GNUC__ >= 3
# define LIKELY(x)	__builtin_expect((x) != 0,1) 
# define UNLIKELY(x)	__builtin_expect((x) != 0,0)
#else
# define LIKELY(x)	(x)
# define UNLIKELY(x)	(x)
#endif

#define MAXPATLEN 64
#define MAXCOLENTS 12
// Coltrafo description	

typedef struct _ct_inf {	
	float val;
	float h, s, i;
} ct_inf;

enum calcmode { FIXX, FIXB, FIXA };

typedef struct _lja_param {
	char magic[4];	/* "LJA0" */
	char rev;	/* revision */
	char cmode;	/* enum calcmode actually */
	unsigned char pass;
	unsigned char patlen;
/*0x08:*/
	double x0;
/*0x10:*/
	double amin, amax, bmin, bmax;
/*0x30:*/
	unsigned int xres, yres, bline, depth;
/*0x40:*/
	char pattern[MAXPATLEN];
/*0x80:*/
	/* Info for lja2bmp */
	float ampln, amplp;
/*0x88:*/
	char res2[55];
	char ctents;
/*0x0c0:*/
	ct_inf coltbl [MAXCOLENTS];
/*0x180:*/
} lja_parm;
	
lja_parm *ljp;
float *results;
unsigned long int nl;
// Count of A's and B's in pattern
static int na, nb;
static char *fullpat;
char quiet;
volatile sig_atomic_t interrupt;


#define INF 20000
#define SML 1e-24
#define DIV 160.0

const char lja_hdr_latest_rev = 2;

void lja_defaults (lja_parm * const ljp)
{
	memset (ljp, 0, sizeof (lja_parm));
	/* Ljapunov settings */
	memcpy (ljp->magic, "LJA0", 4);
	ljp->rev = lja_hdr_latest_rev; ljp->cmode = FIXX;
	ljp->pass = 0; 	ljp->patlen = 2;
	ljp->x0 = 0.4;
	ljp->amin = -2.7; ljp->amax = 4.2;
	ljp->bmin = -2.7; ljp->bmax = 4.2;
	ljp->xres = 1280; ljp->yres = 1024;
	ljp->bline = 0; ljp->depth = 256;
	memset (ljp->pattern, 0, MAXPATLEN); ljp->pattern[1] = 1;
	ljp->ampln = 0.4; ljp->amplp = 1.6;
	ljp->ctents = 0;
}

static //inline 
double extrapol_lja (const double a, const double b, const double l, 
		     const unsigned r, const unsigned dpth)
{
	++nl;
	if (UNLIKELY(l <= 0.0))
		return -INF;
	else
		return ( log(l) + r * 
			(  log (fabs (a)) * na
			 + log (fabs (b)) * nb ) ) / dpth;
}

#if 1 //def __alpha__
/* Exception handling: Will probably not work with multithreading */
jmp_buf lja_fpe_jump;

void fpe_handler (int signum, int info)
{
	if (!info)
		info = -1;
	printf ("\x1b[D!"); fflush (stdout); 
	if (lja_fpe_jump)
		longjmp (lja_fpe_jump, info);
	else
		raise (signum);
}
#endif
	
static inline
void set_pat (double* pat, const int ln, const double a, const double b)
{
	int i;
	for (i = 0; i < ln; ++i)
	       	pat[i] = (fullpat[i%ln]? b: a);
}

double ljapunov (const double a, const double b, const double xi, 
		 const unsigned dpth)
{
	int t = 0; register int p;
	const int exppln = (ljp->patlen > 16?     ljp->patlen: 
			    ljp->patlen >  8? 2 * ljp->patlen: 
			    		      4 * ljp->patlen);
	register LONG_DOUBLE l, r, x = xi;
	const int rep = dpth / exppln;
	
	double pattern[32];
	set_pat (pattern, exppln, a, b);
	//double pattern[2];
	//pattern[0] = a; pattern[1] = b;
	
	l = 1.0;
#if 1 //def __alpha__
	if (UNLIKELY(setjmp (lja_fpe_jump)))
		return extrapol_lja (a, b, l, (dpth-t*exppln)/ljp->patlen, dpth);
#endif
	//asm (".align 32\n");
	for (t = 0; t < rep; ++t) {
		for (p = 0; p < exppln; ++p) {
			//r = (fullpat[p]? b: a);
			r = pattern[p];
			//r = pattern[fullpat[p]];
			x *= r*(1.0-x); 
			l *= fabs(r*(1.0-x*2.0));
		}
		{
		register const double fabsx = fabs(x);
		if ( /*t >= 2 &&*/
		    UNLIKELY(fabsx > DIV || fabsx < SML) )
			return extrapol_lja (a, b, l, (dpth-t*exppln)/ljp->patlen, dpth);

		}
	}
	// Residuals ?
	t = dpth - rep*exppln; 
	for (p = 0; p < t; ++p) {
		//r = (fullpat[p]? b : a);
		r = pattern[p];
		//r = pattern[fullpat[p]];
		x *= r*(1.0-x); 
		l *= fabs (r*(1.0-x*2.0));
	}

	if (LIKELY(l > 0.0)) 
		return log (l) / dpth;
	else 
		return -INF;
/*
 out:
	return extrapol_lja (a, b, l, (dpth-t*exppln)/ljp->patlen, dpth);
 */	
}

void prepare_ljapunov (const int dpth)
{
	// pattern count
	int s; na = 0; nb = 0;
	for (s = 0; s < ljp->patlen; ++s) {
		if (ljp->pattern[s])
			++nb;
		else
			++na;
	}
	for (s = 0; s < dpth; ++s)
		fullpat[s] = ljp->pattern[s%ljp->patlen];
}

void do_calc_lines (const unsigned xstep, const unsigned ystep, 
		    const unsigned depth,
		    const unsigned start, const unsigned end)
{
	unsigned ln, col;
	const double mx = (ljp->amax - ljp->amin)/ljp->xres;
	const double my = (ljp->bmax - ljp->bmin)/ljp->yres;
#ifdef __alpha__
	fenv_t fpu;
	//fegetenv(&fpu);
	//printf("FPU: 0x%016lx", fpu);
	//fesetenv(FE_DFL_ENV);
	fegetenv(&fpu);
	//printf(" 0x%016lx", fpu);
	fpu |= FE_MAP_DMZ | FE_MAP_UMZ;
	fesetenv(&fpu);
	fedisableexcept(FE_ALL_EXCEPT);
	//fegetenv(&fpu);
	//printf(" 0x%016lx\n", fpu);
#endif	
	signal (SIGFPE, (sig_t)fpe_handler);
	
	for (ln = start; ln < end; ln += ystep) {
		double y = (double)ljp->bmax - my * (ln + (double)ystep*0.5);
		float res;
		//printf ("%i (%ix%i)\n", ln, xstep, ystep);
		switch (ljp->cmode) {
		    case FIXX:
			for (col = 0; col < ljp->xres; col += xstep) {
				double x = (double)ljp->amin + mx * (col + (double)xstep*0.5);
				register unsigned ln2, col2;
				res = ljapunov (x, y, ljp->x0, depth);
				for (ln2 = ln; ln2 < MIN(ln+ystep,end); ln2++)
					for (col2 = col; col2 < MIN(col+xstep,ljp->xres-1); ++col2)
						*(results+col2+ljp->xres*ln2) = res;
			}; break;
		    case FIXA:
			for (col = 0; col < ljp->xres; col += xstep) {
				double x = (double)ljp->amin + mx * (col + (double)xstep*0.5);
				register unsigned ln2, col2;
				res = ljapunov (ljp->x0, y, x, depth); 
				for (ln2 = ln; ln2 < MIN(ln+ystep,end); ++ln2)
					for (col2 = col; col2 < MIN(col+xstep,ljp->xres-1); col2++)
						*(results+col2+ljp->xres*ln2) = res;
			}; break;
		    case FIXB:
			for (col = 0; col < ljp->xres; col += xstep) {
				double x = (double)ljp->amin + mx * (col + (double)xstep*0.5);
				register unsigned ln2, col2;
				res = ljapunov (x, ljp->x0, y, depth);
				for (ln2 = ln; ln2 < MIN(ln+ystep,end); ++ln2)
					for (col2 = col; col2 < MIN(col+xstep,ljp->xres-1); col2++)
						*(results+col2+ljp->xres*ln2) = res;
			}; break;
		    default: res = 0; abort ();
		}
	}
	signal (SIGFPE, SIG_DFL);
}

void calcn (const unsigned xstep, const unsigned ystep, 
	    const unsigned dpth)
{
	unsigned int y = 0; int thr;
	const int ychunk = 32;
	unsigned yrmax = ljp->yres;
	fullpat = (char*)malloc (dpth); prepare_ljapunov (dpth);
	for (y = ljp->bline; y < yrmax && !interrupt; y += ychunk) {
		if (!quiet) { 
			printf ("."); fflush (stdout); 
		}
		/* Do the last chunk ourself */
		do_calc_lines (xstep, ystep, dpth, 
				y, y+ychunk);
	}
	if (interrupt) { 
	   	if (!quiet) 
			printf ("%i", y-1); 
	   	ljp->bline = y; 
		free (fullpat);
		return; 
	}
	if (y != yrmax)
		do_calc_lines (xstep, ystep, dpth, y, yrmax);
	//ljp->bline = ljp->yres;
	ljp -> bline = 0;
	free (fullpat);
	if (!quiet) { 
		printf ("%i\n", ljp->yres-1); fflush (stdout); 
	};
}

int main (int argc, char* argv[])
{
	clock_t st, en;
	st = clock();
	ljp = (lja_parm*)malloc(sizeof(lja_parm));
	lja_defaults(ljp);
	results = (float*) malloc(ljp->xres * ljp->yres * sizeof (float));
	calcn (1, 1, 256);
	free(results);
	free(ljp);
	en = clock();
	printf("CPU time: %.3fs\n", (double)(en-st)/CLOCKS_PER_SEC);
	return 0;
}