/* noise function over R3 - implemented by a pseudorandom tricubic spline */
/* Ken Perlin 12/89 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "multicolor.h"


#define B 0x100
#define BM 0xff
#define N 0x10000000
#define NP 12   /* 2^N */
#define NM 0xfff

#define ADOT(a,b) (a[0] * b[0] + a[1] * b[1] + a[2] * b[2])
#define setup(u,b0,b1,r0,r1)\
	t = u + N;\
	b0 = ((int)t) & BM;\
	b1 = (b0+1) & BM;\
	r0 = t - (int)t;\
	r1 = r0 - 1.;


static int	p[B + B + 2];
static double	g[B + B + 2][3];
static int	start = 1;


double Noise3( Vector vec )
{
	int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
	double rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
	register i, j;

	if (start) {
		start = 0;
		init();
	}

	setup(vec.x, bx0,bx1, rx0,rx1);
	setup(vec.y, by0,by1, ry0,ry1);
	setup(vec.z, bz0,bz1, rz0,rz1);

	i = p[ bx0 ];
	j = p[ bx1 ];

	b00 = p[ i + by0 ];
	b10 = p[ j + by0 ];
	b01 = p[ i + by1 ];
	b11 = p[ j + by1 ];

#define at(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] )

#define s_curve(t) ( t * t * (3. - 2. * t) )

#define lerp(t, a, b) ( a + t * (b - a) )

	t  = s_curve(rx0);
	sy = s_curve(ry0);
	sz = s_curve(rz0);

	q = g[ b00 + bz0 ] ; u = at(rx0,ry0,rz0);
	q = g[ b10 + bz0 ] ; v = at(rx1,ry0,rz0);
	a = lerp(t, u, v);

	q = g[ b01 + bz0 ] ; u = at(rx0,ry1,rz0);
	q = g[ b11 + bz0 ] ; v = at(rx1,ry1,rz0);
	b = lerp(t, u, v);

	c = lerp(sy, a, b);		/* interpolate in y at lo x */

	q = g[ b00 + bz1 ] ; u = at(rx0,ry0,rz1);
	q = g[ b10 + bz1 ] ; v = at(rx1,ry0,rz1);
	a = lerp(t, u, v);

	q = g[ b01 + bz1 ] ; u = at(rx0,ry1,rz1);
	q = g[ b11 + bz1 ] ; v = at(rx1,ry1,rz1);
	b = lerp(t, u, v);

	d = lerp(sy, a, b);		/* interpolate in y at hi x */

	return 1.5 * lerp(sz, c, d);	/* interpolate in z */
} /* Noise3() */


void init( void )
{
	long random();
	int i, j, k;
	double v[3], s;

	/* Create an array of random vectors uniformly on the unit sphere */
	srandom(1);
	for (i = 0 ; i < B ; i++) {
		do {				/* Choose uniformly in a cube */
			for (j=0 ; j<3 ; j++)
				v[j] = (double)((random() % (B + B)) - B) / B;
			s = ADOT(v,v);
		} while (s > 1.0);		/* If not in sphere try again */
		s = sqrt(s);
		for (j = 0 ; j < 3 ; j++)	/* Else normalize */
			g[i][j] = v[j] / s;
	}

	/* Create a pseudorandom permutation of [1..B] */
	for (i = 0 ; i < B ; i++)
		p[i] = i;
	for (i = B ; i > 0 ; i -= 2) {
		k = p[i];
		p[i] = p[j = random() % B];
		p[j] = k;
	}

	/* Extend g and p arrays to allow for faster indexing */
	for (i = 0 ; i < B + 2 ; i++) {
		p[B + i] = p[i];
		for (j = 0 ; j < 3 ; j++)
			g[B + i][j] = g[i][j];
	}
} /* init() */


Vector VecNoise3( Vector point )
{
	Vector 	result;

	result.x = Noise3( point );
	point.x += 10.0;
	result.y = Noise3( point );
	point.y += 10.0;
	result.z = Noise3( point );

	return( result );
} /* VecNoise3() */