#include <math.h>

static double 
dot3d(int g[], double x, double y, double z)
{
	return g[0] * x + g[1] * y + g[2] * z;
}

static double 
dot2d(int g[], double x, double y)
{
	return g[0] * x + g[1] * y;
}

static double 
dot4d(int g[], double x, double y, double z, double w)
{
	return g[0] * x + g[1] * y + g[2] * z + g[3] * w;
}

static double 
fade(double t)
{
	return t * t * t * (t * (t * 6 - 15) + 10);
}

static double 
mix(double a, double b, double t)
{
	return (1 - t) * a + t * b;
}

static int      grad3[][3] = {
	{1, 1, 0}, {-1, 1, 0}, {1, -1, 0}, {-1, -1, 0},
	{1, 0, 1}, {-1, 0, 1}, {1, 0, -1}, {-1, 0, -1},
	{0, 1, 1}, {0, -1, 1}, {0, 1, -1}, {0, -1, -1}
};

static int      grad4[][4] = {{0, 1, 1, 1}, {0, 1, 1, -1}, {0, 1, -1, 1}, {0, 1, -1, -1},
{0, -1, 1, 1}, {0, -1, 1, -1}, {0, -1, -1, 1}, {0, -1, -1, -1},
{1, 0, 1, 1}, {1, 0, 1, -1}, {1, 0, -1, 1}, {1, 0, -1, -1},
{-1, 0, 1, 1}, {-1, 0, 1, -1}, {-1, 0, -1, 1}, {-1, 0, -1, -1},
{1, 1, 0, 1}, {1, 1, 0, -1}, {1, -1, 0, 1}, {1, -1, 0, -1},
{-1, 1, 0, 1}, {-1, 1, 0, -1}, {-1, -1, 0, 1}, {-1, -1, 0, -1},
{1, 1, 1, 0}, {1, 1, -1, 0}, {1, -1, 1, 0}, {1, -1, -1, 0},
{-1, 1, 1, 0}, {-1, 1, -1, 0}, {-1, -1, 1, 0}, {-1, -1, -1, 0}
};

/* Classic Perlin noise in 3D, for comparison */

static int      perm[] = {
	151, 160, 137, 91, 90, 15,
	131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
	190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
	88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166,
	77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244,
	102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196,
	135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123,
	5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42,
	223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
	129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228,
	251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107,
	49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
	138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180,

	/*
	 * To remove the need for index wrapping, double the permutation
	 * table length
	 */

	151, 160, 137, 91, 90, 15,
	131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23,
	190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
	88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166,
	77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244,
	102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196,
	135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123,
	5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42,
	223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
	129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228,
	251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107,
	49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254,
	138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
};

/* Classic Perlin noise, 3D version */
static double 
noise3dp(double x, double y, double z)
{
	/* Find unit grid cell containing point */
	int             X = floor(x);
	int             Y = floor(y);
	int             Z = floor(z);

	/* Get relative xyz coordinates of point within that cell */
	x = x - X;
	y = y - Y;
	z = z - Z;

	/*
	 * Wrap the integer cells at 255 (smaller integer period can be
	 * introduced here)
	 */
	X = X & 255;
	Y = Y & 255;
	Z = Z & 255;

	/* Calculate a set of eight hashed gradient indices */
	int             gi000 = perm[X + perm[Y + perm[Z]]] % 12;
	int             gi001 = perm[X + perm[Y + perm[Z + 1]]] % 12;
	int             gi010 = perm[X + perm[Y + 1 + perm[Z]]] % 12;
	int             gi011 = perm[X + perm[Y + 1 + perm[Z + 1]]] % 12;
	int             gi100 = perm[X + 1 + perm[Y + perm[Z]]] % 12;
	int             gi101 = perm[X + 1 + perm[Y + perm[Z + 1]]] % 12;
	int             gi110 = perm[X + 1 + perm[Y + 1 + perm[Z]]] % 12;
	int             gi111 = perm[X + 1 + perm[Y + 1 + perm[Z + 1]]] % 12;

	/*
         * The gradients of each corner are now:
         * g000 = grad3[gi000];
         * g001 = grad3[gi001];
         * g010 = grad3[gi010];
         * g011 = grad3[gi011];
         * g100 = grad3[gi100];
         * g101 = grad3[gi101];
         * g110 = grad3[gi110];
         * g111 = grad3[gi111];
         */

	/* Calculate noise contributions from each of the eight corners */
	double          n000 = dot3d(grad3[gi000], x, y, z);
	double          n100 = dot3d(grad3[gi100], x - 1, y, z);
	double          n010 = dot3d(grad3[gi010], x, y - 1, z);
	double          n110 = dot3d(grad3[gi110], x - 1, y - 1, z);
	double          n001 = dot3d(grad3[gi001], x, y, z - 1);
	double          n101 = dot3d(grad3[gi101], x - 1, y, z - 1);
	double          n011 = dot3d(grad3[gi011], x, y - 1, z - 1);
	double          n111 = dot3d(grad3[gi111], x - 1, y - 1, z - 1);

	/* Compute the fade curve value for each of x, y, z */
	double          u = fade(x);
	double          v = fade(y);
	double          w = fade(z);

	/* Interpolate along x the contributions from each of the corners */
	double          nx00 = mix(n000, n100, u);
	double          nx01 = mix(n001, n101, u);
	double          nx10 = mix(n010, n110, u);
	double          nx11 = mix(n011, n111, u);

	/* Interpolate the four results along y */
	double          nxy0 = mix(nx00, nx10, v);
	double          nxy1 = mix(nx01, nx11, v);

	/* Interpolate the two last results along z */
	double          nxyz = mix(nxy0, nxy1, w);

	return nxyz;
}

