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

#define h 200               // height of the field
#define lw 4                // "long" width, i.e. number of unsigned long longs per row
int const w = 64+(lw-1)*60; // width of the field

/*

This program runs a multispin encoding based algorithm for muziks crystal rule but could be used for other automata
The idea was taken from https://doi.org/10.24355/dbbs.084-200511080100-675

The data is stored in unsigned long longs by overlapping consecutive elements for 4 bits. The first and last two bits
of each long are incorrect after one generation because of edge effects. By overlapping these mistakes are fixed

evol() implements the addition of 12 bits via the algorithm found on
https://en.wikipedia.org/wiki/Wallace_tree#Detailed_explanation

The output 4 bit number is described by (d_4<<3) + (d_3<<2) + (c_2<<1) + f1

This binary number from 0-12 is then evaluated with an additional bit (the cell in question) via the muziks crystal rule
The boolean expression f = c_2|f1; d_4&(l3&f|c_2&f1)|d_3&(l3|~f) was derived via the Quine McCluskey algorithm (then manually compressed)
https://github.com/int-main/Quine-McCluskey

Following terms were used:
minterms = [4, 11, 12, 20, 21, 22, 23, 25, 26, 27, 28]
dontcare = [13, 14, 15, 29, 30, 31]

This is derived from muziks crystal rule which is R2,C0,S4-7,9-12,B4,11-12,NB:
    
  c     0   1
n   0   
    1
    2
    3
    4   1   1
    5       1
    6       1
    7       1
    8
    9       1
    10      1
    11  1   1
    12  1   1
    13  x   x
    14  x   x
    15  x   x

x signifies a "don't care" since more than 12 neighbors is impossible, the cell is bit 5 (number = state<<4 + n)

On my computer this code takes for a random board with h = 200 and lw = 4 and iter = 1000 about 30 ms
The normal implementation takes about 360 ms, which is thus about a 10x speed increase

Currently the edges are set to zero but with some copying it's probably possible to implement toroidal geometry as well

*/

const int BYTES_PER_PIXEL = 3; /// red, green, & blue
const int FILE_HEADER_SIZE = 14;
const int INFO_HEADER_SIZE = 40;

unsigned long long evol(unsigned long long l1, unsigned long long l2, unsigned long long l3, unsigned long long l4, unsigned long long l5) {
    unsigned long long ia_0 = l1<<1;
    unsigned long long ia_1 = l1>>1;
    unsigned long long ia_2 = l2<<2;
    unsigned long long ia_3 = l2;
    unsigned long long ia_4 = l2>>2;
    unsigned long long ia_5 = l3<<1;
    unsigned long long ia_6 = l3>>1;
    unsigned long long ia_7 = l4<<2;
    unsigned long long ia_8 = l4;
    unsigned long long ia_9 = l4>>2;
    unsigned long long ia_10 = l5<<1;
    unsigned long long ia_11 = l5>>1;

    unsigned long long a1 = ia_0^ia_1^ia_2;
    unsigned long long b1 = ia_3^ia_4^ia_5;
    unsigned long long c1 = ia_6^ia_7^ia_8;
    unsigned long long d1 = ia_9^ia_10^ia_11;
    unsigned long long e1 = a1^b1^c1;
    unsigned long long f1 = e1^d1;
    
    unsigned long long a2 = (ia_0&ia_1)|(ia_2&(ia_0^ia_1));
    unsigned long long b2 = (ia_3&ia_4)|(ia_5&(ia_3^ia_4));
    unsigned long long c2 = (ia_6&ia_7)|(ia_8&(ia_6^ia_7));
    unsigned long long d2 = (ia_9&ia_10)|(ia_11&(ia_9^ia_10));
    unsigned long long e2 = (a1&b1)|(c1&(a1^b1));
    unsigned long long f2 = e1&d1;
    
    unsigned long long a_2 = a2^b2^c2;
    unsigned long long b_2 = d2^e2^f2;
    
    unsigned long long c_2 = a_2^b_2;
    unsigned long long c_3 = a_2&b_2;
    
    unsigned long long a_3 = (a2&b2)|(c2&(a2^b2));
    unsigned long long b_3 = (d2&e2)|(f2&(d2^e2));
    
    unsigned long long d_3 = a_3^b_3^c_3;
    unsigned long long d_4 = (a_3&b_3)|(c_3&(a_3^b_3));
    
    unsigned long long f = c_2|f1;
    return d_4&(l3&f|c_2&f1)|d_3&(l3|~f);
}

