#include <math.h>

	/* Perlin Noise Function -- from Graphics Gems II */

#if 0
#define hermite(p0,p1,r0,r1,t) (	p0*((2.*t-3.)*t*t+1.) + \
									p1*(-2.*t+3.)*t*t + \
									r0*((t-2.)*t+1.)*t + \
									r1*(t-1.)*t*t )

#define dhermite(p0,p1,r0,r1,t) (	p0*(t-1.)*6.*t + \
									p1*(-t+1.)*6.*t + \
									r0*((3.*t-4.)*t + 1.) + \
									r1*(3.*t-2.)*t )
#else
#define hermite(p0,p1,r0,r1,t) (	(p0-p1)*((2.*t-3.)*t*t) + p0 \
									r0*((t-2.)*t+1.)*t + \
									r1*(t-1.)*t*t )

#define dhermite(p0,p1,r0,r1,t) (	(p0-p1)*(t-1.)*6.*t + \
									r0*((3.*t-4.)*t + 1.) + \
									r1*(3.*t-2.)*t )
#endif

#define rand3a(x,y,z)	frand(67*(x)+59*(y)+71*(z))
#define rand3b(x,y,z)	frand(73*(x)+79*(y)+83*(z))
#define rand3c(x,y,z)	frand(89*(x)+97*(y)+101*(z))
#define rand3d(x,y,z)	frand(103*(x)+107*(y)+109*(z))

short					xlim[3][2];				/* integer bounds for point		*/
short					xint[3];				/* fractional part				*/
double					xarg[3];				/* fractional part				*/

#if 0
double frand(long s)
{	s = s<<13 ^ s;
	return (1.0 - ((s*(s*s*15731+789221)+1376312589)&0x7fffffff)
					/ 1073741824.0 );
}
#else
double frand(long s)
{	return (float)(rand24(s) - 128) * .0078125;	/* / 128 */
}
#endif

double calc_hermite(double p0,double p1,double r0,double r1,double t)
{	return (p0*((2.*t-3.)*t*t+1.) +
			p1*(-2.*t+3.)*t*t +
			r0*((t-2.)*t+1.)*t +
			r1*(t-1.)*t*t );
}

double calc_dhermite(double p0,double p1,double r0,double r1,double t)
{	return (p0*(t-1.)*6.*t +
			p1*(-t+1.)*6.*t +
			r0*((3.*t-4.)*t + 1.) +
			r1*(3.*t-2.)*t );
}

extern char				rtable0[256],
						rtable1[256],
						rtable2[256],
						rtable3[256];

long					xyzrand16t(short,short,short,char *);

#define HERMITE_TABSIZ	1024
#define HERMITE_SCALE	0x0400
#define HERMITE_UNITY	0x8000

struct HermiteTable {
	short				cp0,cp1,cr0,cr1,		/* hermite coefficients			*/
						dp0,dp1,dr0,dr1;		/* hermite derivative coeff.	*/
	float				ep0,ep1,er0,er1,
						fp0,fp1,fr0,fr1;
	short				t;
};

static struct HermiteTable h_table[HERMITE_TABSIZ];
struct HermiteTable		*h_arg[4];

init_hermite()
{	float				t;
	short				i;
	struct HermiteTable	*ht;

	for (i=0; i<HERMITE_TABSIZ; i++)
	{	ht = &h_table[i];

		t = (float)i / (float)HERMITE_TABSIZ;

		ht->cp0 = ((2.*t-3.)*t*t+1.) * (float)HERMITE_SCALE;
		ht->cp1 = (-2.*t+3.)*t*t * (float)HERMITE_SCALE;
		ht->cr0 = ((t-2.)*t+1.)*t * (float)HERMITE_SCALE;
		ht->cr1 = (t-1.)*t*t * (float)HERMITE_SCALE;

		ht->dp0 = (t-1.)*6.*t * (float)HERMITE_SCALE;
		ht->dp1 = (-t+1.)*6.*t * (float)HERMITE_SCALE;
		ht->dr0 = ((3.*t-4.)*t + 1.) * (float)HERMITE_SCALE;
		ht->dr1 = (3.*t-2.)*t * (float)HERMITE_SCALE;
#if 0
		ht->ep0 = ((2.*t-3.)*t*t+1.);
		ht->ep1 = (-2.*t+3.)*t*t;
		ht->er0 = ((t-2.)*t+1.)*t;
		ht->er1 = (t-1.)*t*t;

		ht->fp0 = (t-1.)*6.*t;
		ht->fp1 = (-t+1.)*6.*t;
		ht->fr0 = ((3.*t-4.)*t);
		ht->fr1 = (3.*t-2.)*t;
#endif
		ht->t	= (int)(t * 4096.0);
	}
}

long lerp_i(short p1, short p2, short t)
{	return p1 + (((p2 - p1) * t) >> 12);
}

