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

#define MAX_X 500
#define MAX_Y 500

#define STEPS_MAX 200
#define UNCHANGED_MAX 10

void main()
{
    int   grid[MAX_X+2][MAX_Y+2];	/* grid of vegetation values */
    int   nx;			/* x dimension of grid */
    int   ny;			/* y dimension of grid */
    int   max_steps;		/* max # timesteps to simulate */
    int   max_unchanged;	/* max # timesteps with no vegetation change */
    int   vegies;		/* amount of stable vegetation */
    int   nsteps;		/* number of steps actually run */
    int   nsims;		/* number of simulations to perform */
    int   ndied;		/* # populations which die out */
    int   nunsettled;		/* # populations which don't stablize */
    int   nstable;		/* # populations which do stablize */
    float tot_steps_stable;	/* total/average steps to stablization */
    float tot_veg_stable;	/* total/average stable vegetation */
    double prob;		/* population probability */
    int   seed, seed0;		/* random number seeds */
    int   i;			/* loop counter */
    void  initialize_grid();
    int   game_o_life();

    /* First read in all the parameters. */

    nx = MAX_X+1;
    while (nx > MAX_X || ny > MAX_Y) {
        printf("Enter X and Y dimensions of wilderness: ");
        scanf("%d%d", &nx, &ny);
    }

    printf("Enter population probability: ");
    scanf("%lf", &prob);

    printf("Enter number of simulations: ");
    scanf("%d", &nsims);

    printf("Enter random number seed: ");
    scanf("%d", &seed0);

    ndied = 0;
    nunsettled = 0;
    nstable = 0;
    tot_steps_stable = 0;
    tot_veg_stable = 0;

    for (i = 1; i <= nsims; i++) {

        /* Initialize the grid values using the given probability. */

	seed = seed0 * i;
        initialize_grid(grid, nx, ny, seed, prob);

        /* Now run the game of life simulation, returning the number of steps. */

        max_steps = STEPS_MAX;
        max_unchanged = UNCHANGED_MAX;
        nsteps = game_o_life(grid, nx, ny, max_steps, max_unchanged, &vegies);

        printf("Number of time steps = %d, Vegetation total = %d\n", nsteps, vegies);
	if (vegies == 0) ndied = ndied + 1;
	else if (nsteps >= max_steps) nunsettled = nunsettled + 1;
	else {
	    nstable = nstable + 1;
	    tot_steps_stable = tot_steps_stable + nsteps;
	    tot_veg_stable = tot_veg_stable + vegies;
	}
    }
    if (nstable > 0) {
        tot_steps_stable = tot_steps_stable / nstable;
        tot_veg_stable = tot_veg_stable / nstable;
    }
    printf("Percentage which died out: %g%%\n", 100.0*ndied/nsims);
    printf("Percentage unsettled:      %g%%\n", 100.0*nunsettled/nsims);
    printf("Percentage stablized:      %g%%\n", 100.0*nstable/nsims);
    printf("  Of which:\n");
    printf("  Average steps:           %g\n", tot_steps_stable);
    printf("  Average vegetation:      %g\n", tot_veg_stable);
}


void initialize_grid(grid, nx, ny, seed, prob)
int   grid[MAX_X+2][MAX_Y+2];	/* grid of vegetation values */
int   nx;			/* x dimension of grid */
int   ny;			/* y dimension of grid */
int   seed;			/* random number seed */
double prob;			/* population probability */
{
    int i, j;			/* loop counters */
    int index;			/* unique value for each grid cell */
    int new_seed;		/* unique seed for each grid point */
    double rand1();

    for (i=1; i<=nx; i++) {
	for (j=1; j<=ny; j++) {
	    index = ny*i + j;
	    new_seed = seed + index;
	    if (rand1(&new_seed) > prob) grid[i][j] = 0;
	    else grid[i][j] = 1;
	}
    }
}


