#include <stdio.h>        /*@\colabel{first_line}\hfill@*//* For standard I/O. */
#include <stdlib.h>       /*@\hfill@*//* For malloc and exit. */
#include <time.h>         /*@\hfill@*//* For time and ctime. */
#include <limits.h>       /*@\hfill@*//* For ULLONG_MAX, if it exists. */
#ifndef ULLONG_MAX
#define ULLONG_MAX (~0ULL)
#endif
#include <fpu_control.h>  /*@\colabel{fpu_control}\hfill@*//* For setting the precision. */
/* More convenient interface to set the precision: */
#define FPSETPREC(a) {fpu_control_t flags; _FPU_GETCW(flags); flags &= ~(_FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE); flags |= (a); _FPU_SETCW(flags);}/*@\colabel{FPSETPREC}@*/
#include <signal.h>       /*@\hfill@*//* For signal handling. */
#include <sys/types.h>    /*@\hfill@*//* For determining the pid of the program. */
#include <unistd.h>       /*@\hfill@*//* For pid and hostname. */
#define LENGTH    (64)        /*@\colabel{length}\hfill@*//* System size. */
#define OBC                   /*@\hfill@*//* Boundary condition. */ 
#define NUM_BOXES (LENGTH*32) /*@\colabel{num_boxes}\hfill@*//* Number of boxes. */
#define ALPHA  (0.07)         /*@\hfill@*//* Level of conservation. */
#define NUM_CHUNKS (5)        /*@\colabel{def_NUM_CHUNKS}\hfill@*//* Number of chunks. */
#define ITERATIONS_PER_CHUNK  (1000000LL) /*@\hfill@*//* ... */
#define SEED (1L)             /*@\colabel{seed_macro}\hfill@*//* Seed of the random number generator. */
/*@\hspace*{\codeheadercommentindent}@*//* 64 bit long long */
#if (1)                  /*@\colabel{comment_in}\colabel{force_ll}\hfill@*//* For quick commenting out. */
typedef unsigned long long int force_type;   /*@\colabel{lng_force_type}\hfill@*//* Type of force. */
/* Name of the force type for result file: */
#define FORCE_TYPE_NAME "unsigned long long"
#define F_SLIP ((force_type)1<<40)           /*@\colabel{lng_f_slip}\hfill@*//* Force threshold. */
/* Minimal force at which slipping takes place: */
#define INIT_OFFSET (ULLONG_MAX-64*F_SLIP)   /*@\colabel{lng_init_offset}@*/
/* Criterion to trigger readjustment of forces: */
#define READJUST_CRITERION(a) ((a)<F_SLIP)   /*@\colabel{lng_READJUST_CRITERION}@*/
#define FORCE_OUT_FMT "%llu"                 /*@\hfill@*//* Output format for forces. */
#define PRECISION _FPU_EXTENDED              /*@\hfill@*//* Precision setting. */
/* Calculate charge in accordance with force_type: */
#define CHARGE(a) (force_type)((long double)ALPHA*(long double)(a)) 
#endif

/*@\hspace*{\codeheadercommentindent}@*//* 64 bit double */
#if (0)                  /*@\colabel{comment_out}\colabel{force_dbl}\hfill@*//* For quick commenting out. */
typedef double force_type;  /*@\colabel{dbl_force_type}@*/
#define FORCE_TYPE_NAME "double"
#define F_SLIP ((force_type)1.0)          /*@\colabel{dbl_f_slip}@*/
#define INIT_OFFSET ((force_type)(1<<24)) /*@\colabel{dbl_init_offset}@*/
#define READJUST_CRITERION(a) ((a)<F_SLIP)/*@\colabel{dbl_READJUST_CRITERION}@*/
#define FORCE_OUT_FMT "%g"                 
#define PRECISION _FPU_DOUBLE             /*@\colabel{dbl_prec_clash}\hfill@*/
/* Note: Precision clashes with MOMENT_TYPE below (long double). */
#define CHARGE(a) ((force_type)ALPHA*(a))  
#endif

