/*
 * Four dimensional noise routines using compressed space.
 * Written by Geoff Wyvill, December 1998
 * Copyright 1998 Geoff Wyvill <geoff@cs.otago.ac.nz>
 * This source may be freely distributed for educational purposes
 * If you have a commercial application, contact the author.
 * The software that generates the tables is not currently available for
 * distribution. There is a known optimisation that will make this code
 * four to five times faster but the implementation is tedious.
 * The improved version will probably be available in October 1999.
 */
 
#include <math.h>

// Congruential sequence constants. Suggested by Neave and Knuth.
#define RANDOM_BEE 3125
#define RANDOM_CEE 49
#define WIDTH 1.4
#define RANGE 4294967296.0

typedef double vecfour[4];
typedef double vector[3];

typedef struct { vecfour point;
                 int mult, add;
               } pivot;

#include "tables.c"
/*
 *  Generated file contains the 4D matrix of data points with
 *  relative co-ordinates and random generation constants:
 *
 *   cubelength = 125000000, layerlength = 250000, rowlength = 500;
 *   static pivot access[] = { generated data };
 */

unsigned int any_random (unsigned int k)
{
unsigned int b = RANDOM_BEE;
unsigned int c = RANDOM_CEE;
unsigned int result = 1; /* have to start somewhere */

for (;k > 0;k>>=1)
	{
	if (k & 1) result = result * b + c;
	c += b * c;
	b *= b;
	}
return result;
}

double four_noise(vector sample, double time)
{ unsigned int seed;
  int ix, iy, iz, iw, index;
  double C, sum = 0.0, weight = 0.0, x, y, z, w, rsquare;
// unskew
  w = time;
  x = sample[0]-w;
  y = sample[1]-w;
  z = sample[2]-w;
  w *= 2.0;
// find base
  ix = floor(x) - 2;
  iy = floor(y) - 2;
  iz = floor(z) - 2;
  iw = floor(w) - 2;
// find seed at base  
  seed = any_random(cubelength*iw + layerlength*iz + rowlength*iy + ix);
// unskewed offsets
  x = x - ix;
  y = y - iy;
  z = z - iz;
  w = w - iw;
//skew
  w *= 0.5;
  x += w;
  y += w;
  z += w;
  for (index = 0; index < TABLESIZE; index++)
  { rsquare = (access[index].point[0]-x)*(access[index].point[0]-x) +
              (access[index].point[1]-y)*(access[index].point[1]-y) +
              (access[index].point[2]-z)*(access[index].point[2]-z) +
              (access[index].point[3]-w)*(access[index].point[3]-w);
//printf("index %d rsquare %f, %d\n", index, rsquare, TABLESIZE);
//printf("at: %f %f %f %f\n", x, y, z, w);
    if (rsquare < WIDTH * WIDTH)
    { C = 1.0 - rsquare/(WIDTH * WIDTH);
      C *= C * C;
      sum += C * (seed * access[index].mult + access[index].add)/RANGE;
      weight += C;
//printf("sum %f, weight %f\n", sum, weight);
    }
  }
//printf("sum %f, weight %f\n", sum, weight);
  return sum/weight;
//printf("returning\n");
}