void set_board(bool f_[w][h]) {     // set a board with the usual array structure
    // fill the board randomly but set the edges to 0

    //srand(time(NULL));    // this line gives "true" random numbers
    srand(12);              // seed random number generator with this

    for (int y = 2; y < h-2; y++) {
        for (int x = 0; x < w; x++) {
            f_[x][y] = 0;//rand()&1;            // comment random back in later
        }
    }

    for (int y = 0; y < 2; y++) {
        for (int x = 0; x < w; x++) {
            f_[x][y] = 0;
            f_[x][h-y-1] = 0;
        }
    }

    for (int y = 2; y < h-2; y++) {
        for (int x = 0; x < 2; x++) {
            f_[x][y] = 0;
            f_[w-x-1][h] = 0;
        }
    }

    f_[w/2+2][h/2] = 1;         // 10 cell exploding seed
    f_[w/2+3][h/2] = 1;
    f_[w/2][h/2+1] = 1;
    f_[w/2+1][h/2+1] = 1;
    f_[w/2+2][h/2+1] = 1;
    f_[w/2][h/2+2] = 1;
    f_[w/2+1][h/2+2] = 1;
    f_[w/2+2][h/2+2] = 1;
    f_[w/2+1][h/2+3] = 1;
    f_[w/2+2][h/2+3] = 1;
}

void conv_to_llu(unsigned long long f[lw][h], bool f_[w][h]) {
    // converts normal 2d array to the uint64 version
    for (int y = 0; y < h; y++) {
        for (int x = 0; x < lw; x++) {
            f[x][y] = 0;
            for (int i = 0; i < 64; i++) {
                f[x][y] = (f[x][y]<<1)|f_[x*60 + i][y];
            }       
        }
    }
}

void p_llu(unsigned long long a) {
    // prints out binary
    for (int i = 63; i >= 0; i--) {
        printf("%i", (a>>i)&1);
    }
}

char bds_kernighan(unsigned long long a) {
    char cnt = 0;
    while (a != 0) {
        cnt++;
        a = a&(a-1);
    }
    return cnt;
}

int count_population(unsigned long long f[lw][h]) {     // this functions determines the population count
    int total = 0;
    unsigned long long mask1 = ~((unsigned long long) 3<<62);
    unsigned long long mask2 = ~((unsigned long long) 3);
    unsigned long long mask3 = mask1&mask2;

    for (int y = 0; y < h; y++) {
        total += bds_kernighan(f[0][y])&mask2;
        for (int x = 1; x < lw-1; x++) {
            total += bds_kernighan(f[x][y]&mask3);
        }
        total += bds_kernighan(f[lw-1][y]&mask1);
    }
    return total;
}

unsigned char* createBitmapFileHeader (int height, int stride) {
    int fileSize = FILE_HEADER_SIZE + INFO_HEADER_SIZE + (stride * height);

    static unsigned char fileHeader[] = {
        0,0,     /// signature
        0,0,0,0, /// image file size in bytes
        0,0,0,0, /// reserved
        0,0,0,0, /// start of pixel array
    };

    fileHeader[ 0] = (unsigned char)('B');
    fileHeader[ 1] = (unsigned char)('M');
    fileHeader[ 2] = (unsigned char)(fileSize      );
    fileHeader[ 3] = (unsigned char)(fileSize >>  8);
    fileHeader[ 4] = (unsigned char)(fileSize >> 16);
    fileHeader[ 5] = (unsigned char)(fileSize >> 24);
    fileHeader[10] = (unsigned char)(FILE_HEADER_SIZE + INFO_HEADER_SIZE);

    return fileHeader;
}

unsigned char* createBitmapInfoHeader (int height, int width) {
    static unsigned char infoHeader[] = {
        0,0,0,0, /// header size
        0,0,0,0, /// image width
        0,0,0,0, /// image height
        0,0,     /// number of color planes
        0,0,     /// bits per pixel
        0,0,0,0, /// compression
        0,0,0,0, /// image size
        0,0,0,0, /// horizontal resolution
        0,0,0,0, /// vertical resolution
        0,0,0,0, /// colors in color table
        0,0,0,0, /// important color count
    };

    infoHeader[ 0] = (unsigned char)(INFO_HEADER_SIZE);
    infoHeader[ 4] = (unsigned char)(width      );
    infoHeader[ 5] = (unsigned char)(width >>  8);
    infoHeader[ 6] = (unsigned char)(width >> 16);
    infoHeader[ 7] = (unsigned char)(width >> 24);
    infoHeader[ 8] = (unsigned char)(height      );
    infoHeader[ 9] = (unsigned char)(height >>  8);
    infoHeader[10] = (unsigned char)(height >> 16);
    infoHeader[11] = (unsigned char)(height >> 24);
    infoHeader[12] = (unsigned char)(1);
    infoHeader[14] = (unsigned char)(BYTES_PER_PIXEL*8);

    return infoHeader;
}