force_type f_threshold; /*@\colabel{f_threshold}@*/
struct lattice {
  force_type f;              /*@\hfill@*//* Local force. */
  force_type charge;         /*@\hfill@*//* Force to be charged. */
#define NUM_NEIGHBOURS (4)
  struct lattice *nb[NUM_NEIGHBOURS];     /*@\colabel{nb_pointer}\hfill@*//* Pointers to neighbours. */
  struct lattice  **bx_left; 
  /* ... the address of (neighbour's) pointer to me. */
  struct lattice  *bx_right; /*@\hfill@*//* Box: pointer to the next site. */
} *site;
#define NUM_SITES (LENGTH*LENGTH)
#define COO2INDEX(x,y) ((x)*LENGTH+(y))/*@\colabel{coo2index}@*/
/*@\hspace*{\codeheadercommentindent}@*//* Stacks */
int stackheightCurr, stackheightNext;
struct lattice **stackCurr, **stackNext;

#define PUSH(a) stackNext[stackheightNext++]=(a)
#define POP(a)  (a)=stackNext[--stackheightNext]
#define MOMENT_TYPE long double   /*@\colabel{moment_type}@*/
#define MOMENT_OUT_FMT "%10.20Lg"

MOMENT_TYPE mom_pow;
int         mom_cnt;

