#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>

/* The cells are stored in a 1-dimensional array of length n*n. 
   The square geometry has to be calculated on the fly:
   x = i%n
   y = i/n
   where i in [0..n*n[.
   There are also macros to find the index of any direct neighbour
   of a cell with index i. These macros incorporate periodic boundary
   conditions on all sides. (The macros take i in [0..n*n[ and n as
   arguments.)
*/
   
/* ------------------------------------------------------------ */

#define up( x, n ) ( ( ((x)/(n)+1  )%(n) )*(n) + (x)%(n) )
#define dn( x, n ) ( ( ((x)/(n)-1+n)%(n) )*(n) + (x)%(n) )
#define lf( x, n ) ( ( (x)/(n) )*(n) + ( (x)-1+(n) )%(n) )
#define rt( x, n ) ( ( (x)/(n) )*(n) + ( (x)+1     )%(n) )


/* ------------------------------------------------------------ */

/* Wrapper function - supply better random number generator if desired */
long rnd( long i ) { 
  return rand()%i;
}

/* Distributes m chips on a square of dimension (n x n) and returns
   the allocated matrix as pointer to long. */
long * init( long m, long n ) {
  long *data = (long *)malloc( n*n*sizeof( long ) );

  long i;
  long ctr = 0;

  if( m > n*n ) { 
    fprintf( stderr, "Too many chips (%ld) for field (%ld)\n", m, n*n );
    exit( EXIT_FAILURE );
  }

  if( !data ) {
    fprintf( stderr, "malloc( %ld ) failed\n", n*n );
    exit( EXIT_FAILURE );
  }

  for( i=0; i<n*n; i++ ) { data[i] = 0; } /* initialize matrix */

  while( ctr < m ) {     /* distribute chips... */
    i = rnd( n*n );
    if( data[i] == 0 ) { /* making sure there is at most one chip per cell */
      data[i] = 1;
      ctr += 1;
    }
  }

  return data;
}

/* ------------------------------------------------------------ */

/* Dumpt a graphic representation of the field to STDOUT.
   Arguments:
     t   : time step
     crr : index of the cell the walker is on
     occ : == 1 if walker is carrying a chip, 0 otherwise
   Output:
     - : empty cell
     # : chip
     X : walker carrying chip
     O : walker not carrying chip */
void dump( long t, long crr, int occ, long n, long *data ) {
  long i;

  printf( "\nt: %ld\n", t );
  for( i=0; i<n*n; i++ ) {
    if( i == crr ) { printf( "%c ", occ     ? 'X' : 'O' ); }
    else           { printf( "%c ", data[i] ? '#' : '-' ); }

    if( i%n == n-1 ) { printf( "\n" ); }
  }

  printf( "\n" );
}

/* Print data as text. */
void output( long t, long n, long *data ) {
  long i;

  printf( "# data t: %ld\n", t );
  for( i=0; i<n*n; i++ ) {
    if( data[i] ) { 
      printf( "%ld %ld %ld\n", i/n, i%n, data[i] ); /* x-pos y-pos data */
    }
  }
  printf( "\n\n" );
}

void print_usage( void ) {
  fprintf( stderr, "Usage: termites <chips m> <size n> <maxitr>\n" );
  exit( EXIT_SUCCESS );  
}

/* ------------------------------------------------------------ */

/* Helper function - used to find the correlation lenght using Newton */
double eqn( double xi, long n, double *c ) {
  long i;
  double z, res = 0.0;

  for( i=0; i<n; i++ ) {
    z = (double)i/xi;
    res += (exp(-z) - c[i])*exp(-z)*z;
  }

  return res;
}

/* Calculate the correlation function
   Arguments:
     t : timestep
     m : number of chips
     n : edge length */   
