//TODO: We could probably get a significant speedup by doing our processing in
//bulk via a lookup table - i.e., running all life cycle ops on a 3x3 or 3x4
//by turning it into binary, using that as a hash table entry for a
//pre-computed calculation for that combination.  3x3 is no problem.  3x4 will
//take 2 gig normally, so we'd have to do some "cheats", taking advantage of
//symmetry to reduce it, since we don't have that much ram and don't want to
//swap out.

#define USE_ZLIB

#define X_RES		40000	//X-res must be divisible by 8 (for ease in saving/loading)
#define Y_RES		40000
#define VERSION		0.1
#define STEPS_PER_SAVE	100
#define ENERGY_PER_TIME	5000
#define ENERGY_MULT	0.3
#define ENERGY_EXP	0.5
#define MAX_PATH	4096
#ifdef USE_ZLIB
  #define OUTFILE 	"mc.sav.gz"
#else
  #define OUTFILE	"mc.sav"
#endif
#define VIEW_WIDTH	75
#define VIEW_HEIGHT	35
#define VIEW_CHANGE	400


#define TRUE 1
#define FALSE 0
#define NULL 0


#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#ifdef USE_ZLIB
  #include <zlib.h>
#endif


const unsigned char ON_BIT[] = { 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80 };
const unsigned char OFF_BIT[] = { 0xFE, 0xFD, 0xFB, 0xF7, 0xEF, 0xDF, 0xBF, 0x7F };

float version=VERSION;
unsigned long x_res=X_RES;
unsigned long y_res=Y_RES;
unsigned long x_resB;
unsigned char* world=NULL;
unsigned char* buffer=NULL;
unsigned long long steps=0;

unsigned long view_x,view_y;

float frand()	//Returns a random number from 0 to 1.
{
    return (((float)random())/((float)RAND_MAX));
}

double drand()	//Returns a random number from 0 to 1.
{
    return (((double)random())/(double)RAND_MAX*0.99999)+
           (((double)random())/(double)RAND_MAX*0.00001);
}

double curve_rand(double center)  //Returns a random number on a centered exponential distribution curve
{
    return center * -log(drand());
}


