module main;
import std.stdio, std.datetime, std.parallelism, std.algorithm, core.memory;

//Adjacency matrix version
//What can we cut from the matrix? First column - I think, first row will
//list all of the connections it makes - actually kept this to keep alignment
//with the sequence for bit hacking. Also the last row can be omitted, every
//connection to the last point will already have been listed. Finally every row can
//start listing points one later as the previous connections will already have been listed
//4.2 seconds SPEEDY!
//3409 matrixes with an unconstrained fold space, yet the constrained gives 3367 folds and
//the correct answer too   
enum LEN = 15;
enum OFFSET = 2;
enum DIM = LEN / 2 + (LEN & 1? 1 : 0);
enum SUBOPTIMAL = 2;
enum TARGET = 263_916;

struct coord { int y, x; }


pure int bitRev(int n)
{   int rev = 0;
    foreach(i;0..LEN)
        if(n & 1 << i)
            rev += 1 << LEN - 1 - i;
    return rev;
}


pure int bitCount(int i) //Slower than lookup
{   i = i - ((i >> 1) & 0x55555555);
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
    return (((i + (i >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
}


//Recursive shape search
void fold(  int position,     // current term of H/P
            ref int optimal,
            ref coord[] path,
            ref int[DIM][DIM] board,
            ref bool[ushort[LEN - 1]] unique_matrix,
            ref ushort[LEN - 1][] matrixes)   // coordinates of groups
{   board[path[position].y][path[position].x] = position + 1; //Block the current square
    if(position == 0) //Reached the end of the sequence
    {   int hh_count = 0;
        foreach(i, c;path)
        {   if(c.x > 0 && board[c.y][c.x - 1]) //Left is full
                ++hh_count;
            if(c.y > 0 && board[c.y - 1][c.x]) //Up is full
                ++hh_count;
        }

        //Check if this is the highest number of bonds arrangement
        if(hh_count > optimal)
            optimal = hh_count;
        if(hh_count >= optimal - SUBOPTIMAL) //Add suboptimal folds to the path list
        {   ushort[LEN - 1] matrix;
            //i iteration number, c coordinate. We don't need to check the last coord
            //as everything linking to it will already have been found
            foreach(i, c;path[0..$ - 1]) 
            {   if(c.x > 0 && board[c.y][c.x - 1] - 1 > i) //Left, test bound and that coord is later in path
                    matrix[i] ^= 1 << (board[c.y][c.x - 1] - 1);
                if(c.x < DIM - 1 && board[c.y][c.x + 1] - 1 > i) //Right
                    matrix[i] ^= 1 << (board[c.y][c.x + 1] - 1);
                if(c.y > 0 && board[c.y - 1][c.x] - 1 > i) //Up
                    matrix[i] ^= 1 << (board[c.y - 1][c.x] - 1);
                if(c.y < DIM - 1 && board[c.y + 1][c.x] - 1 > i) //Down
                    matrix[i] ^= 1 << (board[c.y + 1][c.x] - 1);
            }
            if(unique_matrix.get(matrix, 0) == 0)
            {   unique_matrix[matrix] = true;
                matrixes ~= matrix;
            } 
        }
    }
    else
    {   //Otherwise- try moving in all 4 possible directions
        void advanceProtein(int my, int mx)
        {   if(board[my][mx] == 0)
            {   path[position] = coord(my, mx);
                fold(position, optimal, path, board, unique_matrix, matrixes);
                board[my][mx] = 0;
            }
        }

        const int x = path[position].x;
        const int y = path[position].y;
        --position;
        if(x > 0) advanceProtein(y, x -  1);
        if(x < board.length - 1) advanceProtein(y, x + 1);
        if(y > 0) advanceProtein(y - 1, x);
        if(y < board.length - 1) advanceProtein(y + 1, x);
    }
}

pure int[] solve(const bool[] redundancy, const ushort[LEN - 1][] matrixes, const int[2^^LEN] bits, const int[] numbers)
{   int redun = true;
    static immutable int[] opti_table = [0,0,1,2,4,5,7,8,10,12,13,15,17,18,20,22];  
    int[] max; //Max H-H bonds seen for a given sequence
    max.length = 2^^LEN;
    foreach(sequence;numbers)
    {   int[] filled;
        foreach(i;0..LEN - 1) //-1 because the last adjacency row is not used
            if((sequence & 1 << i))
                filled ~= i;

        foreach(im, matrix;matrixes)
        {   if(redundancy[im]) 
                continue;
            int hh_count = 0;
            foreach(i;filled)
                //Counts the hydrogens present in the sequence which are neighbours to the row
                hh_count += bits[matrix[i] & sequence];

            //Check if this is the highest number of bonds arrangement
            if(hh_count > max[sequence])
            {   max[sequence] = hh_count;
                if(hh_count == opti_table[bits[sequence]]) //If this is the optimal number for n hydrogens break
                    break;
            }
        }
    }
    foreach(sequence;numbers)
        max[bitRev(sequence)] = max[sequence]; //Copy bit reversed optimals
    int temp = max.reduce!("a + b");
    delete max;
    //GC gc;
    //gc.collect;

    if(temp != TARGET) redun = false;
    return [redun, temp];
}


void main()
{   StopWatch sw;
    sw.start;
    
    coord[] path;
    path.length = LEN;
    int[DIM][DIM] board;
    bool[ushort[LEN - 1]] unique_matrix;
    ushort[LEN - 1][] matrixes;
    int optimal = 0;

    path[LEN - 1] = coord(OFFSET, OFFSET); //y, then x
    fold(LEN - 1, optimal, path, board, unique_matrix, matrixes);

    sw.peek.msecs.writeln("msecs after tree search");

    //Strict inferiority filter, removes folds that are strictly worse than others
    ushort[LEN - 1][] matrixes_filtered;
    bool[] superior;
    superior.length = matrixes.length;
    foreach(inum, i;matrixes.parallel)
    {   int j;
        for(j = 0;j < matrixes.length;++j)
        {   if(j == inum)
                continue;
            int k;
            for(k = 0;k < LEN - 1;++k)
                if(i[k] &(~matrixes[j][k]))
                    break;
            if(k == LEN - 1) //Got to the end with all matches, either a copy or inferior
                break;
        }
        if(j == matrixes.length) //Matrix had no matches nor strictly superior
            superior[inum] = true;
    }
    foreach(i, b;superior)
        if(b) matrixes_filtered ~= matrixes[i];

    matrixes = matrixes_filtered;
    
    //Also could keep path information and use that at the end? May be faster
    sw.peek.msecs.writeln("msecs after adjacency filtering");

    int[] numbers; //Remove reverse bits
    foreach(sequence;3..2^^LEN)
        if(sequence <= bitRev(sequence))
            numbers ~= sequence;

    int[2^^LEN] bits; //Set bits for bit counter
    bool[] redun; //Test for useless folds
    redun.length = matrixes.length;
    foreach(i;0..2^^LEN)
        bits[i] = bitCount(i); 

    enum SKIP = 0;
    foreach(i;0..SKIP)
        redun[i] = true;

    /*
    foreach(i, rr;redun)
    {   int[] temp = solve(redun, matrixes, bits, numbers);
        temp[1].writeln;
        //writefln("%.20f", (cast(double) temp[1]) / 2^^LEN);
        //redun[i] = cast(bool) temp[0];
    }
    */
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;
    solve(redun, matrixes, bits, numbers)[1].writeln;

    int used = 0;
    foreach(red;redun)
        if(red == 0) ++used;

    used.writeln;
    sw.peek.msecs.writeln("msecs at the end");
}