static const double noise_d3x[] = { 1.0D, -1.0D, 1.0D, -1.0D, 1.0D, -1.0D, 1.0D, -1.0D, 0.0D, 0.0D, 0.0D, 0.0D, 1.0D, 0.0D, -1.0D, 0.0D };
static const double noise_d3y[] = { 1.0D, 1.0D, -1.0D, -1.0D, 0.0D, 0.0D, 0.0D, 0.0D, 1.0D, -1.0D, 1.0D, -1.0D, 1.0D, -1.0D, 1.0D, -1.0D };
static const double noise_d3z[] = { 0.0D, 0.0D, 0.0D, 0.0D, 1.0D, 1.0D, -1.0D, -1.0D, 1.0D, 1.0D, -1.0D, -1.0D, 0.0D, 1.0D, 0.0D, -1.0D };
static const double noise_d2x[] = { 1.0D, -1.0D, 1.0D, -1.0D, 1.0D, -1.0D, 1.0D, -1.0D, 0.0D, 0.0D, 0.0D, 0.0D, 1.0D, 0.0D, -1.0D, 0.0D };
static const double noise_d2z[] = { 0.0D, 0.0D, 0.0D, 0.0D, 1.0D, 1.0D, -1.0D, -1.0D, 1.0D, 1.0D, -1.0D, -1.0D, 0.0D, 1.0D, 0.0D, -1.0D };

typedef struct {
	int permutations[512];
	double xcoord;
	double ycoord;
	double zcoord;
} octave_t;

typedef struct {
	octave_t *octaves;
	uint n_octaves;
} noisegen_t;

void noise_oct_init(octave_t *oct, ulong *rand) {
	memset(oct->permutations, 0, 512 * sizeof(int));
	oct->xcoord = rng_double(rand) * 256.0D;
	oct->ycoord = rng_double(rand) * 256.0D;
	oct->zcoord = rng_double(rand) * 256.0D;

	for(int n = 0; n < 256; oct->permutations[n] = n++) {
		;
	}

	for(int n = 0; n < 256; ++n) {
		int pos = rng_zrange(rand, 256 - n) + n;
		int value = oct->permutations[n];
		oct->permutations[n] = oct->permutations[pos];
		oct->permutations[pos] = value;
		oct->permutations[n + 256] = oct->permutations[n];
	}
}

void noise_gen_init(noisegen_t *gen, ulong *rand, int octs) {
	gen->n_octaves = octs;
	gen->octaves = mem_alloc(octs * sizeof(octave_t), MEM_MISC);

	for(int n = 0; n < octs; ++n) {
		noise_oct_init(&gen->octaves[n], rand);
	}
}

void noise_gen_free(noisegen_t *gen) {
	if(gen->octaves)
		mem_free(gen->octaves);
	gen->octaves = NULL;
	gen->n_octaves = 0;
}

double noise_lerp(double d, double n, double m) {
	return n + d * (m - n);
}

double noise_grad2d(int n, double xm, double zm) {
	int i = n & 15;
	return noise_d2x[i] * xm + noise_d2z[i] * zm;
}

double noise_grad(int n, double xm, double ym, double zm) {
	int i = n & 15;
	return noise_d3x[i] * xm + noise_d3y[i] * ym + noise_d3z[i] * zm;
}

long noise_floor_double(double value) {
	long i = (long)value;
	return value < (double)i ? i - 1L : i;
}

