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

#ifdef __APPLE__
    #include <GLUT/glut.h>
#else
    #include <GL/glut.h>
#endif

#define WINDOW_WIDTH    500
#define WINDOW_HEIGHT   500

#define LATTICE_WIDTH	32
#define LATTICE_HEIGHT	32
#define LATTICE_DEPTH	32

#define RANDOM_CA_WIDTH	100

#define POINT_SIZE      2

#define DEBUG           1

int lattice[LATTICE_DEPTH][LATTICE_HEIGHT][LATTICE_WIDTH][24];
int tmp_lattice[LATTICE_DEPTH][LATTICE_HEIGHT][LATTICE_WIDTH][24];
int c[24][4] = { 
				{1, 1, 0, 0}, {1, -1, 0, 0}, {-1, 1, 0, 0}, {-1, -1, 0, 0},
				{1, 0, 1, 0}, {1, 0, -1, 0}, {-1, 0, 1, 0}, {-1, 0, -1, 0},
				{1, 0, 0, 1}, {1, 0, 0, -1}, {-1, 0, 0, 1}, {-1, 0, 0, -1},
				{0, 1, 1, 0}, {0, 1, -1, 0}, {0, -1, 1, 0}, {0, -1, -1, 0},
				{0, 1, 0, 1}, {0, 1, 0, -1}, {0, -1, 0, 1}, {0, -1, 0, -1}, 
				{0, 0, 1, 1}, {0, 0, 1, -1}, {0, 0, -1, 1}, {0, 0, -1, -1}
			   };
int q[4];
int applied_normalization_isometries[9];

int random_ca_row[RANDOM_CA_WIDTH];
int tmp_random_ca_row[RANDOM_CA_WIDTH];

// The T# are the sigmas
enum types_of_normalization_isometries {S1, S2, S3, S4, 
                                        P12, P13, P14, P23, P24, P34, 
                                        T1, T2};

void rule30_next( ) {
	int i = 0;
	for(i = 1; i < RANDOM_CA_WIDTH-1; i++) {
		random_ca_row[i] = tmp_random_ca_row[i-1] ^ 
                           (tmp_random_ca_row[i] | tmp_random_ca_row[i+1]);
	}
	for(i = 0; i < RANDOM_CA_WIDTH; i++) {
		tmp_random_ca_row[i] = random_ca_row[i];
	}
}

int rand30(int x) {
	rule30_next( );
	int r = random_ca_row[RANDOM_CA_WIDTH/2];
	return r%x;
}

// Symmetry with respect to x_{alpha}
void s(int x, int y, int z, int alpha) {
	int i;
	alpha--;
	for(i = 0; i < 6; i++) {
		lattice[z][y][x][(alpha * (i+1)) % 24] *= -1;
	}
}

// Permutation of the alpha and beta coordinates
void p(int x, int y, int z, int alpha, int beta) {
	int i, tmp;
	alpha--; beta--;
	for(i = 0; i < 6; i++) {
		tmp = lattice[z][y][x][(alpha * (i+1)) % 24];
		lattice[z][y][x][(alpha * (i+1)) % 24] = lattice[z][y][x][(beta * (i+1)) % 24];
		lattice[z][y][x][(beta * (i+1)) % 24] = tmp;
	}
}

// reflection with respect to the hyperplane x_1 + x_4 = x_2 + x_3
void sigma_1(int x, int y, int z) {
	int i, j, k;
	int tmp[24];
	
	for(i = 0; i < 6; i++) {
		for(j = 0; j < 4; j++) {
			for(k = 0; k < 4; k++) {
				if(k == (i*3 % 4)) {
					tmp[i*4 + j] -= lattice[z][y][x][i*4 + k];
				}
				else {
					tmp[i*4 + j] += lattice[z][y][x][i*4 + k];
				}
			}
		}
	}
	
	for(i = 0; i < 24; i++) {
		lattice[z][y][x][i] = 0.5 * tmp[i];
	}
}

// reflection with respect to the hyperplane x_1 = x_2 + x_3 + x+4
void sigma_2(int x, int y, int z) {
	int i, j, k;
	int tmp[24];
	
	for(i = 0; i < 6; i++) {
		for(j = 0; j < 4; j++) {
			for(k = 0; k < 4; k++) {
				if(!(j == 0 || j == 3 || j == k))  {
					tmp[i*4 + j] -= lattice[z][y][x][i*4 + k];
				}
				else {
					tmp[i*4 + j] += lattice[z][y][x][i*4 + k];
				}
			}
		}
	}
	
	for(i = 0; i < 24; i++) {
		lattice[z][y][x][i] = 0.5 * tmp[i];
	}
}