int main(int argc, char** argv)
{
    unsigned long x,y,x2,y2,x3,y3;
    unsigned long i;
    float f1,f2;
    char neighbors[8];
    
    printf("Initializing...\n");
    
    srand(time(0));
    view_x=(int)(frand()*X_RES);
    view_y=(int)(frand()*Y_RES);
    printf("Setting view to: (%d,%d)-(%d,%d)\n",view_x,view_y,view_x+VIEW_WIDTH,view_y+VIEW_WIDTH);
    
    if (argc==2)
    {
#ifdef USE_ZLIB
        gzFile* fp;

        if (!(fp=gzopen(argv[1],"r")))
        {
            printf("Unable to open file %s\n",argv[1]);
            return 1;
        }
        
        printf("Loading file %s\n",argv[1]);
        
        gzread(fp,&version,4);
        gzread(fp,&x_res,4);
        gzread(fp,&y_res,4);
        gzread(fp,&steps,8);
        x_resB = x_res/8;
        
        world=malloc(x_resB*y_res);
        buffer=malloc(x_resB*2);
        
        for (y=0; y<y_res; y++)
        {
            for (x=0; x<x_resB; x++)
            {
                world[x+y*x_resB]=gzgetc(fp);
            }
        }
        
        gzclose(fp);
#else        
        FILE* fp;

        if (!(fp=fopen(argv[1],"r")))
        {
            printf("Unable to open file %s\n",argv[1]);
            return 1;
        }
        
        printf("Loading file %s\n",argv[1]);
        
        fread(&version,4,1,fp);
        fread(&x_res,4,1,fp);
        fread(&y_res,4,1,fp);
        fread(&steps,8,1,fp);
        x_resB = x_res/8;
        
        world=malloc(x_resB*y_res);
        buffer=malloc(x_resB*2);
        
        for (y=0; y<y_res; y++)
        {
            for (x=0; x<x_resB; x++)
            {
                world[x+y*x_resB]=fgetc(fp);
            }
        }
        
        fclose(fp);
#endif
    }
    else
    {
        x_resB = x_res/8;
        
        world=malloc(x_resB*y_res);
        buffer=malloc(x_resB*2);
        
        for (y=0; y<y_res; y++)
        {
            for (x=0; x<x_resB; x++)
            {
                world[x+y*x_resB]=random()&0xFF;
            }
        }
    }

    printf("Done.  Starting.\n");

    view_x=(int)(frand()*X_RES);
    view_y=(int)(frand()*Y_RES);
    printf("Setting view to: (%d,%d)-(%d,%d)\n\n",view_x,view_y,view_x+VIEW_WIDTH,view_y+VIEW_WIDTH);
    
    
    //Loop forever
    while (TRUE)
    {
        unsigned long num_inputs;
        int count;
        
        //Dump energy input into the system
        for (i=0; i<ENERGY_PER_TIME; i++)
        {
            double loc_x,loc_y;
            double radius;
            double dist;
            int bit;
            unsigned long rnd;	//We want to, for speed, minimize the number of times that we call random().  Here, we make sure to use every bit of the number before calling random again.
            
            loc_x=frand()*x_res;
            loc_y=frand()*y_res;
            radius=sqrt(ENERGY_MULT/pow(frand(),ENERGY_EXP))+1;
            for (y=floor(loc_y-radius/2); y<floor(loc_y+radius/2); y++)
            {
                y2=y%y_res;
                rnd=random();
                f1=floor(loc_x-radius/2);
                f2=floor(loc_x+radius/2);
                for (x=f1, bit=x%8, x2=x/8; x<f2; x++)
                {
                    world[x2+y2*x_resB]=(world[x2+y2*x_resB] & OFF_BIT[bit]) | (rnd & ON_BIT[bit]);
                    
                    bit++;
                    if (bit==8)
                    {
                        bit=0;
                        x2++;
                        if (x2==x_resB) //I'll trust the compiler to speed this up  :)
                           x2=0;
                        rnd=rnd>>8;
                        if (rnd==0)
                            rnd=random();
                    }
                }
            }
        }
    
    
        //Process using Conway's rules.
        for (y=0; y<y_res; y++)
        {
            if (y>0)
                y2=y-1;
            else
                y2=y_res-1;
            if (y!=y_res-1)
                y3=y+1;
            else
                y3=0;

            //Fill in a row of the buffer.
            for (x=0; x<x_resB; x++)	//We process 8 cells at once.
            {
                if (x>0)
                    x2=x-1;
                else
                    x2=x_resB-1;
                if (x!=x_resB-1)
                    x3=x+1;
                else
                    x3=0;
                neighbors[0]=((world[x2+y2*x_resB]&ON_BIT[7])>0)+
                             ((world[x+y2*x_resB]&ON_BIT[0])>0)+
                             ((world[x+y2*x_resB]&ON_BIT[1])>0)+
                             ((world[x2+y*x_resB]&ON_BIT[7])>0)+
                             ((world[x+y*x_resB]&ON_BIT[1])>0)+
                             ((world[x2+y3*x_resB]&ON_BIT[7])>0)+
                             ((world[x+y3*x_resB]&ON_BIT[0])>0)+
                             ((world[x+y3*x_resB]&ON_BIT[1])>0);
                for (i=1; i<7; i++)
                {
                    neighbors[i]=((world[x+y2*x_resB]&ON_BIT[i-1])>0)+
                                 ((world[x+y2*x_resB]&ON_BIT[i])>0)+
                                 ((world[x+y2*x_resB]&ON_BIT[i+1])>0)+
                                 ((world[x+y*x_resB]&ON_BIT[i-1])>0)+
                                 ((world[x+y*x_resB]&ON_BIT[i+1])>0)+
                                 ((world[x+y3*x_resB]&ON_BIT[i-1])>0)+
                                 ((world[x+y3*x_resB]&ON_BIT[i])>0)+
                                 ((world[x+y3*x_resB]&ON_BIT[i+1])>0);
                }
                neighbors[7]=((world[x+y2*x_resB]&ON_BIT[6])>0)+
                             ((world[x+y2*x_resB]&ON_BIT[7])>0)+
                             ((world[x3+y2*x_resB]&ON_BIT[0])>0)+
                             ((world[x+y*x_resB]&ON_BIT[6])>0)+
                             ((world[x3+y*x_resB]&ON_BIT[0])>0)+
                             ((world[x+y3*x_resB]&ON_BIT[6])>0)+
                             ((world[x+y3*x_resB]&ON_BIT[7])>0)+
                             ((world[x3+y3*x_resB]&ON_BIT[0])>0);
                buffer[x+x_resB]=0;	//Set them all dead
                for (i=0; i<8; i++)
                {
                    if (neighbors[i]==3)	//Create it no matter what
                        buffer[x+x_resB]|=1<<i;
                    else if (neighbors[i]==2)	//Leave it as it was.
                        buffer[x+x_resB]|=world[x+y*x_resB]&ON_BIT[i];
	            else
	            	; //Do nothing - leave it a dead cell.                        
                }
            }

            //Write out a row of the buffer.
            if (y>=1)
            {
                for (x=0; x<x_resB; x++)
                {
                    world[x+(y-1)*x_resB]=buffer[x];
                    buffer[x]=buffer[x+x_resB];
                }
            }
        }
        
        //Write out that last buffer row.
        for (x=0; x<x_resB; x++)
        {
            world[x+(y_res-1)*x_resB]=buffer[x];
        }
        
        steps++;
        
        //Check to see if we need to save.
        if (!(steps % STEPS_PER_SAVE))
        {
#ifdef USE_ZLIB
            gzFile* fp;
            char outfile[MAX_PATH];
            
            sprintf(outfile,OUTFILE,steps);
            printf("Writing file %s\n",outfile);
            
            if (!(fp=gzopen(outfile,"w")))
            {
                printf("Unable to open file %s\n",argv[1]);
                return 1;
            }
        
            gzwrite(fp,&version,4);
            gzwrite(fp,&x_res,4);
            gzwrite(fp,&y_res,4);
            gzwrite(fp,&steps,8);

            for (y=0; y<y_res; y++)
            {
                for (x=0; x<x_resB; x++)
                {
                    gzputc(fp,world[x+y*x_resB]);
                }
            }
            
            gzclose(fp);
#else
            FILE* fp;
            char outfile[MAX_PATH];
            
            sprintf(outfile,OUTFILE,steps);
            printf("Writing file %s\n",outfile);
            
            if (!(fp=fopen(outfile,"w")))
            {
                printf("Unable to open file %s\n",argv[1]);
                return 1;
            }
        
            fwrite(&version,4,1,fp);
            fwrite(&x_res,4,1,fp);
            fwrite(&y_res,4,1,fp);
            fwrite(&steps,8,1,fp);

            for (y=0; y<y_res; y++)
            {
                for (x=0; x<x_resB; x++)
                {
                    fputc(world[x+y*x_resB],fp);
                }
            }
            
            fclose(fp);
#endif
        }
        
        //Draw output from our viewing window
        printf("\n\n");
        printf("Step: %d\n",steps);
        for (y=view_y; y<view_y+VIEW_HEIGHT; y++)
        {
            y2=y%y_res;
            for (x=view_x; x<view_x+VIEW_WIDTH; x++)
            {
                x2=x%x_res;
                printf("%c",((world[x2/8+y2*x_resB]&ON_BIT[x2%8])?'*':' '));
            }
            printf("\n");
        }
        printf("\n");
        
        //Change the view if needed.
        if (!(steps%VIEW_CHANGE))
        {
            view_x=(int)(frand()*X_RES);
            view_y=(int)(frand()*Y_RES);
            printf("Changing view to: (%d,%d)-(%d,%d)\n",view_x,view_y,view_x+VIEW_WIDTH,view_y+VIEW_WIDTH);
        }
    }
    return 0;
}