void noise_oct_gen(octave_t *oct, double *noise, double xoff, double yoff, double zoff, int xsize, int ysize, int zsize, 
		double xscale, double yscale, double zscale, double scale) {
	int *permutations = oct->permutations;
	double xcoord = oct->xcoord;
	double ycoord = oct->ycoord;
	double zcoord = oct->zcoord;
	if(ysize == 1) {
		int i5 = 0;
		int j5 = 0;
		int j = 0;
		int k5 = 0;
		double d14 = 0.0D;
		double d15 = 0.0D;
		int l5 = 0;
		double d16 = 1.0D / scale;

		for(int j2 = 0; j2 < xsize; ++j2) {
			double d17 = xoff + (double)j2 * xscale + xcoord;
			int i6 = (int)d17;

			if(d17 < (double)i6) {
				--i6;
			}

			int k2 = i6 & 255;
			d17 = d17 - (double)i6;
			double d18 = d17 * d17 * d17 * (d17 * (d17 * 6.0D - 15.0D) + 10.0D);

			for(int j6 = 0; j6 < zsize; ++j6) {
				double d19 = zoff + (double)j6 * zscale + zcoord;
				int k6 = (int)d19;

				if(d19 < (double)k6) {
					--k6;
				}

				int l6 = k6 & 255;
				d19 = d19 - (double)k6;
				double d20 = d19 * d19 * d19 * (d19 * (d19 * 6.0D - 15.0D) + 10.0D);
				i5 = permutations[k2] + 0;
				j5 = permutations[i5] + l6;
				j = permutations[k2 + 1] + 0;
				k5 = permutations[j] + l6;
				d14 = noise_lerp(d18, noise_grad2d(permutations[j5], d17, d19), noise_grad(permutations[k5], d17 - 1.0D, 0.0D, d19));
				d15 = noise_lerp(d18, noise_grad(permutations[j5 + 1], d17, 0.0D, d19 - 1.0D),
						noise_grad(permutations[k5 + 1], d17 - 1.0D, 0.0D, d19 - 1.0D));
				double d21 = noise_lerp(d20, d14, d15);
				int i7 = l5++;
				noise[i7] += d21 * d16;
			}
		}
	}
	else {
		int i = 0;
		double d0 = 1.0D / scale;
		int k = -1;
		int l = 0;
		int i1 = 0;
		int j1 = 0;
		int k1 = 0;
		int l1 = 0;
		int i2 = 0;
		double d1 = 0.0D;
		double d2 = 0.0D;
		double d3 = 0.0D;
		double d4 = 0.0D;

		for(int l2 = 0; l2 < xsize; ++l2) {
			double d5 = xoff + (double)l2 * xscale + xcoord;
			int i3 = (int)d5;

			if(d5 < (double)i3) {
				--i3;
			}

			int j3 = i3 & 255;
			d5 = d5 - (double)i3;
			double d6 = d5 * d5 * d5 * (d5 * (d5 * 6.0D - 15.0D) + 10.0D);

			for(int k3 = 0; k3 < zsize; ++k3) {
				double d7 = zoff + (double)k3 * zscale + zcoord;
				int l3 = (int)d7;

				if(d7 < (double)l3) {
					--l3;
				}

				int i4 = l3 & 255;
				d7 = d7 - (double)l3;
				double d8 = d7 * d7 * d7 * (d7 * (d7 * 6.0D - 15.0D) + 10.0D);

				for(int j4 = 0; j4 < ysize; ++j4) {
					double d9 = yoff + (double)j4 * yscale + ycoord;
					int k4 = (int)d9;

					if(d9 < (double)k4) {
						--k4;
					}

					int l4 = k4 & 255;
					d9 = d9 - (double)k4;
					double d10 = d9 * d9 * d9 * (d9 * (d9 * 6.0D - 15.0D) + 10.0D);

					if(j4 == 0 || l4 != k) {
						k = l4;
						l = permutations[j3] + l4;
						i1 = permutations[l] + i4;
						j1 = permutations[l + 1] + i4;
						k1 = permutations[j3 + 1] + l4;
						l1 = permutations[k1] + i4;
						i2 = permutations[k1 + 1] + i4;
						d1 = noise_lerp(d6, noise_grad(permutations[i1], d5, d9, d7), noise_grad(permutations[l1], d5 - 1.0D, d9, d7));
						d2 = noise_lerp(d6, noise_grad(permutations[j1], d5, d9 - 1.0D, d7),
								noise_grad(permutations[i2], d5 - 1.0D, d9 - 1.0D, d7));
						d3 = noise_lerp(d6, noise_grad(permutations[i1 + 1], d5, d9, d7 - 1.0D),
								noise_grad(permutations[l1 + 1], d5 - 1.0D, d9, d7 - 1.0D));
						d4 = noise_lerp(d6, noise_grad(permutations[j1 + 1], d5, d9 - 1.0D, d7 - 1.0D),
								noise_grad(permutations[i2 + 1], d5 - 1.0D, d9 - 1.0D, d7 - 1.0D));
					}

					double d11 = noise_lerp(d10, d1, d2);
					double d12 = noise_lerp(d10, d3, d4);
					double d13 = noise_lerp(d8, d11, d12);
					int j7 = i++;
					noise[j7] += d13 * d0;
				}
			}
		}
	}
}

void noise_gen(noisegen_t *gen, double *noise, double xoff, double yoff, double zoff, int xsize, int ysize, int zsize, 
		double xscale, double yscale, double zscale) {
	memset(noise, 0, sizeof(double) * xsize * ysize * zsize);
	uint octaves = gen->n_octaves;
	double factor = 1.0D;
	for(int n = 0; n < octaves; ++n) {
		double xo = xoff * factor * xscale;
		double yo = yoff * factor * yscale;
		double zo = zoff * factor * zscale;
		long xn = noise_floor_double(xo);
		long zn = noise_floor_double(zo);
		xo = xo - (double)xn;
		zo = zo - (double)zn;
		xn = xn % 16777216L;
		zn = zn % 16777216L;
		xo = xo + (double)xn;
		zo = zo + (double)zn;
		noise_oct_gen(&gen->octaves[n], noise, xo, yo, zo, xsize, ysize, zsize, xscale * factor, yscale * factor, zscale * factor, factor);
		factor /= 2.0D;
	}
}