int game_o_life(grid, nx, ny, max_steps, max_unchanged, pvegies)
int   grid[MAX_X+2][MAX_Y+2];	/* grid of vegetation values */
int   nx;			/* x dimension of grid */
int   ny;			/* y dimension of grid */
int   max_steps;		/* max # timesteps to simulate */
int   max_unchanged;		/* max # timesteps with no vegetation change */
int  *pvegies;			/* final vegetation level */
{
    int   step;			/* counts the time steps */
    int   converged;		/* has the vegetation stablized? */
    int   n_unchanged;		/* # timesteps with no vegetation change */
    int   old_vegies;		/* previous level of vegetation */
    int   old2_vegies;		/* previous level of vegetation */
    int   old3_vegies;		/* previous level of vegetation */
    int   vegies;		/* total amount of vegetation */
    int   neighbors;		/* quantity of neighboring vegetation */
    int   temp_grid[MAX_X][MAX_Y];	/* grid to hold updated values */
    int i, j;			/* loop counters */

    step = 1;
    vegies = 1;
    old_vegies = -1;
    old2_vegies = -1;
    old3_vegies = -1;
    n_unchanged = 0;
    converged = 0;

    while (!converged && vegies > 0 && step < max_steps) {

	/* Count the total amount of vegetation. */

	vegies = 0;
	for (i=1; i<=nx; i++) {
	    for (j=1; j<=ny; j++) {
		vegies = vegies + grid[i][j];
	    }
	}
	if (vegies == old_vegies || vegies == old2_vegies
	    || vegies == old3_vegies) {
	    n_unchanged = n_unchanged + 1;
	    if (n_unchanged >= max_unchanged) converged = 1;
	}
	else {
	    n_unchanged = 0;
	}
	old3_vegies = old2_vegies;
	old2_vegies = old_vegies;
	old_vegies = vegies;

/* printf(" step %d: vegies = %d\n", step, vegies);*/

	if (!converged) {
            /* Copy the sides of the grid to make torus simple. */
            for (i = 1; i <= nx; i++) {
	        grid[i][0] = grid[i][ny];
	        grid[i][ny+1] = grid[i][1];
            }

            for (j = 0; j <= ny + 1; j++) {
	        grid[0][j] = grid[nx][j];
	        grid[nx+1][j] = grid[1][j];
            }


            /* Now run one time step, putting result in temp_grid. */

	    for (i=1; i<=nx; i++) {
	        for (j=1; j<=ny; j++) {
		    neighbors = grid[i-1][j-1] + grid[i-1][j] + grid[i-1][j+1]
		              + grid[i][j-1] + grid[i][j+1]
		              + grid[i+1][j-1] + grid[i+1][j] + grid[i+1][j+1];
		    temp_grid[i][j] = grid[i][j];
		    if (neighbors >= 25 || neighbors <= 3) {
		        temp_grid[i][j] = temp_grid[i][j] - 1;
			if (temp_grid[i][j] < 0) temp_grid[i][j] = 0;
		    }
		    else if (neighbors <= 15) {
		        temp_grid[i][j] = temp_grid[i][j] + 1;
			if (temp_grid[i][j] > 10) temp_grid[i][j] = 10;
		    }
	        }
	    }


	    /* Now copy temp_grid back to grid. */

	    for (i=1; i<=nx; i++) {
	        for (j=1; j<=ny; j++) {
		    grid[i][j] = temp_grid[i][j];
	        }
	    }
	    step = step + 1;
	}
    }

    *pvegies = vegies;
    return(step);
}


double    rand1(iseed)
int      *iseed;
{
    double aa=16807.0;
    double mm=2147483647.0;
    double sseed;
    int    jseed;
    int    i;

    jseed = *iseed;
    for (i = 1; i <= 5; i++) {
       sseed = jseed;
       jseed = aa*sseed/mm;
       sseed = aa*sseed - mm*jseed;
       jseed = sseed;
    }

    *iseed = jseed;

    return(sseed/mm);
}