void generateBitmapImage (unsigned char* image, int height, int width, char* imageFileName) {
    int widthInBytes = width * BYTES_PER_PIXEL;

    unsigned char padding[3] = {0, 0, 0};
    int paddingSize = (4 - (widthInBytes) % 4) % 4;

    int stride = (widthInBytes) + paddingSize;

    FILE* imageFile = fopen(imageFileName, "wb");

    unsigned char* fileHeader = createBitmapFileHeader(height, stride);
    fwrite(fileHeader, 1, FILE_HEADER_SIZE, imageFile);

    unsigned char* infoHeader = createBitmapInfoHeader(height, width);
    fwrite(infoHeader, 1, INFO_HEADER_SIZE, imageFile);

    int i;
    for (i = 0; i < height; i++) {
        fwrite(image + (i*widthInBytes), BYTES_PER_PIXEL, width, imageFile);
        fwrite(padding, 1, paddingSize, imageFile);
    }

    fclose(imageFile);
}

int main() {
    int iter;

    // number of iterations of the board
    iter = 500;

    bool f_[w][h];
    set_board(f_);      // fill conventional board

    unsigned long long f[lw][h];
    unsigned long long f2[lw][h];

    conv_to_llu(f, f_);     // convert conventional board to custom format

    unsigned long long n15 = ~((unsigned long long) 15);        // bitmasks
    unsigned long long n15s = ~((unsigned long long) 15<<60);
    unsigned long long n3 = ~((unsigned long long) 3);
    unsigned long long n3s = ~((unsigned long long) 3<<62);

    printf("Simulation started!\n");
    clock_t begin = clock();

    for (int k = 0; k < iter; k++) {
        for (int y = 2; y < h-2; y++) {
            for (int x = 0; x < lw; x++) {
                f2[x][y] = evol(f[x][y-2], f[x][y-1], f[x][y], f[x][y+1], f[x][y+2]);

                // clean up edges between longs!
                if (x > 0) {
                    f2[x-1][y] = f2[x-1][y]&n3|((f2[x][y]>>60)&3);
                    f2[x][y] = f2[x][y]&n3s|((f2[x-1][y]&12)<<60);
                }
            }
            // set left and right edges to 00
            f2[0][y] = f2[0][y]&n3s;
            f2[lw-1][y] = f2[lw-1][y]&n3;
        }

        /*for (int y = 0; y < 2; y++) { // unneccessary if edges get set to zero at the beginning since these rows never change anyways
            for (int x = 0; x < 4; x++) {
                f2[x][y] = 0;
                f2[x][h-y-1] = 0;
            }
        }*/

        // copy the new field into the old field
        memcpy(f, f2, lw*h*sizeof(unsigned long long));
    }

    clock_t end = clock();
	double time_spent = (double)(end - begin) / CLOCKS_PER_SEC *1000;
    printf("That took %.0f ms\n", time_spent);


    /*printf("Saving image!\n");                      // this part of the code is for saving the image
    begin = clock();

    unsigned char image[h][w][BYTES_PER_PIXEL];
    
    char* imageFileName;
    char t[100];
    sprintf(t, "muzik_%i.bmp", iter);
    imageFileName = (char *) t;

    for (int y = 0; y < h; y++) {
        for (int i = 0; i < 62; i++) {
            image[h-y-1][i][2] = (unsigned char) ((f[0][y]>>(63-i)&1) * 255);      // red
            image[h-y-1][i][1] = (unsigned char) ((f[0][y]>>(63-i)&1) * 255);      // green
            image[h-y-1][i][0] = (unsigned char) ((f[0][y]>>(63-i)&1) * 255);      // blue
        }
        for (int x = 1; x < lw-1; x++) {
            for (int i = 2; i < 62; i++) {
                image[h-y-1][x*60+i][2] = (unsigned char) ((f[x][y]>>(63-i)&1) * 255);      // red
                image[h-y-1][x*60+i][1] = (unsigned char) ((f[x][y]>>(63-i)&1) * 255);      // green
                image[h-y-1][x*60+i][0] = (unsigned char) ((f[x][y]>>(63-i)&1) * 255);      // blue
            }
        }
        for (int i = 2; i < 64; i++) {
            image[h-y-1][(lw-1)*60+i][2] = (unsigned char) ((f[lw-1][y]>>(63-i)&1) * 255);      // red
            image[h-y-1][(lw-1)*60+i][1] = (unsigned char) ((f[lw-1][y]>>(63-i)&1) * 255);      // green
            image[h-y-1][(lw-1)*60+i][0] = (unsigned char) ((f[lw-1][y]>>(63-i)&1) * 255);      // blue
        }
    }

    generateBitmapImage((unsigned char*) image, h, w, imageFileName);

    end = clock();
    time_spent = (double)(end - begin) / CLOCKS_PER_SEC *1000;
    printf("That took %.0f ms\n", time_spent);*/

    //printf("%i\n", count_population(f));

	return 0;
}