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

/******* RANDOM NUMBER GENERATOR BY ZIFF **********/
#define A 471
#define B 1586
#define C 6988
#define D 9689
#define M 16383
#define RIMAX 2147483648.0        /* = 2^31 */
#define RandomInteger (++nd, ra[nd & M] = ra[(nd-A) & M] ^ ra[(nd-B) & M] ^ ra[(nd-C) & M] ^ ra[(nd-D) & M])
void seed(long seed);
static long ra[M+1], nd;
/**************************************************/

#define STEPNUMBER 600 /* number of Monte-Carlo steps */
#define SIDE 400           
#define SIZE (SIDE*SIDE)   /* lattice size */

#define a 0.02             /* addition mutation probability */#define b 0.02             /* deletion mutation probability */
#define fgvT_A 0.1         /* template function parameters  */
#define fgvT_B 200.0001
#define fgvT_H 3.000001
#define fgvR_A 0.1         /* replicase function parameters */
#define fgvR_B 200.0001
#define fgvR_H 3.000001
#define fgvF_A 0.9         /* fidelity function parameters  */#define fgvF_B 5.0001
#define fgvF_H 2.0001#define NUKLPAR 6.0000001  /* monomerdensity dependence parameter */#define DEATH 0.001        /* probability of decay                */
#define DIFFRATE 0.0       /* diffusion rate                      */
#define INIC_N 5
#define CATCH 5

#define SAVEL 60000         /**************************/
#define RESULTL1 60000    /* intervals of data-saving */
#define RESULTL2 600000    /*                          */
#define RESULTL3 6000000    /**************************/

#define MAXL_N 70 /* in the program to give a maximum possible replicator length, */
                  /* before defining arrays                                       */

typedef char matrix1[SIDE+1][SIDE+1];                /* array definitions */
typedef char matrix2[SIDE+1][SIDE+1][MAXL_N+1];
typedef long int vector1[SIZE+1];
typedef long int vector2[SIDE+1];
typedef long int vector3[9];
typedef long int vector4[5];
typedef double   vector5[MAXL_N+1];typedef long int vector6[MAXL_N+1];
typedef double restomb1[MAXL_N+1]; typedef double restomb2[MAXL_N+1][9]; 
void seed(long seed);  /* random number initialisation */
void trq_values(void); /* evaluates the template, replicase and fidelity values of replicators on each cell */
                       /* at initialisation */
void trq_update(void); /* overwrites template, replicase, fidelity values */
void inic(void);       /* inoculation of the square lattice with short replicators */
void mc_stepi(void);
void xy_choose(void);  /* picks a site at random */
void neighbourhood(void); /* counts the number and position of neighbouring empty and occupied sites among the*/
                          /* four nearest neighbours */
long randl(long);      /* random number between 0 and parameter */
double randd(void);    /* random double between 0 and 1         */
void toff_marg(void);  /* Toffoli-Margolus algorithm            */
void mc_step(void);    /* performs one Monte-Carlo step, by picking each site of the grid in a random sequence,*/
                       /* and applying the update algorithm on them */
void save(void);        /***************************/void result1(void);    /* data recording procedures */void result2(void);    /*                           */void result3(void);     /***************************/matrix1 t_n, r_n, q_n, sq_length;matrix2 sq_seq;
vector5 fgv_template, fgv_replicase, fgv_fidelity;

double t, r, q;        /* template, replicase activity, fidelity */
double n_denz;         /* density of inbuilt (not free!) monomers */

long int length, length_o, length_n;        /*variables, containing temporarily the length of replicators */
long int x, y, xs, ys;                      /* coordinates */
long int xs1, ys1, xs2, ys2;                /* coordinates of replicase and empty site */
long int max_A, max_B, max_C, max_D;        /* maximal contiguous sequnce of A,B,C,D monomers respectively */

long int i, j, k, l, m;                     /* cycle variables */
double rd;
long int rl;
double n_factor;
vector6 data;char prev_n, nukl;vector4 max_N, len_N;
vector1 availabl;
vector2 left, right, up, down;
long int available;
long int s;
int life, empty;
vector3 lifes, empties;
long int position;
int s1, s2;
FILE * tf;FILE * tf2;/***********************************************************
 * Random generator initialization                         *
 *                      by a simple Congruential generator *
 ***********************************************************/

