typedef double number; // used in macros to indicate "any numeric type"
#define pi 3.14159265358979
#undef HUGE
#define HUGE 1e30f
#define nodata HUGE

macro number max (number a, number b) { return a > b ? a : b; }
macro number min (number a, number b) { return a < b ? a : b; }
macro number sq (number x) { return x*x; }
macro number cube (number x) { return x*x*x; }
macro int sign (number x) { return (int)(x > 0 ? 1 : -1); }
macro int sign2 (number x) { return (int)(x > 0 ? 1 : x < 0 ? -1 : 0); }
macro number clamp (number x, number a, number b) {
  return x < a ? a : x > b ? b : x;
}

#define swap(type,a,b) do { type _tmp_ = a; a = b; b = _tmp_; } while(false)

#include "grid/config.h"

static inline double noise() { return 1. - 2.*rand()/(double)RAND_MAX; }

// the grid
typedef struct {
  long n;       // number of (leaf) cells for this process
  long tn;      // number of (leaf) cells for all processes
  int depth;    // the depth for this process
  int maxdepth; // the maximum depth for all processes
} Grid;
Grid * grid = NULL;
// coordinates of the lower-left corner of the box
double X0 = 0., Y0 = 0., Z0 = 0.;
// size of the box
double L0 = 1. [1];
// number of grid points
#if dimension <= 2
int N = 64;
#else
int N = 16;
#endif

typedef struct { int i; } scalar;

typedef struct {
  scalar x;
#if dimension > 1
  scalar y;
#endif
#if dimension > 2
  scalar z;
#endif
} vector;

typedef struct {
  scalar * x;
#if dimension > 1
  scalar * y;
#endif
#if dimension > 2
  scalar * z;
#endif
} vectorl;

typedef struct {
  vector x;
#if dimension > 1
  vector y;
#endif
#if dimension > 2
  vector z;
#endif
} tensor;

struct { int x, y, z; } Period = {false, false, false};

typedef struct {
  double x, y, z;
} coord;

OMP(omp declare reduction (+ : coord :
			   omp_out.x += omp_in.x,
			   omp_out.y += omp_in.y,
			   omp_out.z += omp_in.z))

#if dimension == 1
# define norm(v) fabs(v.x[])
# define dv() (Delta*cm[])
#elif dimension == 2
# define norm(v) (sqrt(sq(v.x[]) + sq(v.y[])))
# define dv() (sq(Delta)*cm[])
#else // dimension == 3
# define norm(v) (sqrt(sq(v.x[]) + sq(v.y[]) + sq(v.z[])))
# define dv() (cube(Delta)*cm[])
#endif

void normalize (coord * n)
{
  double norm = 0.;
  foreach_dimension()
    norm += sq(n->x);
  norm = sqrt(norm);
  foreach_dimension()
    n->x /= norm;
}

void origin (double x = 0., double y = 0., double z = 0.) {
  X0 = x; Y0 = y; Z0 = z;
}

void size (double L) {
  L0 = L;
}

double zero (double s0, double s1, double s2) { return 0.; }

// boundary conditions for each direction/variable

#if dimension == 1
  enum { right, left };
#elif dimension == 2
  enum { right, left, top, bottom };
#else
  enum { right, left, top, bottom, front, back };
#endif
int nboundary = 2*dimension;

#define none -1

double  * _constant = NULL;
size_t datasize = 0;
typedef struct _Point Point;

#include "grid/boundaries.h"

// attributes for each scalar

typedef struct {
  int x;
#if dimension > 1
  int y;
#endif
#if dimension > 2
  int z;
#endif
} ivec;
typedef double (* BoundaryFunc) (Point, Point, scalar, bool *);
typedef struct {
  BoundaryFunc * boundary;
  BoundaryFunc * boundary_homogeneous;
  double (* gradient)              (double, double, double);
  void   (* delete)                (scalar);
  char * name;
  ivec d; // staggering
  vector v;
  int face;
  bool   nodump, freed;
  int    block;
  scalar * depends; // boundary conditions depend on other fields
} _Attributes;

static _Attributes * _attribute = NULL;

#define foreach_block() // treatment of block data is disabled by default
#define foreach_blockf(s)

#if dimension == 1
ivec Dimensions = {1};
#elif dimension == 2
ivec Dimensions = {1,1};
#elif dimension == 3
ivec Dimensions = {1,1,1};
#endif

// lists

int list_len (scalar * list)
{
  if (!list) return 0;
  int ns = 0;
  for (scalar s in list) ns++;
  return ns;
}

scalar * list_append (scalar * list, scalar s)
{
  int len = list_len (list);
  qrealloc (list, len + 2, scalar);
  list[len] = s;
  list[len + 1].i = -1;
  return list;
}

scalar * list_prepend (scalar * list, scalar s)
{
  int len = list_len (list);
  qrealloc (list, len + 2, scalar);
  for (int i = len; i >= 1; i--)
    list[i] = list[i-1];
  list[0] = s;
  list[len + 1].i = -1;
  return list;
}