void propagation( ) {
    int x,y,z, i;
    int toroidal_x, toroidal_y, toroidal_z;
    for(x = 0; x < LATTICE_WIDTH; x++) {
        for(y = 0; y < LATTICE_HEIGHT; y++) {
            for(z = 0; z < LATTICE_DEPTH; z++) {
                for(i = 0; i < 24; i++) {

                    toroidal_x = (x + c[i][0] + LATTICE_WIDTH)  % LATTICE_WIDTH;
                    toroidal_y = (y + c[i][1] + LATTICE_HEIGHT) % LATTICE_HEIGHT;
                    toroidal_z = (z + c[i][2] + LATTICE_DEPTH)  % LATTICE_DEPTH;

                    tmp_lattice[toroidal_z][toroidal_y][toroidal_x][i] = 
                                lattice[z][y][x][i];
                }
            }
        }
    }
    for(x = 0; x < LATTICE_WIDTH; x++) {
        for(y = 0; y < LATTICE_HEIGHT; y++) {
            for(z = 0; z < LATTICE_DEPTH; z++) {
                for(i = 0; i < 24; i++) {
                    lattice[z][y][x][i] = tmp_lattice[z][y][x][i];
                }
            }
        }
    }

}

void collision( ) {
	int x, y, z, alpha, i;
	int t = 0;
	int r = 0;
	
	for(x = 0; x < 9; x++) {
		applied_normalization_isometries[x] = 0;
	}
	
	for(x = 0; x < LATTICE_WIDTH; x++) {
		for(y = 0; y < LATTICE_HEIGHT; y++) {
			for(z = 0; z < LATTICE_DEPTH; z++) {
				q[0] = q[1] = q[2] = q[3] = 0;

				for(alpha = 0; alpha < 4; alpha++) {
					for(i = 0; i < 24; i++) {
						q[alpha] += lattice[z][y][x][i] * c[i][alpha];
					}
				}
				t = 0;
				
				if(q[0] < 0) {
					s(x, y, z, 1);
					applied_normalization_isometries[t++] = S1;
				}
				if(q[1] < 0) {
					s(x, y, z, 2);
					applied_normalization_isometries[t++] = S2;
				}
				if(q[2] < 0) {
					s(x, y, z, 3);
					applied_normalization_isometries[t++] = S3;
				}
				if(q[3] < 0) {
					s(x, y, z, 4);
					applied_normalization_isometries[t++] = S4;
				}
				if(q[0] < q[1]) {
					p(x, y, z, 1, 2);
					applied_normalization_isometries[t++] = P12;
				}
				if(q[0] < q[2]) {
					p(x, y, z, 1, 3);
					applied_normalization_isometries[t++] = P13;
				}
				if(q[0] < q[3]) {
					p(x, y, z, 1, 4);
					applied_normalization_isometries[t++] = P14;
				}
				if(q[1] < q[2]) {
					p(x, y, z, 2, 3);
					applied_normalization_isometries[t++] = P23;
				}
				if(q[1] < q[3]) {
					p(x, y, z, 2, 4);
					applied_normalization_isometries[t++] = P24;
				}
				if(q[2] < q[3]) {
					p(x, y, z, 3, 4);
					applied_normalization_isometries[t++] = P34;
				}
				if(q[3] > 0 && (q[0] + q[3] == q[1] + q[2])) {
					sigma_2(x, y, z);
					applied_normalization_isometries[t++] = T2;
				}
				if(q[3] > 0 && (q[0] + q[3] > q[1] + q[2])) {
					sigma_1(x, y, z);
					applied_normalization_isometries[t++] = T1;
					if(q[0] - q[3] > q[1] + q[2]) {
						s(x,y,z, 4);
						applied_normalization_isometries[t++] = S4;
					}
				}
				
				
				for(alpha = 0; alpha < 4; alpha++) {
					for(i = 0; i < 24; i++) {
						q[alpha] += lattice[z][y][x][i] * c[i][alpha];
					}
				}
				
				
				if(q[0] == q[1] && q[1] > q[2] && q[2] > q[3] && q[3] > 0) {
					p(x,y,z, 1,2);
				}
				else if(q[0] == q[1] && q[1] == q[2] && q[2] > q[3] && q[3] > 0) {
					r = rand30(2);
					if(r == 0) {
						p(x,y,z, 1,2);
						p(x,y,z, 2,3);
					}
					else {
						p(x,y,z, 1,3);
						p(x,y,z, 2,3);
					}
				}
				else if((q[0] > q[1] && q[1] > q[2] && q[2] > q[3] && q[3] == 0) &&
						(q[0] == q[1] + q[2])) {
					r = rand30(2);
					if(r == 0) {
						sigma_1(x,y,z);
						s(x,y,z, 4);
					}
					else {
						sigma_2(x,y,z);
						s(x,y,z, 4);
					}
				}
				else if((q[0] > q[1] && q[1] > q[2] && q[2] > q[3] && q[3] == 0) &&
						(q[0] != q[1] + q[2])) {
					s(x,y,z, 4);
				}
				else if(q[0] == q[1] && q[1] > q[2] && q[2] > q[3] && q[3] == 0) {
					p(x,y,z, 1,2);
					s(x,y,z, 4);
				}
				else if((q[0] > q[1] && q[1] == q[2] && q[2] > q[3] && q[3] == 0) &&
						(q[0] == 2*q[1])) {
					r = rand30(4);
					if(r == 0) {
						sigma_1(x,y,z);
						s(x,y,z, 4);
					}
					else if(r == 1) {
						sigma_2(x,y,z);
						s(x,y,z, 4);
					}
					else if(r == 2) {
						sigma_1(x,y,z);
						p(x,y,z, 2,3);
						s(x,y,z, 4);
					}
					else {
						sigma_2(x,y,z);
						p(x,y,z, 2,3);
						s(x,y,z, 4);
					}
				}
				else if((q[0] > q[1] && q[1] == q[2] && q[2] > q[3] && q[3] == 0) &&
						(q[0] != 2*q[1])) {
					p(x,y,z, 2,3);
					s(x,y,z, 4);
				}
				else if(q[0] == q[1] && q[1] == q[2] && q[2] > q[3] && q[3] == 0) {
					r = rand30(4);
					if(r == 0) {
						p(x,y,z, 1,2);
						p(x,y,z, 2,3);
					}
					else if(r == 1) {
						p(x,y,z, 1,3);
						p(x,y,z, 2,3);
					}
					else if(r == 2) {
						p(x,y,z, 1,2);
						p(x,y,z, 2,3);
						s(x,y,z, 4);
					}
					else {
						p(x,y,z, 1,3);
						p(x,y,z, 2,3);
						s(x,y,z, 4);
					}
				}
				else if(q[0] > q[1] && q[1] > q[2] && q[2] == q[3] && q[3] == 0) {
					r = rand30(3);
					if(r == 0) {
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else if(r == 1) {
						p(x,y,z, 3,4);
						s(x,y,z, 3);
					}
					else {
						p(x,y,z, 3,4);
						s(x,y,z, 4);
					}
				}
				else if(q[0] == q[1] && q[1] > q[2] && q[2] == q[3] && q[3] == 0) {
					r = rand30(6);
					if(r == 0) {
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 3);
					}
					else if(r == 1) {
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 4);
					}
					else if(r == 2) {
						sigma_1(x,y,z);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else if(r == 3) {
						sigma_1(x,y,z);
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else if(r == 4) {
						sigma_2(x,y,z);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else {
						sigma_2(x,y,z);
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
					}
				}
				else if(q[0] > q[1] && q[1] == q[2] && q[2] == q[3] && q[3] == 0) {
					r = rand30(6);
					if(r == 0) {
						p(x,y,z, 2,3);
						s(x,y,z, 2);
						s(x,y,z, 4);
					}
					else if(r == 1) {
						p(x,y,z, 2,3);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else if(r == 2) {
						p(x,y,z, 2,4);
						s(x,y,z, 2);
						s(x,y,z, 3);
					}
					else if(r == 3) {
						p(x,y,z, 2,4);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else if(r == 4) {
						p(x,y,z, 3,4);
						s(x,y,z, 2);
						s(x,y,z, 3);
					}
					else {
						p(x,y,z, 3,4);
						s(x,y,z, 2);
						s(x,y,z, 4);
					}
				}
				else if(q[0] == q[1] && q[1] == q[2] && q[2] == q[3] && q[3] == 0) {
					r = rand30(12);
					if(r == 0) {
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 1);
						s(x,y,z, 3);
					}
					else if(r == 1) {
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 1);
						s(x,y,z, 4);
					}
					else if(r == 2) {
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 2);
						s(x,y,z, 3);
					}
					else if(r == 3) {
						p(x,y,z, 1,2);
						p(x,y,z, 3,4);
						s(x,y,z, 2);
						s(x,y,z, 4);
					}
					else if(r == 4) {
						p(x,y,z, 1,3);
						p(x,y,z, 2,4);
						s(x,y,z, 1);
						s(x,y,z, 2);
					}
					else if(r == 5) {
						p(x,y,z, 1,3);
						p(x,y,z, 2,4);
						s(x,y,z, 1);
						s(x,y,z, 4);
					}
					else if(r == 6) {
						p(x,y,z, 1,3);
						p(x,y,z, 2,4);
						s(x,y,z, 2);
						s(x,y,z, 3);
					}
					else if(r == 7) {
						p(x,y,z, 1,3);
						p(x,y,z, 2,4);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
					else if(r == 8) {
						p(x,y,z, 1,4);
						p(x,y,z, 2,3);
						s(x,y,z, 1);
						s(x,y,z, 2);
					}
					else if(r == 9) {
						p(x,y,z, 1,4);
						p(x,y,z, 2,3);
						s(x,y,z, 1);
						s(x,y,z, 3);
					}
					else {
						p(x,y,z, 1,4);
						p(x,y,z, 2,3);
						s(x,y,z, 3);
						s(x,y,z, 4);
					}
				}
				
				while(t-- != -1) {
					if(applied_normalization_isometries[t] == S1) {
						s(x,y,z, 1);
					}
					else if(applied_normalization_isometries[t] == S2) {
						s(x,y,z, 2);
					}
					else if(applied_normalization_isometries[t] == S3) {
						s(x,y,z, 3);
					}
					else if(applied_normalization_isometries[t] == S4) {
						s(x,y,z, 4);
					}
					else if(applied_normalization_isometries[t] == P12) {
						p(x,y,z, 1,2);
					}
					else if(applied_normalization_isometries[t] == P13) {
						p(x,y,z, 1,3);
					}
					else if(applied_normalization_isometries[t] == P14) {
						p(x,y,z, 1,4);
					}
					else if(applied_normalization_isometries[t] == P23) {
						p(x,y,z, 2,3);
					}
					else if(applied_normalization_isometries[t] == P24) {
						p(x,y,z, 2,4);
					}
					else if(applied_normalization_isometries[t] == P34) {
						p(x,y,z, 3,4);
					}
					else if(applied_normalization_isometries[t] == T1) {
						sigma_1(x,y,z);
					}
					else if(applied_normalization_isometries[t] == T2) {
						sigma_2(x,y,z);
					}
				}
				
			}
		}
	}
}

void next( ) {
    collision( );
    propagation( );
}

void display( ) {
    int x, y, z, i, t;

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
    glBegin(GL_POINTS);
    for(x = 0; x < LATTICE_WIDTH; x++) {
        for(y = 0; y < LATTICE_HEIGHT; y++) {
            for(z = 0; z < LATTICE_DEPTH; z++) {
                for(i = 0; i < 24; i++) {
                    t = lattice[z][y][x][i];
					if(t) {
						glColor4f(t, t, t, 0.2);
					}
					else {
						glColor4f(0.0, 0.0, 0.5, 0.05);
					}
                    glVertex3f((x + i)*0.01, (y + i)*0.01, (z + i)*0.01);
                }
            }
        }
    }
    glEnd( );

    glutSwapBuffers( );
}


void idle( ) {
	next( );
    glutPostRedisplay( );
}

void init_vars( ) {
	int x, y, z, i;
	for(x = 0; x < LATTICE_WIDTH; x++) {
		for(y = 0; y < LATTICE_HEIGHT; y++) {
			for(z = 0; z < LATTICE_DEPTH; z++) {
				for(i = 0; i < 24; i++) {
					lattice[z][x][y][i] = 0;
					tmp_lattice[z][x][y][i] = 0;
				}
			}
		}
	}

	
    for(x = 0; x < RANDOM_CA_WIDTH; x++) {
        tmp_random_ca_row[x] = 0;
        random_ca_row[x]     = 0;
    }
    tmp_random_ca_row[0] = 1;
	tmp_random_ca_row[RANDOM_CA_WIDTH/2] = 1;

	for(x = 0; x < LATTICE_WIDTH/3; x++) {
		for(y = 0; y < LATTICE_HEIGHT/3; y++) {
			for(z = 0; z < LATTICE_DEPTH/3; z++) {
				for(i = 0; i < 24; i++) {
					if(rand30(64564) == 0) {
						lattice[z][x][y][i] = 0;
					}
					else {
						lattice[z][x][y][i] = 1;
					}
				}
			}
		}
	}
}

void init_gl( ) {
    glClearColor(0.0, 0.0, 0.0, 0.0);

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glFrustum(0, 2*WINDOW_WIDTH, 0, 2*WINDOW_HEIGHT, -80.0, 80.0);

    glMatrixMode(GL_MODELVIEW);

    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
    glShadeModel(GL_SMOOTH);
    glPointSize(POINT_SIZE);
	
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA,GL_ONE);
    glEnable ( GL_ALPHA_TEST ) ;

}

int main(int argc, char* argv[]) {
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
    glutInitWindowSize(WINDOW_WIDTH, WINDOW_HEIGHT);
    glutCreateWindow("fhp");
    glEnable(GL_DEPTH_TEST);
    init_vars( );
    init_gl();
    glutDisplayFunc(display);
    glutIdleFunc(idle);
    glutMainLoop();
	return 0;
}