#if 0
interpolate(double f[4], short e[4], int i, int n)
{
	if (n == 0)
	{	long			ix = xlim[0][i&1],
						iy = xlim[1][i>>1&1],
						iz = xlim[2][i>>2];

#if 0
		f[0] = (float)(xyzrand16t(ix,iy,iz,rtable0) - 128);
		f[1] = (float)(xyzrand16t(ix,iy,iz,rtable1) - 128);
		f[2] = (float)(xyzrand16t(ix,iy,iz,rtable2) - 128);
		f[3] = (float)(xyzrand16t(ix,iy,iz,rtable3) - 128);
#endif

		e[0] = ((xyzrand16t(ix,iy,iz,rtable0) - 128) << 5) + (1<<4);
		e[1] = ((xyzrand16t(ix,iy,iz,rtable1) - 128) << 5) + (1<<4);
		e[2] = ((xyzrand16t(ix,iy,iz,rtable2) - 128) << 5) + (1<<4);
		e[3] = ((xyzrand16t(ix,iy,iz,rtable3) - 128) << 5) + (1<<4);

#if 0
		f[0] = rand3a(ix,iy,iz);
		f[1] = rand3b(ix,iy,iz);
		f[2] = rand3c(ix,iy,iz);
		f[3] = rand3d(ix,iy,iz);
#endif
		return;
	}
	else
	{	double			f0[4],					/* result for first half		*/
						f1[4];					/* result for second half		*/
		double			t,
						nt;
		short			e0[4],
						e1[4];
		struct HermiteTable	*ht;

		n--;
		interpolate(f0,e0,i,n);
		interpolate(f1,e1,i|1<<n,n);

		ht = h_arg[n];

		e[0] = lerp_i( e0[0], e1[0], ht->t );
		e[1] = lerp_i( e0[1], e1[1], ht->t );
		e[2] = lerp_i( e0[2], e1[2], ht->t );

		e[3] = (e0[3] * ht->cp0 +
				e1[3] * ht->cp1 +
				e0[n] * ht->cr0 +
				e1[n] * ht->cr1) >> 10;

/*		e[3] = (int)(f[3] * 32.0); */

#if 0
		t = xarg[n];
		nt = 1.0 - t;

		f[0] = nt*f0[0] + t*f1[0];
		f[1] = nt*f0[1] + t*f1[1];
		f[2] = nt*f0[2] + t*f1[2];

/*		f[3] = calc_hermite(f0[3], f1[3], f0[n], f1[n], t); */
/*		f[n] = calc_dhermite(f0[3], f1[3], f0[n], f1[n], t); */
#endif
	}
}
#else
interpolate(short e[4], int i, int n)
{
	if (n == 0)
	{	long			ix = xlim[0][i&1],
						iy = xlim[1][i>>1&1],
						iz = xlim[2][i>>2];

		e[0] = ((xyzrand16t(ix,iy,iz,rtable0) - 128) << 5) + (1<<4);
		e[1] = ((xyzrand16t(ix,iy,iz,rtable1) - 128) << 5) + (1<<4);
		e[2] = ((xyzrand16t(ix,iy,iz,rtable2) - 128) << 5) + (1<<4);
		e[3] = ((xyzrand16t(ix,iy,iz,rtable3) - 128) << 5) + (1<<4);
		return;
	}
	else
	{	short			e0[4],
						e1[4];
		struct HermiteTable	*ht;

		n--;
		interpolate(e0,i,n);
		interpolate(e1,i|1<<n,n);

		ht = h_arg[n];

		e[0] = lerp_i( e0[0], e1[0], ht->t );
		e[1] = lerp_i( e0[1], e1[1], ht->t );
		e[2] = lerp_i( e0[2], e1[2], ht->t );

		e[3] = (e0[3] * ht->cp0 +
				e1[3] * ht->cp1 +
				e0[n] * ht->cr0 +
				e1[n] * ht->cr1) >> 10;
	}
}
#endif

double *noise3(double x[3])
{	static double		f[4];
	short				e[4];
	long				l;

	l = x[0] * (float)HERMITE_TABSIZ;
	xlim[0][0] = l >> 10; xlim[0][1] = xlim[0][0] + 1;
	h_arg[0] = &h_table[l & (HERMITE_TABSIZ - 1)];

	l = x[1] * (float)HERMITE_TABSIZ;
	xlim[1][0] = l >> 10; xlim[1][1] = xlim[1][0] + 1;
	h_arg[1] = &h_table[l & (HERMITE_TABSIZ - 1)];

	l = x[2] * (float)HERMITE_TABSIZ;
	xlim[2][0] = l >> 10; xlim[2][1] = xlim[2][0] + 1;
	h_arg[2] = &h_table[l & (HERMITE_TABSIZ - 1)];

#if 0
	xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
	xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
	xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1;
	xarg[0] = x[0] - xlim[0][0];
	xarg[1] = x[1] - xlim[1][0];
	xarg[2] = x[2] - xlim[2][0];
	h_arg[0] = &h_table[ (int)(xarg[0] * HERMITE_TABSIZ) ];
	h_arg[1] = &h_table[ (int)(xarg[1] * HERMITE_TABSIZ) ];
	h_arg[2] = &h_table[ (int)(xarg[2] * HERMITE_TABSIZ) ];
#endif
	interpolate(e,0,3);
/*	f[0] /= 128.0;
	f[1] /= 128.0;
	f[2] /= 128.0;
	f[3] /= 128.0;
*/
	f[3] = (float)e[3] / (32.0 * 128.0);
	return (f);
}

double *noise2(double x[2])
{	static double		f[4];
	short				e[4];

	xlim[0][0] = floor(x[0]); xlim[0][1] = xlim[0][0] + 1;
	xlim[1][0] = floor(x[1]); xlim[1][1] = xlim[1][0] + 1;
/*	xlim[2][0] = floor(x[2]); xlim[2][1] = xlim[2][0] + 1; */
	xarg[0] = x[0] - xlim[0][0];
	xarg[1] = x[1] - xlim[1][0];
/*	xarg[2] = x[2] - xlim[2][0]; */
	h_arg[0] = &h_table[ (int)(xarg[0] * HERMITE_TABSIZ) ];
	h_arg[1] = &h_table[ (int)(xarg[1] * HERMITE_TABSIZ) ];
	h_arg[2] = &h_table[ (int)(xarg[2] * HERMITE_TABSIZ) ];
	interpolate(e,0,2);
	f[0] /= 128.0;
	f[1] /= 128.0;
	f[2] /= 128.0;
	f[3] /= 128.0;
	return (f);
}