void seed(long seed)      
{
  int  i;
 
  if(seed<0) { puts("SEED error."); exit(1); }
  ra[0]= (long) fmod(16807.0*(double)seed, 2147483647.0);
  for(i=1; i<=M; i++)
  {
    ra[i] = (long)fmod( 16807.0 * (double) ra[i-1], 2147483647.0);
  }
}

 long randl(long num)      /* random integer number between 0 and num-1 */
 {
  return(RandomInteger % num);
 }

 double randd(void)        /* random real number between 0 and 1 */
  {
   return((double) RandomInteger / RIMAX);
  }

 void trq_values(void) /* evaluates the template, replicase and fidelity values of replicators on each cell */ {                     /* at initialisation */  for(i=1; i<=SIDE; i++)  {   for(j=1; j<=SIDE; j++)   {    length=sq_length[i][j];
    if(!length)    {     t_n[i][j]=0;     r_n[i][j]=0;     q_n[i][j]=0;     continue;    }    max_N[1]=0; max_N[2]=0; max_N[3]=0; max_N[4]=0;
    len_N[0]=0; len_N[1]=0; len_N[2]=0; len_N[3]=0; len_N[4]=0;
    prev_n=0;
    for(k=1; k<=length; k++)
    {
     if(sq_seq[i][j][k]==prev_n)
     {
      len_N[prev_n]++;
     }
     else
     {
      if(len_N[prev_n]>max_N[prev_n])
      {
       max_N[prev_n]=len_N[prev_n];
      }
      len_N[prev_n]=0;      prev_n=sq_seq[i][j][k];
      len_N[prev_n]++;
     }
    }
    t_n[i][j]=max_N[1];
    r_n[i][j]=max_N[2];
    q_n[i][j]=max_N[3];
   }  } } void trq_update(void) /* overwrites template, replicase, fidelity values */ {  length=sq_length[xs2][ys2];
  max_N[1]=0; max_N[2]=0; max_N[3]=0; max_N[4]=0;
  len_N[0]=0; len_N[1]=0; len_N[2]=0; len_N[3]=0; len_N[4]=0;
  prev_n=0;
  for(k=1; k<=length; k++)
  {
   if(sq_seq[xs2][ys2][k]==prev_n)
   {
    len_N[prev_n]++;
   }
   else
   {
    if(len_N[prev_n]>max_N[prev_n])
    {
     max_N[prev_n]=len_N[prev_n];
    }
    len_N[prev_n]=0;    prev_n=sq_seq[xs2][ys2][k];
    len_N[prev_n]++;
   }
  }
  t_n[xs2][ys2]=max_N[1];
  r_n[xs2][ys2]=max_N[2];
  q_n[xs2][ys2]=max_N[3];
 } void inic(void) 
 {
  long int i;

  seed(555);

  position=0;  for(i=1; i<=SIDE; i++)
  {
   left[i]=i-1;
   right[i]=i+1;
   up[i]=i-1;
   down[i]=i+1;
  }
  left[1]=SIDE;
  up[1]=SIDE;
  right[SIDE]=1;
  down[SIDE]=1;
/* inoculation of the square lattice with short replicators */
   n_denz=0;   for(i=1; i<=SIDE; i++) 
   {
    for(j=1; j<=SIDE; j++)
    {
     sq_length[i][j]=0;
     for(k=1; k<=MAXL_N; k++)
     {
      sq_seq[i][j][k]=0;
     }
     if(randd()<0.5)
     {
      sq_length[i][j]=INIC_N;
      length=sq_length[i][j];
      n_denz+=length;      for(k=1; k<=length; k++)
      {
       sq_seq[i][j][k]=randl(4)+1;
      }
     }
    }
   }
   n_denz=n_denz/SIZE;   trq_values();  for(i=0; i<=MAXL_N; i++)  /* sets the values of template, replicase, fidelity functions */  {   fgv_template[i]=fgvT_A+(1-fgvT_A)*pow(i,fgvT_H)/(fgvT_B+pow(i,fgvT_H));
   fgv_replicase[i]=fgvR_A+(1-fgvR_A)*pow(i,fgvR_H)/(fgvR_B+pow(i,fgvR_H));
   fgv_fidelity[i]=fgvF_A+(1-fgvF_A)*pow(i,fgvF_H)/(fgvF_B+pow(i,fgvF_H));
  }
  fgv_template[0]=0;  fgv_replicase[0]=0; }

 void mc_stepi(void)
  {
   for (j=1; j<=SIZE; j++)
   {
    availabl[j]=j;
   }
   available=SIZE;
  }

 void xy_choose(void)
  {
   rl=randl(available)+1;
   s=availabl[rl];
   x=fmod((s-1),SIDE)+1;
   y=(s-1)/SIDE+1;

   availabl[rl]=availabl[available];

   available--;
  }

 void neighbourhood(void)
  {
   life=0;
   empty=0;

   for(k=1; k<=4; k++)
   {
    if      (k==1) { xs=left[x];  ys=y;       }
    else if (k==2) { xs=right[x]; ys=y;       }
    else if (k==3) { xs=x;        ys=up[y];   }
    else           { xs=x;        ys=down[y]; }

    length_n=sq_length[xs][ys];
    if(length_n>=5)
    {
     life++;
     lifes[life]=k;
    }
    else    if(length_n==0)    {
     empty++;
     empties[empty]=k;
    }
   }
  }

 void toff_marg(void)
  {
   long int xd1,xd2,xd3,xd4,xd5,xd6,xd7;
   long int yd1,yd2,yd3,yd4,yd5,yd6,yd7;
   long int temp_l;
   double temp_t, temp_r, temp_q;
   char temp_seq[MAXL_N+1];

   rl=randl(8)+1;
   switch(rl)
   {
    case 1: case 2:
     xd1=randl(SIDE)+1; yd1=randl(SIDE)+1;
     xd2=xd1;           yd2=down[yd1];
     xd3=left[xd1];     yd3=down[yd1];
     xd4=left[xd1];     yd4=yd1;
     xd5=xd1;           yd5=up[yd1];
     xd6=right[xd1];    yd6=up[yd1];
     xd7=right[xd1];    yd7=yd1;
    break;
    case 3: case 4:
     xd1=randl(SIDE)+1; yd1=randl(SIDE)+1;
     xd2=left[xd1];     yd2=yd1;
     xd3=left[xd1];     yd3=up[yd1];
     xd4=xd1;           yd4=up[yd1];
     xd5=right[xd1];    yd5=yd1;
     xd6=right[xd1];    yd6=down[yd1];
     xd7=xd1;           yd7=down[yd1];
    break;
    case 5: case 6:
     xd1=randl(SIDE)+1; yd1=randl(SIDE)+1;
     xd2=xd1;           yd2=down[yd1];
     xd3=left[xd1];     yd3=down[yd1];
     xd4=left[xd1];     yd4=yd1;
     xd5=xd1;           yd5=up[yd1];
     xd6=right[xd1];    yd6=up[yd1];
     xd7=right[xd1];    yd7=yd1;
    break;
    case 7: case 8:
     xd1=randl(SIDE)+1; yd1=randl(SIDE)+1;
     xd2=left[xd1];     yd2=yd1;
     xd3=left[xd1];     yd3=up[yd1];
     xd4=xd1;           yd4=up[yd1];
     xd5=right[xd1];    yd5=yd1;
     xd6=right[xd1];    yd6=down[yd1];
     xd7=xd1;           yd7=down[yd1];
    break;
   }

   switch(rl)
   {
    case 2: case 4: case 6: case 8: /* rotate left */
     length=sq_length[xd1][yd1];
     temp_l=length;
     temp_t=t_n[xd1][yd1]; temp_r=r_n[xd1][yd1]; temp_q=q_n[xd1][yd1];
     for(k=1; k<=length; k++)
     {
      temp_seq[k]=sq_seq[xd1][yd1][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }

     length=sq_length[xd2][yd2]; /* 1. exchange */
     sq_length[xd1][yd1]=sq_length[xd2][yd2];
     t_n[xd1][yd1]=t_n[xd2][yd2]; r_n[xd1][yd1]=r_n[xd2][yd2]; q_n[xd1][yd1]=q_n[xd2][yd2];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd1][yd1][k]=sq_seq[xd2][yd2][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd1][yd1][k]=0;
     }

     length=sq_length[xd3][yd3]; /* 2. exchange */
     sq_length[xd2][yd2]=sq_length[xd3][yd3];
     t_n[xd2][yd2]=t_n[xd3][yd3]; r_n[xd2][yd2]=r_n[xd3][yd3]; q_n[xd2][yd2]=q_n[xd3][yd3];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd2][yd2][k]=sq_seq[xd3][yd3][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd2][yd2][k]=0;
     }

     length=sq_length[xd4][yd4]; /* 3. exchange */
     sq_length[xd3][yd3]=sq_length[xd4][yd4];
     t_n[xd3][yd3]=t_n[xd4][yd4]; r_n[xd3][yd3]=r_n[xd4][yd4]; q_n[xd3][yd3]=q_n[xd4][yd4];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd3][yd3][k]=sq_seq[xd4][yd4][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd3][yd3][k]=0;
     }

     length=temp_l;            /* 4. exchange */
     sq_length[xd4][yd4]=temp_l;
     t_n[xd4][yd4]=temp_t; r_n[xd4][yd4]=temp_r; q_n[xd4][yd4]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd4][yd4][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd4][yd4][k]=0;
     }
     break;

    case 1: case 3: case 5: case 7: /* rotate right */
     length=sq_length[xd1][yd1];
     temp_l=length;
     temp_t=t_n[xd1][yd1]; temp_r=r_n[xd1][yd1]; temp_q=q_n[xd1][yd1];
     for(k=1; k<=length; k++)
     {
      temp_seq[k]=sq_seq[xd1][yd1][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }

     length=sq_length[xd4][yd4]; /* 1. exchange */
     sq_length[xd1][yd1]=sq_length[xd4][yd4];
     t_n[xd1][yd1]=t_n[xd4][yd4]; r_n[xd1][yd1]=r_n[xd4][yd4]; q_n[xd1][yd1]=q_n[xd4][yd4];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd1][yd1][k]=sq_seq[xd4][yd4][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd1][yd1][k]=0;
     }

     length=sq_length[xd3][yd3]; /* 2. exchange */
     sq_length[xd4][yd4]=sq_length[xd3][yd3];
     t_n[xd4][yd4]=t_n[xd3][yd3]; r_n[xd4][yd4]=r_n[xd3][yd3]; q_n[xd4][yd4]=q_n[xd3][yd3];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd4][yd4][k]=sq_seq[xd3][yd3][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd4][yd4][k]=0;
     }

     length=sq_length[xd2][yd2]; /* 3. exchange */
     sq_length[xd3][yd3]=sq_length[xd2][yd2];
     t_n[xd3][yd3]=t_n[xd2][yd2]; r_n[xd3][yd3]=r_n[xd2][yd2]; q_n[xd3][yd3]=q_n[xd2][yd2];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd3][yd3][k]=sq_seq[xd2][yd2][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd3][yd3][k]=0;
     }

     length=temp_l;            /* 4. exchange */
     sq_length[xd2][yd2]=temp_l;
     t_n[xd2][yd2]=temp_t; r_n[xd2][yd2]=temp_r; q_n[xd2][yd2]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd2][yd2][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd2][yd2][k]=0;
     }
     break;
   }
     length=sq_length[xd2][yd2];
     temp_l=length;
     temp_t=t_n[xd2][yd2]; temp_r=r_n[xd2][yd2]; temp_q=q_n[xd2][yd2];
     for(k=1; k<=length; k++)
     {
      temp_seq[k]=sq_seq[xd2][yd2][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      temp_seq[k]=0;
     }
     length=sq_length[xd7][yd7];
     sq_length[xd2][yd2]=sq_length[xd7][yd7];
     t_n[xd2][yd2]=t_n[xd7][yd7]; r_n[xd2][yd2]=r_n[xd7][yd7]; q_n[xd2][yd2]=q_n[xd7][yd7];
     for(k=1; k<=length; k++)
     {
      sq_seq[xd2][yd2][k]=sq_seq[xd7][yd7][k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd2][yd2][k]=0;
     }
     length=temp_l;
     sq_length[xd7][yd7]=temp_l;
     t_n[xd7][yd7]=temp_t; r_n[xd7][yd7]=temp_r; q_n[xd7][yd7]=temp_q;
     for(k=1; k<=length; k++)
     {
      sq_seq[xd7][yd7][k]=temp_seq[k];
     }
     for(k=(length+1); k<=MAXL_N; k++)
     {
      sq_seq[xd7][yd7][k]=0;
     }
  
     length=sq_length[xd4][yd4];
     temp_l=length;
     temp_t=t_n[xd4][yd4]; temp_r=r_n[xd4][yd4]; temp_q=q_n[xd4][yd4];
     for(k=1; k<=length; k++)
     {
      temp_seq[k]=sq_seq[xd4][yd4][k];