void corr( long t, long m, long n, long *data ) {
  int itr = 0;
  double x, err, de; 
  double g = (double)m/(double)(n*n); /* uniform density */

  long i, j;
  int x0, y0, x1, y1, dx, dy, r;

  long   *c  = (long *)malloc( n*sizeof(long) ); /* histogram */
  double *dc = (double *)malloc( n*sizeof(double) ); /* normalized corr fct */
  if( !c  ) { fprintf( stderr, "malloc 1 failed in corr()\n" ); exit( 1 ); }
  if( !dc ) { fprintf( stderr, "malloc 2 failed in corr()\n" ); exit( 1 ); }

  for( i=0; i<n; c[i++] = 0 ); /* initialize */


  /* Double-loop over all data cells to build up correlation histogram */
  for( i=0; i<n*n; i++ ) {
    if( data[i] == 0 ) { continue; } /* skip, if "origin" is empty */

    y0 = i/n; /* coordinates of "origin" */
    x0 = i%n;

    for( j=i; j<n*n; j++ ) {
      if( data[j] == 0 ) { continue; } /* skip, if "target" is empty */

      y1 = j/n; /* coordinates of "target" */
      x1 = j%n;

      dx = abs(x0-x1) < n/2 ? abs(x0-x1) : n - abs(x0-x1); /* periodic bdry */
      dy = abs(y0-y1) < n/2 ? abs(y0-y1) : n - abs(y0-y1);
      r = (int)sqrt( dx*dx + dy*dy ); /* distance, as integer */
      
      c[r] += 2; /* add 2 particles, for both pairs: (i,j) AND (j,i) */
    }
  }


  /* Normalize to obtain correlation function */
  for( i=0; i<n/2; i++ ) {
    double f = (double)c[i]/(double)m; /* normalize over number of "origins" */
    f /= M_PI*(double)(2*i+1); /* divide by area of ring: pi*((r+dr)^2 - r^2)*/
    dc[i] = (f-g)/(M_2_PI-g); /* subtrct uniform dens, rescale by val at orig*/
  }
  

  /* Newton-Raphson to fit the correlation length */
  x   = 1.0;
  err = eqn( x, n/2, dc );
  while( fabs(err) > 1.0e-4 && itr++ < 15 ) {
    de = ( eqn( x+1.0e-6, n/2, dc ) - err )/1.0e-6;
    x -= err/de;
    err = eqn( x, n/2, dc );
  }

  /* Finally: output */
  printf( "# corr t: %ld xi: %f\n", t, x );
  for( i=0; i<n/2; i++ ) {
    double f = (double)c[i]/( M_PI*(double)(m*(2*i+1) ) );
    printf( "%ld\t%ld %f\t%ld\t%f %f\n",t,i,(double)i/(double)n,c[i],f,dc[i]);
  }
  printf( "\n\n" );
}

/* ------------------------------------------------------------ */

int main( int argc, char *argv[] ) {
  int  occ;

  long m, n;
  long *data;
  long crr, tgt;
  long x, y;
  long t, maxitr;

  srand( time( NULL ) );
  if( argc != 4 ) { print_usage(); }

  m = (int)strtol( argv[1], NULL, 10 ); /* number of chips */
  n = (int)strtol( argv[2], NULL, 10 ); /* cells along one edge */
  maxitr = strtol( argv[3], NULL, 10 ); /* number of iteration steps */

  data = init( m, n );

  x = 0;    /* "true" position of walker */
  y = 0;    /*  ditto */
  occ = 0;  /* flag to indicate walker empty (0) or carrying chip (else) */
  crr = 0;  /* current position of walker as index in [0..n*n[ */

  for( t=0; t<maxitr; t++ ) {
    int r = rnd( 4 );

    switch( r ) {               /* select target cell to move to */
    case 0: tgt = dn( crr, n ); break; 
    case 1: tgt = up( crr, n ); break;
    case 2: tgt = lf( crr, n ); break;
    case 3: tgt = rt( crr, n ); break;
    default: 
      fprintf( stderr, "Illegal value %d in switch\n", r );
      exit( EXIT_FAILURE );
    }    

    if( data[tgt] == 0 ) { /* if target cell empty: move there */
      crr = tgt;
      switch( r ) {
      case 0: y--; break;
      case 1: y++; break;
      case 2: x--; break;
      case 3: x++; break;
      }
    } else {               /* if cell occupied... */
      if( occ == 1 ) {     /* ... and walker is empty */
	if( data[crr] == 0 ) { occ = 0; data[crr] = t; } /* pick up chip */
	else { fprintf( stderr, "Cannot happen!\n" ); exit( 1 ); }
      } else {             /* ... and walker is carrying */
	if( data[crr] == 0 ) { occ = 1; data[tgt] = 0; } /* drop of chip */
	else { /* no-op */ }
      }
    }

    // printf( "%ld %ld %ld %ld\n", t, x*x+y*y, x, y );
    if( t%(maxitr/10L) == 0 ) { /* generate 10 equally spaced data dumps */
      // dump( t, crr, occ, n, data );

      output( t, n, data );
      corr( t, m, n, data );
    }
  }

  return EXIT_SUCCESS;
}