/*
 * A lookup table to traverse the simplex around a given point in 4D.
 * Details can be found where this table is used, in the 4D noise method.
 */
int             simplex[][4] = {
	{0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0},
	{0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0},
	{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
	{1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0},
	{1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0},
	{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
	{2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0},
{2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}};

double grad1(int hash) {
	double g;
	int h = hash & 15;
	g = 1.0 + (h & 7);   /* Gradient value is one of 1.0, 2.0, ..., 8.0 */
	if (h&8) g = -g;   /* Make half of the gradients negative */
	return g;
}

double noise1d(double x, double *dx)
{
	int i0 = floor(x);
	int i1 = i0 + 1;
	double x0 = x - i0;
	double x1 = x0 - 1.0;

	double gx0, gx1;
	double n0, n1;
	double t20, t40, t21, t41;

	double x20 = x0*x0;
	double t0 = 1.0 - x20;

	/*  if(t0 < 0.0f) t0 = 0.0f; /* Never happens for 1D: x0<=1 always */
	t20 = t0 * t0;
	t40 = t20 * t20;
	gx0 = grad1(perm[i0 & 0xff]);
	n0 = t40 * gx0 * x0;

	double x21 = x1*x1;
	double t1 = 1.0 - x21;
	/*  if(t1 < 0.0f) t1 = 0.0f; /* Never happens for 1D: |x1|<=1 always */
	t21 = t1 * t1;
	t41 = t21 * t21;
	gx1 = grad1(perm[i1 & 0xff]);
	n1 = t41 * gx1 * x1;

	/*
	 * Compute derivative according to:
	 * *dnoise_dx = -8.0f * t20 * t0 * x0 * (gx0 * x0) + t40 * gx0;
	 * *dnoise_dx += -8.0f * t21 * t1 * x1 * (gx1 * x1) + t41 * gx1;
	 */

	if (dx) {
		*dx = t20 * t0 * gx0 * x20;
		*dx += t21 * t1 * gx1 * x21;
		*dx *= -8.0;
		*dx += t40 * gx0 + t41 * gx1;
		*dx *= 0.25; /* Scale derivative to match the noise scaling */
	}

	/*
	 * The maximum value of this noise is 8*(3/4)^4 = 2.53125
	 * A factor of 0.395 would scale to fit exactly within [-1,1], but
	 * to better match classic Perlin noise, we scale it down some more.
	 * scale = 1.0 / 2.53125
	 */
	return (n0 + n1)/2.53125;
}

/* 2D simplex noise */
static double 
noise2d(double xin, double yin)
{
	double          n0, n1, n2;	/* Noise contributions from the three
					 * corners */

	/* Skew the input space to determine which simplex cell we're in */
	double          F2 = 0.5 * (sqrt(3.0) - 1.0);
	double          s = (xin + yin) * F2;	/* Hairy factor for 2D */

	int             i = floor(xin + s);
	int             j = floor(yin + s);

	double          G2 = (3.0 - sqrt(3.0)) / 6.0;
	double          t = (i + j) * G2;
	double          X0 = i - t;	/* Unskew the cell origin back to
					 * (x,y) space */
	double          Y0 = j - t;
	double          x0 = xin - X0;	/* The x,y distances from the cell
					 * origin */
	double          y0 = yin - Y0;

	/* For the 2D case, the simplex shape is an equilateral triangle. */

	/* Determine which simplex we are in. */
	int             i1, j1;	/* Offsets for second (middle) corner of
				 * simplex in (i,j) coords */
	if (x0 > y0) {
		/* lower triangle, XY order: (0,0)->(1,0)->(1,1) */
		i1 = 1;
		j1 = 0;
	} else {
		/* upper triangle, YX order: (0,0)->(0,1)->(1,1) */
		i1 = 0;
		j1 = 1;
	}
	/*
         * A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
         * a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where
         * c = (3-sqrt(3))/6
         */
	double          x1 = x0 - i1 + G2;	/* Offsets for middle corner
						 * in (x,y) unskewed coords */
	double          y1 = y0 - j1 + G2;
	double          x2 = x0 - 1.0 + 2.0 * G2;	/* Offsets for last
							 * corner in (x,y)
							 * unskewed coords */
	double          y2 = y0 - 1.0 + 2.0 * G2;

	/* Work out the hashed gradient indices of the three simplex corners */
	int             ii = i & 255;
	int             jj = j & 255;
	int             gi0 = perm[ii + perm[jj]] % 12;
	int             gi1 = perm[ii + i1 + perm[jj + j1]] % 12;
	int             gi2 = perm[ii + 1 + perm[jj + 1]] % 12;

	/* Calculate the contribution from the three corners */
	double          t0 = 0.5 - x0 * x0 - y0 * y0;
	if (t0 < 0) {
		n0 = 0.0;
	} else {
		t0 *= t0;
		n0 = t0 * t0 * dot2d(grad3[gi0], x0, y0);	/* (x,y) of grad3 used
								 * for 2D gradient */
	}

	double          t1 = 0.5 - x1 * x1 - y1 * y1;
	if (t1 < 0) {
		n1 = 0.0;
	} else {
		t1 *= t1;
		n1 = t1 * t1 * dot2d(grad3[gi1], x1, y1);
	}

	double          t2 = 0.5 - x2 * x2 - y2 * y2;
	if (t2 < 0) {
		n2 = 0.0;
	} else {
		t2 *= t2;
		n2 = t2 * t2 * dot2d(grad3[gi2], x2, y2);
	}

	/*
	 * Add contributions from each corner to get the final noise value.
	 * The result is scaled to return values in the interval [-1,1].
	 */
	return 70.0 * (n0 + n1 + n2);
}

/* 3D simplex noise */
static double 
noise3d(double xin, double yin, double zin)
{
	double          n0, n1, n2, n3;	/* Noise contributions from the four
					 * corners */
	/* Skew the input space to determine which simplex cell we're in */
	double          F3 = 1.0 / 3.0;
	double          s = (xin + yin + zin) * F3;	/* Very nice and simple
							 * skew factor for 3D */
	int             i = floor(xin + s);
	int             j = floor(yin + s);
	int             k = floor(zin + s);
	double          G3 = 1.0 / 6.0;	/* Very nice and simple unskew
					 * factor, too */
	double          t = (i + j + k) * G3;
	double          X0 = i - t;	/* Unskew the cell origin back to
					 * (x,y,z) space */
	double          Y0 = j - t;
	double          Z0 = k - t;
	double          x0 = xin - X0;	/* The x,y,z distances from the cell
					 * origin */
	double          y0 = yin - Y0;
	double          z0 = zin - Z0;

	/*
	 * For the 3D case, the simplex shape is a slightly irregular
	 * tetrahedron.
	 */

	/* Determine which simplex we are in. */

	int             i1, j1, k1;	/* Offsets for second corner of
					 * simplex in (i,j,k) coords */

	int             i2, j2, k2;	/* Offsets for third corner of
					 * simplex in (i,j,k) coords */
	if (x0 >= y0) {
		if (y0 >= z0) {
			i1 = 1;
			j1 = 0;
			k1 = 0;
			i2 = 1;
			j2 = 1;
			k2 = 0;
		}
		 /* X Y Z order */ 
		else if (x0 >= z0) {
			i1 = 1;
			j1 = 0;
			k1 = 0;
			i2 = 1;
			j2 = 0;
			k2 = 1;
		}
		 /* X Z Y order */ 
		else {
			i1 = 0;
			j1 = 0;
			k1 = 1;
			i2 = 1;
			j2 = 0;
			k2 = 1;
		}		/* Z X Y order */
	} else {		/* x0<y0 */
		if (y0 < z0) {
			i1 = 0;
			j1 = 0;
			k1 = 1;
			i2 = 0;
			j2 = 1;
			k2 = 1;
		}
		 /* Z Y X order */ 
		else if (x0 < z0) {
			i1 = 0;
			j1 = 1;
			k1 = 0;
			i2 = 0;
			j2 = 1;
			k2 = 1;
		}
		 /* Y Z X order */ 
		else {
			i1 = 0;
			j1 = 1;
			k1 = 0;
			i2 = 1;
			j2 = 1;
			k2 = 0;
		}		/* Y X Z order */
	}

	/*
	 * A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in
	 * (x,y,z),
	 */

	/*
	 * a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in
	 * (x,y,z), and
	 */

	/*
	 * a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in
	 * (x,y,z), where
	 */

	/* c = 1/6. */
	double          x1 = x0 - i1 + G3;	/* Offsets for second corner
						 * in (x,y,z) coords */
	double          y1 = y0 - j1 + G3;
	double          z1 = z0 - k1 + G3;
	double          x2 = x0 - i2 + 2.0 * G3;	/* Offsets for third
							 * corner in (x,y,z)
							 * coords */
	double          y2 = y0 - j2 + 2.0 * G3;
	double          z2 = z0 - k2 + 2.0 * G3;
	double          x3 = x0 - 1.0 + 3.0 * G3;	/* Offsets for last
							 * corner in (x,y,z)
							 * coords */
	double          y3 = y0 - 1.0 + 3.0 * G3;
	double          z3 = z0 - 1.0 + 3.0 * G3;

	/* Work out the hashed gradient indices of the four simplex corners */
	int             ii = i & 255;
	int             jj = j & 255;
	int             kk = k & 255;
	int             gi0 = perm[ii + perm[jj + perm[kk]]] % 12;
	int             gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12;
	int             gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12;
	int             gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12;

	/* Calculate the contribution from the four corners */
	double          t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
	if (t0 < 0)
		n0 = 0.0;
	else {
		t0 *= t0;
		n0 = t0 * t0 * dot3d(grad3[gi0], x0, y0, z0);
	}

	double          t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
	if (t1 < 0) {
		n1 = 0.0;
	}
	else {
		t1 *= t1;
		n1 = t1 * t1 * dot3d(grad3[gi1], x1, y1, z1);
	}

	double          t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
	if (t2 < 0) {
		n2 = 0.0;
	}
	else {
		t2 *= t2;
		n2 = t2 * t2 * dot3d(grad3[gi2], x2, y2, z2);
	}

	double          t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
	if (t3 < 0) {
		n3 = 0.0;
	}
	else {
		t3 *= t3;
		n3 = t3 * t3 * dot3d(grad3[gi3], x3, y3, z3);
	}
	/* Add contributions from each corner to get the final noise value. */
	/* The result is scaled to stay just inside [-1,1] */
	return 32.0 * (n0 + n1 + n2 + n3);
}

/* 4D simplex noise */
double 
noise4d(double x, double y, double z, double w)
{

	/* The skewing and unskewing factors are hairy again for the 4D case */
	double          F4 = (sqrt(5.0) - 1.0) / 4.0;
	double          G4 = (5.0 - sqrt(5.0)) / 20.0;
	double          n0, n1, n2, n3, n4;	/* Noise contributions from
						 * the five corners */
	/*
	 * Skew the (x,y,z,w) space to determine which cell of 24 simplices
	 * we're in
	 */
	double          s = (x + y + z + w) * F4;	/* Factor for 4D skewing */
	int             i = floor(x + s);
	int             j = floor(y + s);
	int             k = floor(z + s);
	int             l = floor(w + s);
	double          t = (i + j + k + l) * G4;	/* Factor for 4D
							 * unskewing */
	double          X0 = i - t;	/* Unskew the cell origin back to
					 * (x,y,z,w) space */
	double          Y0 = j - t;
	double          Z0 = k - t;
	double          W0 = l - t;
	double          x0 = x - X0;	/* The x,y,z,w distances from the
					 * cell origin */
	double          y0 = y - Y0;
	double          z0 = z - Z0;
	double          w0 = w - W0;
	/*
	 * For the 4D case, the simplex is a 4D shape I won't even try to
	 * describe.
	 */
	/*
	 * To find out which of the 24 possible simplices we're in, we need
	 * to
	 */
	/* determine the magnitude ordering of x0, y0, z0 and w0. */
	/*
	 * The method below is a good way of finding the ordering of x,y,z,w
	 * and
	 */
	/* then find the correct traversal order for the simplex we’re in. */
	/*
	 * First, six pair-wise comparisons are performed between each
	 * possible pair
	 */
	/*
	 * of the four coordinates, and the results are used to add up binary
	 * bits
	 */
	/* for an integer index. */
	int             c1 = (x0 > y0) ? 32 : 0;
	int             c2 = (x0 > z0) ? 16 : 0;
	int             c3 = (y0 > z0) ? 8 : 0;
	int             c4 = (x0 > w0) ? 4 : 0;
	int             c5 = (y0 > w0) ? 2 : 0;
	int             c6 = (z0 > w0) ? 1 : 0;
	int             c = c1 + c2 + c3 + c4 + c5 + c6;
	int             i1, j1, k1, l1;	/* The integer offsets for the second
					 * simplex corner */
	int             i2, j2, k2, l2;	/* The integer offsets for the third
					 * simplex corner */
	int             i3, j3, k3, l3;	/* The integer offsets for the fourth
					 * simplex corner */
	/*
	 * simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some
	 * order.
	 */
	/*
	 * Many values of c will never occur, since e.g. x>y>z>w makes x<z,
	 * y<w and x<w
	 */
	/*
	 * impossible. Only the 24 indices which have non-zero entries make
	 * any sense.
	 */
	/*
	 * We use a thresholding to set the coordinates in turn from the
	 * largest magnitude.
	 */
	/*
	 * The number 3 in the "simplex" array is at the position of the
	 * largest coordinate.
	 */
	i1 = simplex[c][0] >= 3 ? 1 : 0;
	j1 = simplex[c][1] >= 3 ? 1 : 0;
	k1 = simplex[c][2] >= 3 ? 1 : 0;
	l1 = simplex[c][3] >= 3 ? 1 : 0;
	/*
	 * The number 2 in the "simplex" array is at the second largest
	 * coordinate.
	 */
	i2 = simplex[c][0] >= 2 ? 1 : 0;
	j2 = simplex[c][1] >= 2 ? 1 : 0;
	k2 = simplex[c][2] >= 2 ? 1 : 0;
	l2 = simplex[c][3] >= 2 ? 1 : 0;
	/*
	 * The number 1 in the "simplex" array is at the second smallest
	 * coordinate.
	 */
	i3 = simplex[c][0] >= 1 ? 1 : 0;
	j3 = simplex[c][1] >= 1 ? 1 : 0;
	k3 = simplex[c][2] >= 1 ? 1 : 0;
	l3 = simplex[c][3] >= 1 ? 1 : 0;
	/*
	 * The fifth corner has all coordinate offsets = 1, so no need to
	 * look that up.
	 */
	double          x1 = x0 - i1 + G4;	/* Offsets for second corner
						 * in (x,y,z,w) coords */
	double          y1 = y0 - j1 + G4;
	double          z1 = z0 - k1 + G4;
	double          w1 = w0 - l1 + G4;
	double          x2 = x0 - i2 + 2.0 * G4;	/* Offsets for third
							 * corner in (x,y,z,w)
							 * coords */
	double          y2 = y0 - j2 + 2.0 * G4;
	double          z2 = z0 - k2 + 2.0 * G4;
	double          w2 = w0 - l2 + 2.0 * G4;
	double          x3 = x0 - i3 + 3.0 * G4;	/* Offsets for fourth
							 * corner in (x,y,z,w)
							 * coords */
	double          y3 = y0 - j3 + 3.0 * G4;
	double          z3 = z0 - k3 + 3.0 * G4;
	double          w3 = w0 - l3 + 3.0 * G4;
	double          x4 = x0 - 1.0 + 4.0 * G4;	/* Offsets for last
							 * corner in (x,y,z,w)
							 * coords */
	double          y4 = y0 - 1.0 + 4.0 * G4;
	double          z4 = z0 - 1.0 + 4.0 * G4;
	double          w4 = w0 - 1.0 + 4.0 * G4;
	/* Work out the hashed gradient indices of the five simplex corners */
	int             ii = i & 255;
	int             jj = j & 255;
	int             kk = k & 255;
	int             ll = l & 255;
	int             gi0 = perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32;
	int             gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32;
	int             gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32;
	int             gi3 = perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32;
	int             gi4 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32;
	/* Calculate the contribution from the five corners */
	double          t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
	if (t0 < 0) {
		n0 = 0.0;
	}
	else {
		t0 *= t0;
		n0 = t0 * t0 * dot4d(grad4[gi0], x0, y0, z0, w0);
	}
	double          t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
	if (t1 < 0) {
		n1 = 0.0;
	}
	else {
		t1 *= t1;
		n1 = t1 * t1 * dot4d(grad4[gi1], x1, y1, z1, w1);
	}
	double          t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
	if (t2 < 0) {
		n2 = 0.0;
	}
	else {
		t2 *= t2;
		n2 = t2 * t2 * dot4d(grad4[gi2], x2, y2, z2, w2);
	}
	double          t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
	if (t3 < 0) {
		n3 = 0.0;
	}
	else {
		t3 *= t3;
		n3 = t3 * t3 * dot4d(grad4[gi3], x3, y3, z3, w3);
	}
	double          t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
	if (t4 < 0) {
		n4 = 0.0;
	}
	else {
		t4 *= t4;
		n4 = t4 * t4 * dot4d(grad4[gi4], x4, y4, z4, w4);
	}

	/* Sum up and scale the result to cover the range [-1,1] */
	return 27.0 * (n0 + n1 + n2 + n3 + n4);
}

double noise_simplex(int dim, double scale, double *pt) {
	switch (dim) {
		case 1:
			return noise1d(pt[0]*scale,0);
			break;
		case 2:
			return noise2d(pt[0]*scale,pt[1]*scale);
			break;
		case 3:
			return noise3d(pt[0]*scale,pt[1]*scale,pt[2]*scale);
			break;
		case 4:
			return noise4d(pt[0]*scale,pt[1]*scale,pt[2]*scale,pt[3]*scale);
			break;
		default:
			break;
	}
	return 0.0;
}

double noise_fractal(int dim, double *pt, int octaves, double lacunarity, double persistence) {
	double noise,scale,lac;
	noise = 0.0;
	scale = 1.0;
	lac = 1.0;

	for (scale=1.0; octaves--;scale*=persistence) {
		noise += noise_simplex(dim, lac, pt) * scale;
		scale *= persistence;
		lac *= lacunarity;
	}
	return noise;
}

#ifdef MAIN

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

int 
main(int argc, char *argv[])
{
	double	r;
	double	base;
	double	s[2];

	int x,y;

	int octaves;
	int pos = 0, neg = 0;
	double max, min;
#ifdef INFINITY
	max = -INFINITY;
	min = +INFINITY;
#else
	max = DBL_MIN;
	min = DBL_MAX;
#endif

	octaves = argc > 1 ? atoi(argv[1]) : 1;
	base = 40.0;

	for (x=0;x<256;x++) {
		s[0] = x/base;
		for (y=0;y<256;y++) {
			s[1] = y/base;
			r = noise_fractal(2, s, octaves, 2.0, 0.7);
			r /= 1.5;

			/* grain */
			r = noise_simplex(2,0.20,s)*10;
			r = (r - (int)r);
			/* fine */
			r += noise_simplex(2, 75.0, s) * 0.20;
#if 1
			/* streak */
			r += noise2d(s[0]*0.45, s[1]*200.0) * 0.3;
#endif

			r /= 2.5;
			printf("%f\n", r);
			if (r > 1.0) {
				pos++;
			} else if (r < -1.0) {
				neg++;
			}
			if (r > max) max = r;
			if (r < min) min = r;

		}
	}
	fprintf(stderr, "max value: %f\n", max);
	fprintf(stderr, "min value: %f\n", min);
	if (pos) fprintf(stderr, "%d values > +1.0\n", pos);
	if (neg) fprintf(stderr, "%d values < -1.0\n", neg);
	
	return 0;

	switch (--argc) {
		case 1:
		r = noise1d(atof(argv[1]),0);
		break;

	case 2:
		r = noise2d(atof(argv[1]), atof(argv[2]));
		break;
	case 3:
		r = noise3d(atof(argv[1]), atof(argv[2]), atof(argv[3]));
		break;
	case 4:
		r = noise4d(atof(argv[1]), atof(argv[2]), atof(argv[3]), atof(argv[4]));
		break;
	default:
		fprintf(stderr, "Wrong number of arguments: %d, must be 2, 3, or 4\n", argc);
		exit(EXIT_FAILURE);
		break;
	};

	printf("%f\n", r);

	return 0;
}

#endif