#define MOMENTS_DECL(n,m) MOMENT_TYPE mom_ ##n[(m) + 1]; const int mom_max_ ##n = m;
#define MOMENTS_INIT(n) {for (mom_cnt=0; mom_cnt<=mom_max_ ##n; mom_cnt++) mom_ ##n [mom_cnt]=0;}
#define MOMENTS(n,x) {mom_ ##n [0]++; for (mom_pow=1, mom_cnt=1; mom_cnt<=mom_max_ ##n; mom_cnt++) {mom_pow*=((MOMENT_TYPE)(x)); mom_ ##n[mom_cnt]+=mom_pow;}}
#define MOMENTS_OUT(n)  {force_type chunk_time; \
  chunk_time = total_time_readjusted - chunk_start_readjusted; \
  chunk_time += (total_time - chunk_start); \
  printf("#M_%s %i %i %g %lli " FORCE_OUT_FMT " " FORCE_OUT_FMT " " FORCE_OUT_FMT " %i " MOMENT_OUT_FMT " ", #n, chunk, LENGTH, ALPHA, ITERATIONS_PER_CHUNK, chunk_time, total_time + total_time_readjusted, f_threshold, mom_max_ ##n, mom_ ##n [0]); \
  for (mom_cnt=0; mom_cnt<=mom_max_ ##n; mom_cnt++) \
  printf(" " MOMENT_OUT_FMT, mom_ ##n [mom_cnt]/((mom_ ##n [0]) ? mom_ ##n [0] : -1)); fputc('\n', stdout);}

MOMENTS_DECL(siz, 5); /*@\hfill@*//* Avalanche size. */
MOMENTS_DECL(dur, 5); /*@\hfill@*//* Avalanche duration. */
MOMENTS_DECL(deg, 2); /*@\hfill@*//* Degeneracy of the maximum. */
/*@\hspace*{\codeheadercommentindent}@*//* Signalling */
void signal_handler(int signo);
int global_stop_flg=0;
#define BOX_POS(f) ((int)(((long long unsigned)((double)(NUM_BOXES-1)*(double)(f)/(double)F_SLIP))%(long long unsigned)NUM_BOXES))  /*@\colabel{BOX_POS}@*/

struct lattice *box[NUM_BOXES]; /*@\colabel{def_box}@*/

int box_insert(struct lattice *a);/*@\colabel{box_insert}@*/
int box_delete(struct lattice *a);
int box_getmaxima(struct lattice **list, int *num);/*@\colabel{box_getmaxima}@*/
int box_insert(struct lattice *a)
{
  int b;

  b=BOX_POS(a->f); /*@\colabel{get_box_pos}@*/
  if ((a->bx_right=box[b])!=NULL) box[b]->bx_left=&(a->bx_right);/*@\colabel{connect_to_a}@*/
  box[b]=a;
  a->bx_left=&box[b];
return(0);     /*@\hfill@*//* Success. */
}
int box_delete(struct lattice *a) /*@\colabel{box_delete_start}@*/
{
  if ((*(a->bx_left)=a->bx_right)!=NULL)
    a->bx_right->bx_left=a->bx_left; /*@\colabel{delete_back_ref}@*/
return(0);     /*@\hfill@*//* Success. */
}
int box_getmaxima(struct lattice **list, int *num)
{
  int b;
  struct lattice *p;

  for (b=BOX_POS(f_threshold); ; b=(b-1+NUM_BOXES)%NUM_BOXES) {
    if (box[b]==NULL) continue;
    list[0]=box[b]; /*@\colabel{first_candidate}@*/
    *num=1;
    for (p=box[b]->bx_right; (p!=NULL); p=p->bx_right) {
      if (p->f>list[0]->f) list[0]=p, *num=1;      /*@\hfill@*//* New max found. */
      else if (p->f==list[0]->f) list[(*num)++]=p; /*@\hfill@*//* Another, same...*/
                                                   /*@\hfill@*//* ... max found. */
    }
    return(0);
  }
  /*@\colabel{example_assertion}\hfill@*//* Example for an assertion: Code cannot reach this point. */
return(-1);
}
int main(int argc, char *argv[])
{
long long int it;                               /*@\hfill@*//* Iteration counter. */
int chunk;                                      /*@\hfill@*//* Chunk counter. */
int i, j;                                       /*@\hfill@*//* Auxiliary variables. */
struct lattice *nb;                             /*@\hfill@*//* Pointer to neighbour. */
force_type d, charge;                           /*@\hfill@*//* Derivation of charges. */
int s, t;                                       /*@\hfill@*//* Avalanche size and duration. */
force_type f_threshold_start;                   /*@\hfill@*//* Force threshold at start of time. */
force_type total_time;                          /*@\hfill@*//* Total time counter. */
force_type total_time_readjusted;               /*@\hfill@*//* Total at readjustment. */
force_type chunk_start, chunk_start_readjusted; /*@\hfill@*//* Chunk start time. */
sigset_t sigmask;                               /*@\hfill@*//* Signalling. */


{ char hostname[256]; /*@\hfill@*//* HOST_NAME_MAX is 255. */
  gethostname(hostname, sizeof(hostname)-1);
  hostname[sizeof(hostname)-1]=0;
  printf("#Info Running on: %s\n", hostname);
}
{ time_t tm;
  time(&tm);
  printf("#Info Started at time %s", ctime(&tm));
}
printf("#Info PID: %i\n", (int)getpid());/*@\colabel{getpid}@*/
#define PRINT_PARAM(a,f) printf("#Info " #a ": " f "\n", a)
PRINT_PARAM(LENGTH, "%i");
PRINT_PARAM(NUM_BOXES, "%i");
PRINT_PARAM(ALPHA, "%g");
PRINT_PARAM(ITERATIONS_PER_CHUNK, "%lli");
PRINT_PARAM(SEED, "%li");
PRINT_PARAM(F_SLIP, FORCE_OUT_FMT);
PRINT_PARAM(INIT_OFFSET, FORCE_OUT_FMT);
PRINT_PARAM(FORCE_OUT_FMT, "%s");
PRINT_PARAM(FORCE_TYPE_NAME, "%s");
PRINT_PARAM(PRECISION, "%i");
printf("#Info force_type has %i bytes\n", (int)sizeof(force_type));

setlinebuf(stdout);    /*@\colabel{setlinebuf}\hfill@*//* Avoid buffering of output. */

FPSETPREC(PRECISION);  /*@\colabel{FPSETPREC_call}\hfill@*//* Set precision. */

srand48(SEED);         /*@\hfill@*//* Initialise random number generator. ONLY for */
/*@\hfill@*//* illustration purposes, the library functions are used. To be REPLACED. */
signal(SIGHUP, signal_handler); /*@\colabel{signal_SIGHUP}\hfill@*//* Register signal handler. */
sigemptyset(&sigmask);          /*@\hfill@*//* Prepare signal mask. */
sigaddset(&sigmask, SIGHUP);
if ((long double)(NUM_BOXES-1)*(long double)(INIT_OFFSET+F_SLIP)/*@\colabel{test_init_offset}@*/
    /(long double)F_SLIP > (long double)ULLONG_MAX) {
  fprintf(stderr, "Box position overspills. Reduce INIT_OFFSET.\n");/*@\colabel{stderr}@*/
  fprintf(stdout, "#ERROR "
                  "Box position overspills. Reduce INIT_OFFSET.\n");
  exit(0);
}
if (BOX_POS(INIT_OFFSET+F_SLIP)==BOX_POS(INIT_OFFSET)) { /*@\colabel{test_ambigious}@*/
  fprintf(stderr, "Box position ambiguous. Reduce INIT_OFFSET.\n");
  fprintf(stdout, "#ERROR "
                  "Box position ambiguous. Reduce INIT_OFFSET.\n");
  exit(0);
}


if ((site=malloc(sizeof(*site)*NUM_SITES))==NULL)/*@\colabel{malloc_site}@*/
  exit(EXIT_FAILURE); /*@\hfill@*//* Malloc has failed. */
if ((stackNext=malloc(sizeof(*stackNext)*NUM_SITES))==NULL)/*@\colabel{malloc_stackNext}@*/
  exit(EXIT_FAILURE); /*@\hfill@*//* Malloc has failed. */
if ((stackCurr=malloc(sizeof(*stackCurr)*NUM_SITES))==NULL)/*@\colabel{malloc_stackCurr}@*/
  exit(EXIT_FAILURE); /*@\hfill@*//* Malloc has failed. */

/* Initialise boxes (although that should have happened automatically,
 * as it is a global variable): */
for (i=0; i<NUM_BOXES; i++) box[i]=NULL; /*@\colabel{init_box}@*/

total_time=0;/*@\colabel{total_time_init}@*/
total_time_readjusted=0;/*@\colabel{total_time_readjusted_init}@*/
f_threshold=INIT_OFFSET;/*@\colabel{f_threshold_INIT_OFFSET}@*/

/* Initialise stacks: */
stackheightCurr=0;
stackheightNext=0;

/* Initialise lattice: */
for (i=0; i<LENGTH; i++) for (j=0; j<LENGTH; j++) { /*@\colabel{init_loop}@*/
  site[COO2INDEX(i,j)].f=INIT_OFFSET+(force_type)((double)F_SLIP*drand48());  /*@\colabel{initialise_forces}@*/
  if (site[COO2INDEX(i,j)].f>f_threshold) f_threshold=site[COO2INDEX(i,j)].f;
  site[COO2INDEX(i,j)].charge=0;
  
#ifdef OBC /*@\colabel{nb}\colabel{calc_neigbours_start}@*/
  /* Left neighbour: */
  if (i>0) site[COO2INDEX(i,j)].nb[0]=&site[COO2INDEX(i-1,j)];
  else site[COO2INDEX(i,j)].nb[0]=NULL;
  /* Right neighbour: */
  if (i<LENGTH-1) site[COO2INDEX(i,j)].nb[1]=&site[COO2INDEX(i+1,j)];
  else site[COO2INDEX(i,j)].nb[1]=NULL;
  /* Lower neighbour: */
  if (j>0) site[COO2INDEX(i,j)].nb[2]=&site[COO2INDEX(i,j-1)];
  else site[COO2INDEX(i,j)].nb[2]=NULL;
  /* Upper neighbour: */
  if (j<LENGTH-1) site[COO2INDEX(i,j)].nb[3]=&site[COO2INDEX(i,j+1)];
  else site[COO2INDEX(i,j)].nb[3]=NULL;
#define BC_NAME "OBC"
#elif defined(PBC)
  /* Left neighbour: */
  if (i>0) site[COO2INDEX(i,j)].nb[0]=&site[COO2INDEX(i-1,j)];
  else site[COO2INDEX(i,j)].nb[0]=&site[COO2INDEX(LENGTH-1,j)];
  /* Right neighbour: */
  if (i<LENGTH-1) site[COO2INDEX(i,j)].nb[1]=&site[COO2INDEX(i+1,j)];
  else site[COO2INDEX(i,j)].nb[1]=&site[COO2INDEX(0,j)];;
  /* Lower neighbour: */
  if (j>0) site[COO2INDEX(i,j)].nb[2]=&site[COO2INDEX(i,j-1)];
  else site[COO2INDEX(i,j)].nb[2]=&site[COO2INDEX(i,LENGTH-1)];
  /* Upper neighbour: */
  if (j<LENGTH-1) site[COO2INDEX(i,j)].nb[3]=&site[COO2INDEX(i,j+1)];
  else site[COO2INDEX(i,j)].nb[3]=&site[COO2INDEX(i,0)];/*@\colabel{calc_neigbours_stop}@*/
#define BC_NAME "PBC"
#else
#error "No boundary condition defined"
#endif
  
  /* bx_right and bx_left are initialised automatically with box_insert: */
  box_insert(&site[COO2INDEX(i,j)]);
}
printf("#Info Boundary condition: %s\n", BC_NAME); /*@\colabel{output_BC_NAME}@*/
f_threshold_start=f_threshold; /*@\colabel{f_threshold_start}@*/
for (chunk=1; ((chunk<=NUM_CHUNKS) || (NUM_CHUNKS<1)); chunk++) { /*@\colabel{indefinite_run}@*/
  MOMENTS_INIT(siz); /*@\colabel{moments_init_start}@*/
  MOMENTS_INIT(dur);
  MOMENTS_INIT(deg); /*@\colabel{moments_init_stop}@*/
  chunk_start=total_time; /*@\colabel{chunk_start}@*/
  chunk_start_readjusted=total_time_readjusted; 
  for (it=1; it<=ITERATIONS_PER_CHUNK; it++) {/*@\colabel{start_of_it_loop}@*/
    /* Find sites with maximum force: */
    box_getmaxima(stackCurr, &stackheightCurr); /*@\colabel{getmaxima}@*/
    MOMENTS(deg, stackheightCurr); /*@\colabel{moment_deg}@*/
    for (i=0; i<stackheightCurr; i++) 
      /* Delete from boxes (sites in stack are not in boxes): */
      box_delete(stackCurr[i]);
    /*@\colabel{trigger_criterion_start}@*//* Adjust force threshold of OFC model; every site with that force
     * or greater should topple. */
    f_threshold=stackCurr[0]->f;/*@\colabel{trigger_criterion_end}@*/
    /* Total macroscopic time without readjustment: */
    total_time=(f_threshold_start-f_threshold);
    if (READJUST_CRITERION(f_threshold)) { /*@\hfill@*/
      /* Readjust all forces, as they have fallen below the threshold: */
      printf("#Info Readjusting.\n");       /*@\colabel{status_readjusting}@*/
      /* Charge to all sites that pushes minimum to INIT_OFFSET: */
      charge=INIT_OFFSET+F_SLIP-f_threshold;/*@\colabel{ramp_up}\hfill@*/
      /* Reset force threshold, using the same arithmetic operation 
       * as for all forces: */
      f_threshold+=charge;
      f_threshold_start=f_threshold;    /*@\hfill@*//* Reset time keeping. */
      total_time_readjusted+=total_time;/*@\colabel{total_time_readjusted}@*/
      total_time=0;                     /*@\colabel{total_time_reset}@*/
      for (i=0; i<NUM_BOXES; i++) box[i]=NULL;/*@\colabel{flush_boxes}\hfill@*//* Flush boxes. */
      for (i=0; i<NUM_SITES; i++) {
	site[i].f+=charge;                      
	site[i].bx_right=NULL;
	site[i].bx_left=NULL;
	box_insert(&site[i]);           /*@\colabel{re-insert}\hfill@*//* Insert all sites into boxes. */
      }
      for (i=0; i<stackheightCurr; i++)
	box_delete(stackCurr[i]);       /*@\colabel{delete_again}\hfill@*//* Delete from boxes toppling sites. */
    }
    s=t=0;                                          /*@\colabel{init_s_and_t}\hfill@*//* Avalanche size and durations. */
    do {                                            /*@\colabel{start_do_loop}@*/
      t++; 					    /*@\colabel{time_increment}\hfill@*//* Time ticks away. */
      for (i=0; i<stackheightCurr; i++) {           /*@\colabel{stackCurr_loop}@*/
	/* Force at toppling site above current effective 0: */
	d=stackCurr[i]->f - (f_threshold - F_SLIP); /*@\colabel{force_diff}@*/
	stackCurr[i]->f=f_threshold - F_SLIP;/*@\colabel{force_reset}\hfill@*//* Reset to current 0. */
	/* Insert toppled site back into boxes, as site was taken out
	 * of box when placed on stack: */
	box_insert(stackCurr[i]);
	charge=CHARGE(d);                           /*@\colabel{neighbour_charge}\hfill@*//* To be added later. */
	for (j=0; j<NUM_NEIGHBOURS; j++)            /*@\colabel{loop_neighbour_for}\hfill@*//* Visit all neighbours. */
	  if (stackCurr[i]->nb[j]!=NULL)            /*@\colabel{loop_neighbour}\hfill@*//* Neighbour exists. */
	    stackCurr[i]->nb[j]->charge+=charge;    /*@\colabel{charge_force}@*/
      } /*@\hfill@*//* End of loop over stackCurr. */
/*@\hspace*{0.7\codeheadercommentindent}@*//* All sites in boxes at this stage. Now charge. */
      for (i=0; i<stackheightCurr; i++) { 
        /* Visit all neighbours: */
	for (j=0; j<NUM_NEIGHBOURS; j++)       
	  if (stackCurr[i]->nb[j]!=NULL) {     /*@\hfill@*//* Neighbour exists. */
	    nb=stackCurr[i]->nb[j];            /*@\hfill@*//* To ease notation. */
	    /* Ignore sites visited already: */
	    if (nb->charge==0) continue;/*@\colabel{check_charge}\hfill@*/
	    box_delete(nb);                    /*@\colabel{box_delete_at_charge}\hfill@*//* Delete from box before update. */
	    /* On the stack goes >= not just >, because f=f_threshold has
	     * initiated the avalanche: */
	    if ((nb->f+=nb->charge)>=f_threshold) PUSH(nb); /*@\colabel{add_charge}@*/
	    else box_insert(nb);                            /*@\colabel{box_insert_at_charge}@*/
	    nb->charge=0;   /*@\colabel{charge_reset}@*/
	  }
      }
      s+=stackheightCurr;                     /*@\colabel{add_to_s}\hfill@*//* Avalanche size. */
/*@\hspace*{\codeheadercommentindent}@*//* Swap stacks. */
      { struct lattice **stackTemp;
	stackheightCurr=stackheightNext;  /*@\colabel{stack_swap_start}@*/
	stackheightNext=0;
	stackTemp=stackCurr;
	stackCurr=stackNext;
	stackNext=stackTemp; /*@\colabel{stack_swap_end}@*/
      }
    } while (stackheightCurr);
/*@\hspace*{\codeheadercommentindent}@*//* Avalanche has finished, collect statistics. */
    MOMENTS(siz, s);
    MOMENTS(dur, t);
  } /*@\hfill@*//* End of iterations over s aingle chunk; output follows. */
  sigprocmask(SIG_BLOCK, &sigmask, NULL);    /*@\colabel{sigprocmask_block}@*/
  { time_t tm;
    time(&tm);
    printf("#Info Chunk %i at total time " FORCE_OUT_FMT 
           " (%Lg F_SLIP) finished at %s", 
      chunk, total_time+total_time_readjusted, 
      ((long double)(total_time+total_time_readjusted))/((long double)F_SLIP),
      ctime(&tm));
  }
  MOMENTS_OUT(siz);
  MOMENTS_OUT(dur);
  MOMENTS_OUT(deg);
  sigprocmask(SIG_UNBLOCK, &sigmask, NULL);    /*@\colabel{sigprocmask_unblock}@*/
  /*@\hfill@*//* End of output. */
  if (global_stop_flg) { /*@\colabel{global_stop}@*/
    printf("#Info Termination because of global_stop_flg.\n");
    exit(0);
  }
} /*@\hfill@*//* End of chunk loop. */
printf("#Info Termination at the end of the loop.\n");
return(0);
}


/* For documentation reasons this is at the end. */
void signal_handler(int signo) /*@\colabel{signal_handler}@*/
{
  global_stop_flg=1;   /*@\colabel{global_stop_flg_setting}@*/
}/*@\colabel{last_line}@*/