scalar * list_add (scalar * list, scalar s)
{
  for (scalar t in list)
    if (t.i == s.i)
      return list;
  return list_append (list, s);
}

int list_lookup (scalar * l, scalar s)
{
  if (l != NULL)
    for (scalar s1 in l)
      if (s1.i == s.i)
	return true;
  return false;
}

scalar * list_copy (scalar * l)
{
  scalar * list = NULL;
  if (l != NULL)
    for (scalar s in l)
      list = list_append (list, s);
  return list;
}

scalar * list_concat (scalar * l1, scalar * l2)
{
  scalar * l3 = list_copy (l1);
  for (scalar s in l2)
    l3 = list_append (l3, s);
  return l3;
}

void list_print (scalar * l, FILE * fp)
{
  int i = 0;
  for (scalar s in l)
    fprintf (fp, "%s%s", i++ == 0 ? "{" : ",", s.name);
  fputs (i > 0 ? "}\n" : "{}\n", fp);
}

int vectors_len (vector * list)
{
  if (!list) return 0;
  int nv = 0;
  for (vector v in list) nv++;
  return nv;
}

vector * vectors_append (vector * list, vector v)
{
  int len = vectors_len (list);
  qrealloc (list, len + 2, vector);
  list[len] = v;
  list[len + 1] = (vector){{-1}};
  return list;
}

vector * vectors_add (vector * list, vector v)
{
  for (vector w in list) {
    bool id = true;
    foreach_dimension()
      if (w.x.i != v.x.i)
	id = false;
    if (id)
      return list;
  }
  return vectors_append (list, v);
}

vector * vectors_copy (vector * l)
{
  vector * list = NULL;
  if (l != NULL)
    for (vector v in l)
      list = vectors_append (list, v);
  return list;
}

vector * vectors_from_scalars (scalar * s)
{
  vector * list = NULL;
  while (s->i >= 0) {
    vector v;
    foreach_dimension() {
      assert (s->i >= 0);
      v.x = *s++;
    }
    list = vectors_append (list, v);
  }
  return list;
}

int tensors_len (tensor * list)
{
  if (!list) return 0;
  int nt = 0;
  for (tensor t in list) nt++;
  return nt;
}

tensor * tensors_append (tensor * list, tensor t)
{
  int len = tensors_len (list);
  qrealloc (list, len + 2, tensor);
  list[len] = t;
  list[len + 1] = (tensor){{{-1}}};
  return list;
}

tensor * tensors_from_vectors (vector * v)
{
  tensor * list = NULL;
  while (v->x.i >= 0) {
    tensor t;
    foreach_dimension() {
      assert (v->x.i >= 0);
      t.x = *v++;
    }
    list = tensors_append (list, t);
  }
  return list;
}

static inline bool is_vertex_scalar (scalar s)
{
  foreach_dimension()
    if (s.d.x != -1)
      return false;
  return true;
}

scalar * all = NULL; // all the fields
scalar * baseblock = NULL; // base block fields

// basic methods

scalar (* init_scalar)        (scalar, const char *);
scalar (* init_vertex_scalar) (scalar, const char *);
vector (* init_vector)        (vector, const char *);
vector (* init_face_vector)   (vector, const char *);
tensor (* init_tensor)        (tensor, const char *);
void   (* scalar_clone)       (scalar, scalar);

#define vector(x) (*((vector *)&(x)))

// timers

#if _MPI
static double mpi_time = 0.;
#endif

typedef struct {
  clock_t c;
  struct timeval tv;
  double tm;
} timer;

timer timer_start (void)
{
  timer t;
  t.c = clock();
  gettimeofday (&t.tv, NULL);
#if _MPI
  t.tm = mpi_time;
#endif
  return t;
}

double timer_elapsed (timer t)
{
  struct timeval tvend;
  gettimeofday (&tvend, NULL);
  return ((tvend.tv_sec - t.tv.tv_sec) + 
	  (tvend.tv_usec - t.tv.tv_usec)/1e6);
}

// Constant fields

const face vector zerof[] = {0.,0.,0.};
const face vector unityf[] = {1.,1.,1.};
const scalar unity[] = 1.;
const scalar zeroc[] = 0.;

// Metric

(const) face vector fm[] = {1.[0],1.[0],1.[0]};
(const) scalar cm[] = 1.[0];

// Embedded boundaries
// these macros are overloaded in embed.h

#define SEPS 0.
#define face_gradient_x(a,i) ((a[i] - a[i-1])/Delta)
#define face_gradient_y(a,i) ((a[0,i] - a[0,i-1])/Delta)
#define face_gradient_z(a,i) ((a[0,0,i] - a[0,0,i-1])/Delta)
#define face_value(a,i)      ((a[i] + a[i-1])/2.)
#define center_gradient(a)   ((a[1] - a[-1])/(2.*Delta))

// matrices

void * matrix_new (int n, int p, size_t size)
{
  void ** m = qmalloc (n, void *);
  char * a = qmalloc (n*p*size, char);
  for (int i = 0; i < n; i++)
    m[i] = a + i*p*size;
  return m;
}

double matrix_inverse (double ** m, int n, double pivmin)
{
  int indxc[n], indxr[n], ipiv[n];
  int i, icol = 0, irow = 0, j, k, l, ll;
  double big, dum, pivinv, minpiv = HUGE;

  for (j = 0; j < n; j++)
    ipiv[j] = -1;

  for (i = 0; i < n; i++) {
    big = 0.0;
    for (j = 0; j < n; j++)
      if (ipiv[j] != 0)
	for (k = 0; k < n; k++) {
	  if (ipiv[k] == -1) {
	    if (fabs (m[j][k]) >= big) {
	      big = fabs (m[j][k]);
	      irow = j;
	      icol = k;
	    }
	  }
	}
    ipiv[icol]++;
    if (irow != icol)
      for (l = 0; l < n; l++) 
	swap (double, m[irow][l], m[icol][l]);
    indxr[i] = irow;
    indxc[i] = icol;
    if (fabs (m[icol][icol]) <= pivmin)
      return 0.;
    if (fabs (m[icol][icol]) < minpiv)
      minpiv = fabs (m[icol][icol]);
    pivinv = 1.0/m[icol][icol];
    m[icol][icol] = 1.0;
    for (l = 0; l < n; l++) m[icol][l] *= pivinv;
    for (ll = 0; ll < n; ll++)
      if (ll != icol) {
	dum = m[ll][icol];
	m[ll][icol] = 0.0;
	for (l = 0; l < n; l++)
	  m[ll][l] -= m[icol][l]*dum;
      }
  }
  for (l = n - 1; l >= 0; l--) {
    if (indxr[l] != indxc[l])
      for (k = 0; k < n; k++)
	swap (double, m[k][indxr[l]], m[k][indxc[l]]);
  }
  return minpiv;
}

void matrix_free (void * m)
{
  free (((void **) m)[0]);
  free (m);
}

// Solver cleanup

typedef void (* free_solver_func) (void);

static Array * free_solver_funcs = NULL;

void free_solver_func_add (free_solver_func func)
{
  if (!free_solver_funcs)
    free_solver_funcs = array_new();
  array_append (free_solver_funcs, &func, sizeof(free_solver_func));
}

// Default objects to display

static char * display_defaults = NULL;

static void free_display_defaults() {
  free (display_defaults);
}

void display (const char * commands, bool overwrite = false)
{
  if (display_defaults == NULL)
    free_solver_func_add (free_display_defaults);
  if (overwrite) {
    free (display_defaults);
    display_defaults = malloc (strlen(commands) + 2);
    strcpy (display_defaults, "@");
    strcat (display_defaults, commands);
  }
  else {
    if (!display_defaults)
      display_defaults = strdup ("@");
    display_defaults =
      realloc (display_defaults,
	       strlen(display_defaults) + strlen(commands) + 1);
    strcat (display_defaults, commands);
  }
}

#define display_control(val, ...)

typedef struct {
  double x;
#if dimension > 1
  double y;
#endif
#if dimension > 2
  double z;
#endif
} _coord;

// Types and macros for compatibility with GLSL

typedef struct {
  float r, g, b, a;
} vec4;

#if dimension == 1
# define avector(x, ...)    {x}
#elif dimension == 2
# define avector(x, y, ...) {x, y}
#else // dimension == 3
# define avector(x, y, z)   {x, y, z}
#endif

typedef struct {
  coord x, y, z;
} mat3;

OMP(omp declare reduction (+ : mat3 :
			   omp_out.x.x += omp_in.x.x,
			   omp_out.x.y += omp_in.x.y,
			   omp_out.x.z += omp_in.x.z,
			   omp_out.y.x += omp_in.y.x,
			   omp_out.y.y += omp_in.y.y,
			   omp_out.y.z += omp_in.y.z,
			   omp_out.z.x += omp_in.z.x,
			   omp_out.z.y += omp_in.z.y,
			   omp_out.z.z += omp_in.z.z
			   ))

typedef struct {
  uint32_t s;
} Adler32Hash;

static
inline void a32_hash_init (Adler32Hash * hash)
{
  hash->s = 0;
}

static
inline void a32_hash_add (Adler32Hash * hash, const void * data, size_t size)
{
  const uint8_t * buffer = (const uint8_t*) data;
  for (size_t n = 0; n < size; n++, buffer++)
    hash->s = *buffer + (hash->s << 6) + (hash->s << 16) - hash->s;
}

static
inline uint32_t a32_hash (const Adler32Hash * hash)
{
  return hash->s;
}