/* PARR -- counting PARR configurations   Copyright (C) 2015 Frans Faase

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

GNU General Public License:
   http://www.iwriteiam.nl/GNU.txt

Details:
   http://www.iwriteiam.nl/D1509.html#26

*/

#include <stdio.h>
#include <string.h>
#include <time.h>


#define MAX_N 5
#define MAX_M 4
#define MAX_POINTS MAX_M*MAX_N
#define MAX_LINES (MAX_M-1)*MAX_N + MAX_M*(MAX_N-1) + 2*(MAX_M-1)*(MAX_N-1)
#define MAX_DIAG (MAX_M-1)*(MAX_N-1)

//#define USE_DOUBLE
#ifdef USE_DOUBLE
typedef double Counter;
#define FORMAT "%.0lf"
#define DFORMAT "%*.0lf"
#else
typedef unsigned long long Counter;
#define FORMAT "%lld"
#define DFORMAT "%*lld"
#endif

void countWithTransitions(bool all_lines, bool no_cross, bool connected, bool verify, bool save_slices);


Counter counters[MAX_POINTS+1][MAX_LINES+1][MAX_DIAG+1];

void init_counters(int m, int n)
{
    const int points = n*m;
    const int lines = m*(n-1) + (m-1)*n + 2*(m-1)*(n-1);
    const int diags = (m-1)*(n-1);
    for (int p = 0; p <= points; p++)
        for (int l = 0; l <= lines; l++)
            for (int d = 0; d <= diags; d++)
                counters[p][l][d] = 0;
}

void print_counters_by_lines(int m, int n, int dia = -1)
{
    const int points = n*m;
    const int lines = m*(n-1) + (m-1)*n + 2*(m-1)*(n-1);
    const int diags = (m-1)*(n-1);
    int filled_lines = 0;
    int filled_points = 0;
    // determine formating of table
    int digits[MAX_POINTS+1];
    for (int p = 1; p <= points; p++)
        digits[p] = p < 10 ? 1 : 2;
    for (int l = 0; l <= lines; l++)
        for (int p = 1; p <= points; p++)
        {
            Counter counter = 0;
            if (dia == -1)
            {
                int f = 1;
                for (int d = 0; d <= diags; d++, f*=2)
                    counter += counters[p][l][d] * f;
            }
            else
                counter = counters[p][l][dia];
            if (counter != 0)
            {
                filled_lines = l+1;
                if (p > filled_points)
                    filled_points = p;
            }
            char buffer[30];
            sprintf(buffer, FORMAT, counter);
            int d = strlen(buffer);
            if (d > digits[p])
                digits[p] = d;
        }
    printf("<P>\n<PRE>\n    ");
    for (int p = 1; p <= filled_points; p++)
        printf(" %*d", digits[p], p);
    printf("\n   +");
    for (int p = 1; p <= filled_points; p++)
        for (int i = 0; i <= digits[p]; i++)
            printf("-");
    printf("\n");
    Counter total = counters[0][0][0];
    for (int l = 0; l < filled_lines; l++)
    {
        printf("%3d:", l);
        for (int p = 1; p <= filled_points; p++)
        {
            Counter counter = 0;
            if (dia == -1)
            {
                int f = 1;
                for (int d = 0; d <= diags; d++, f*=2)
                    counter += counters[p][l][d] * f;
            }
            else
                counter = counters[p][l][dia];
            printf(" "DFORMAT, digits[p], counter);
            total += counter;
        }
        printf("\n");
    }
    printf("</PRE>\n<P>\n");

    printf("The total number is  "FORMAT". <!-- (%ld)-->\n", total, (1L << 20));
}

long long over(int m, int n)
{
    long long result = 1;
    for (int i = 1; i <= n; i++, m--)
        result = result * m / i;
    return result;
}

void all(int m, int n, bool no_cross, const char* anchor)
{
    init_counters(m, n);

    const int points = n*m;
    
    for (long long c = 0; c < (1L << points); c++)
    {
        bool field[MAX_M][MAX_N];
        int k = 0;
        int d = 0;
        int p = 0;
        for (int i = 0; i < m; i++)
            for (int j = 0; j < n; j++, k++)
            {
                field[i][j] = (c & (1L << k)) != 0;
                if (field[i][j])
                    p++;
            }
        int l = 0;
        for (int i = 0; i < m; i++)
            for (int j = 0; j < n-1; j++, k++)
                if (field[i][j] && field[i][j+1])
                    l++;
        for (int i = 0; i < m-1; i++)
            for (int j = 0; j < n; j++, k++)
                if (field[i][j] && field[i+1][j])
                    l++;
        for (int i = 0; i < m-1; i++)
            for (int j = 0; j < n-1; j++, k++)
            {
                bool a = (field[i][j] && field[i+1][j+1]);
                if (a)
                    l++;
                bool b = (field[i][j+1] && field[i+1][j]);
                if (b)
                    l++;
                if (no_cross && a && b)
                {
                    l--;
                    d++;
                }
            }
        counters[p][l][d]++;    
    }
        
    print_counters_by_lines(m, n);
    
    countWithTransitions(true, no_cross, false, true, false);


    const int lines = m*(n-1) + (m-1)*n + 2*(m-1)*(n-1);
    const int diags = (m-1)*(n-1);
    for (int l = 0; l <= lines; l++)
        for (int p = 0; p <= points; p++)
            for (int l2 = l+1; l2 <= lines; l2++)
            {
                for (int d2 = 0; d2 <= diags; d2++)
                    for (int d = 0; d <= d2; d++)
                    {
                        int hvl = l - d;
                        int hvl2 = l2 - d2;
                        if (hvl >= 0 && hvl2 >= hvl)
                            counters[p][l][d] += counters[p][l2][d2] * over(d2, d) * over(hvl2, hvl);
                    }
            }

    printf("<P>\n<A NAME=\"%s\">\n", anchor);
    printf("And when we add all the PARR's where not all neighbour points are connected,\n");
    printf("we get the following table:\n");

    print_counters_by_lines(4, 5);

    countWithTransitions(false, no_cross, false, true, false);
}

bool trace = false;

#define WIDTH 4
#define MAX_LINES_WITH_WIDTH 4*WIDTH-3
#define VLINE(X) (X)
#define D1LINE(X) (WIDTH+(X))
#define D2LINE(X) (2*WIDTH-1+(X))
#define HLINE(X) (3*WIDTH-2+(X))

class Slice
{
public:
    Slice(int n_points, int n_lines, int n_min, int n_max, Slice* n_next)
      : points(n_points), lines(n_lines), min(n_min), max(n_max), next(n_next) {}
    ~Slice()
    {
        delete next;
    }  
    
    int points;
    int lines;
    int min;
    int max;
    Slice *next;
};

class State;

class Transition
{
public:
    Transition(State *n_to)
    {
        to = n_to;
        next = 0;
        slices = 0;
        for (int p = 0; p <= WIDTH; p++)
            for (int l = 0; l <= MAX_LINES_WITH_WIDTH; l++)
                counts[p][l] = 0L;
    }
    ~Transition()
    {
        delete slices;
        delete next;
    }
    long counts[WIDTH+1][MAX_LINES_WITH_WIDTH+1];
    State *to;
    Slice *slices;
    Transition *next;
};

class State
{
public:
    State(int *n_value, bool n_first)
    {
        for (int i = 0; i < WIDTH; i++)
            value[i] = n_value[i];
        first = n_first;
        transitions = 0;
        next = 0;
        for (int p = 0; p <= MAX_POINTS; p++)
            for (int l = 0; l <= MAX_LINES; l++)
                old_counters[p][l] = 0L;
        if (first)
            old_counters[0][0] = 1;
                
    }
    ~State()
    {
        delete next;
        delete transitions;
    }
    int value[WIDTH];
    bool first;
    Counter old_counters[MAX_POINTS+1][MAX_LINES+1];
    Counter new_counters[MAX_POINTS+1][MAX_LINES+1];
    Transition *transitions;
    State *next;
};

State *all_states = 0;

State* get_state(int *value, bool first)
{
    State** ref_state = &all_states;
    for (; *ref_state != 0; ref_state = &(*ref_state)->next)
    {
        bool equal = (*ref_state)->first == first;
        for (int i = 0; equal && i < WIDTH; i++)
            equal = (*ref_state)->value[i] == value[i];
        if (equal)
            return *ref_state;
    }
    if (trace) 
    {
        printf("add state: ");
        for (int i = 0; i < WIDTH; i++)
            printf("%2d", value[i]);
        printf("\n");
    }
    return *ref_state = new State(value, first);
}

Transition* get_transition(State *from, State* to)
{
    Transition** ref_trans = &from->transitions;
    for (; *ref_trans != 0; ref_trans = &(*ref_trans)->next)
        if ((*ref_trans)->to == to)
            return *ref_trans;
    if (trace) printf("add transition\n");
    return *ref_trans = new Transition(to);
}

void combine(int v1, int v2, int *old_value, int *new_value)
{
    if (v1 == 0 || v2 == 0)
        return;
    if (v2 < v1)
    {
        int h = v2;
        v2 = v1;
        v1 = h;
    }
    if (v1 < v2)
    {
        for (int i = 0; i < WIDTH; i++)
        {
            if (old_value[i] == v2)
                old_value[i] = v1;
            if (new_value[i] == v2)
                new_value[i] = v1;
        }
    }
}

void normalize(int *new_value)
{
    int mapping[4*WIDTH];
    for (int i = 0; i < 4*WIDTH; i++)
        mapping[i] = 0;
    int next_value = 1;
    for (int i = 0; i < WIDTH; i++)
    {
        int v = new_value[i];
        if (v > 0)
        {
            if (mapping[v] == 0)
                mapping[v] = next_value++;
            new_value[i] = mapping[v];
        }
    }
}

void countWithTransitions(bool all_lines, bool no_cross, bool connected, bool verify, bool save_slices)
{
    delete all_states;
    all_states = 0;

    // Fill initial element
    int value[WIDTH];
    for (int i = 0; i < WIDTH; i++)
        value[i] = 0;
    get_state(value, true);
    
    if (trace) printf("<!--\n");
    
    // expand states
    for (State *state = all_states; state != 0; state = state->next)
    {
        if (trace)
        {
            printf("Processing state:");
            for (int i = 0; i < WIDTH; i++)
                printf("%2d", state->value[i]);
            printf("\n");
        }
        bool final = !state->first;
        for (int i = 0; final && i < WIDTH; i++)
            final = state->value[i] == 0;
        for (int ps = 0; ps < (1 << WIDTH); ps++)
        {
            if (trace) printf(" points %d\n", ps);
            bool new_points[WIDTH];
            int nr_points = 0;
            for (int i = 0; i < WIDTH; i++)
            {
                new_points[i] = (ps & (1 << i)) != 0;
                if (new_points[i])
                    nr_points++;
            }
            int max_nr_lines = 0;
            for (int i = 0; i < WIDTH; i++)
                if (state->value[i] != 0 && new_points[i])
                    max_nr_lines++;
            for (int i = 0; i < WIDTH-1; i++)
            {
                bool diag1 = (state->value[i] != 0 && new_points[i+1]);
                bool diag2 = (state->value[i+1] != 0 && new_points[i]);
                if (diag1)
                    max_nr_lines++;
                if (diag2)
                    max_nr_lines++;
                if (no_cross && diag1 && diag2)
                    max_nr_lines--;
                    
                if (new_points[i] && new_points[i+1])
                    max_nr_lines++;
            }
            for (long ls = 0L; ls < (1L << MAX_LINES_WITH_WIDTH); ls++)
            {
                //fprintf("  lines %d\n", ls);
                bool new_lines[MAX_LINES_WITH_WIDTH];
                int nr_lines = 0;
                for (int i = 0; i < MAX_LINES_WITH_WIDTH; i++)
                {
                    new_lines[i] = (ls & (1L << i)) != 0;
                    if (new_lines[i])
                        nr_lines++;
                }
                bool okay = !all_lines || nr_lines == max_nr_lines;
                if (no_cross)
                    for (int i = 0; okay && i < WIDTH-1; i++)
                        okay = !new_lines[D1LINE(i)] || !new_lines[D2LINE(i)];
                for (int i = 0; okay && i < WIDTH; i++)
                    okay = !new_lines[VLINE(i)] || (state->value[i] != 0 && new_points[i]);
                for (int i = 0; okay && i < WIDTH-1; i++)
                    okay = !new_lines[D1LINE(i)] || (state->value[i] != 0 && new_points[i+1]);
                for (int i = 0; okay && i < WIDTH-1; i++)
                    okay = !new_lines[D2LINE(i)] || (state->value[i+1] != 0 && new_points[i]);
                for (int i = 0; okay && i < WIDTH-1; i++)
                    okay = !new_lines[HLINE(i)] || (new_points[i] && new_points[i+1]);
                if (okay)
                {
                    State *new_state = 0;
                    if (!connected)
                    {
                        int new_value1[WIDTH];
                        for (int i = 0; i < WIDTH; i++)
                            new_value1[i] = new_points[i] ? 1 : 0;
                        int new_value2[WIDTH];
                        for (int i = 0; i < WIDTH; i++)
                            new_value2[WIDTH-1-i] = new_points[i] ? 1 : 0;
                        int diff = 0;
                        for (int i = 0; diff == 0 && i < WIDTH; i++)
                            if (new_value1[i] < new_value2[i])
                                diff = -1;
                            else if (new_value1[i] > new_value2[i])
                                diff = 1;
                        new_state = get_state(diff < 0 ? new_value1 : new_value2, false);
                    }
                    else if (ps == 0)
                    {
                        for (int i = 0; okay && i < WIDTH; i++)
                            okay = state->value[i] < 2;
                        if (okay)
                        {
                            int new_value[WIDTH];
                            for (int i = 0; i < WIDTH; i++)
                                new_value[i] = 0;
                            new_state = get_state(new_value, nr_points == 0 && state->first);
                        }
                    }
                    else if (!final)
                    {
                        int old_value[WIDTH];
                        int new_value[WIDTH];
                        for (int i = 0; i < WIDTH; i++)
                        {
                            old_value[i] = state->value[i];
                            new_value[i] = new_points[i] != 0 ? 2*WIDTH+i : 0;
                        }
                        for (int i = 0; i < WIDTH; i++)
                            if (new_lines[VLINE(i)])
                                combine(old_value[i], new_value[i], old_value, new_value);
                        for (int i = 0; i < WIDTH-1; i++)
                            if (new_lines[D1LINE(i)])
                                combine(old_value[i], new_value[i+1], old_value, new_value);
                        for (int i = 0; i < WIDTH-1; i++)
                            if (new_lines[D2LINE(i)])
                                combine(old_value[i+1], new_value[i], old_value, new_value);
                        for (int i = 0; i < WIDTH-1; i++)
                            if (new_lines[HLINE(i)])
                                combine(new_value[i], new_value[i+1], old_value, new_value);
                        for (int i = 0; okay && i < WIDTH; i++)
                            if (old_value[i] != 0)
                            {
                                okay = false;
                                for (int j = 0; !okay && j < WIDTH; j++)
                                    okay = new_value[j]    == old_value[i];
                            }
                        if (okay)
                        {
                            int new_value1[WIDTH];
                            for (int i = 0; i < WIDTH; i++)
                                new_value1[i] = new_value[i];
                            normalize(new_value1);
                            new_state = get_state(new_value1, false);
                        }
                    }
                    if (new_state != 0)
                    {
                        Transition* new_transition = get_transition(state, new_state);
                        new_transition->counts[nr_points][nr_lines]++;
                    
                        if (save_slices)
                        {
                            int min = 100;
                            int max = 0;
                            for (int i = 0; i < WIDTH; i++)
                                if (new_points[i])
                                {
                                    max = i;
                                    if (i < min)
                                        min = i;
                                }
                            
                            new_transition->slices = new Slice(ps, ls, min, max, new_transition->slices);
                        }
                    }
                }
            }
        }
    }
    if (trace) printf("-->\n");
    
    for (int loop = 0; loop < 6; loop++)
    {
        for (State *state = all_states; state != 0; state = state->next)
            for (int p = 0; p <= MAX_POINTS; p ++)
                for (int l = 0; l <= MAX_LINES; l++)
                    state->new_counters[p][l] = 0L;
        for (State *state = all_states; state != 0; state = state->next)
            for (Transition *transition = state->transitions; transition != 0; transition = transition->next)
            {
                State *to_state = transition->to;
                for (int p = 0; p <= MAX_POINTS; p++)
                    for (int l = 0; l <= MAX_LINES; l++)
                    {
                        Counter old_counter = state->old_counters[p][l];
                        for (int pc = 0; pc <= WIDTH && p+pc <= MAX_POINTS; pc++)
                            for (int lc = 0; lc <= MAX_LINES_WITH_WIDTH && l+lc <= MAX_LINES; lc++)
                                to_state->new_counters[p+pc][l+lc] += transition->counts[pc][lc] * old_counter;
                    }
            }
        for (State *state = all_states; state != 0; state = state->next)
            for (int p = 0; p <= MAX_POINTS; p ++)
                for (int l = 0; l <= MAX_LINES; l++)
                    state->old_counters[p][l] = state->new_counters[p][l];
    }
    State* final_state = 0;
    for (State *state = all_states; state != 0; state = state->next)
    {
        bool final = !state->first;
        for (int i = 0; final && i < WIDTH; i++)
            final = state->value[i] == 0;
        if (final)
        {
            final_state = state;
            break;
        }
    }
    if (final_state == 0)
    {
        printf("No final state\n");
        return;
    }
    
    if (verify)
    {
        const int diags = (5-1)*(4-1);
        bool equal = true;
        for (int l = 0; equal && l <= MAX_LINES; l++)
            for (int p = 0; equal && p <= MAX_POINTS; p ++)
            {
                Counter counter = 0;
                int f = 1;
                for (int d = 0; d <= diags; d++, f*=2)
                    counter += counters[p][l][d] * f;
                equal = (counter == final_state->new_counters[p][l]);
            }
        if (equal)
        {
            printf("Result match with <A HREF=\"#sec\">second method</A>.\n");
            return;
        }
        else
            printf("Result does not match with second method. Result of second method is:\n");
            
    }

    init_counters(4, 5);
    for (int l = 0; l <= MAX_LINES; l++)
        for (int p = 0; p <= MAX_POINTS; p ++)
            counters[p][l][0] = final_state->new_counters[p][l];
    print_counters_by_lines(4, 5);
}

class OutputBitStream
{
public:
    OutputBitStream(const char *filename)
    {
        _f = fopen(filename, "wb");
        _bit = 1;
        _value = 0;
    }
    ~OutputBitStream() { flush(); fclose(_f); }
    void push(bool bit)
    {
        if (bit)
            _value |= _bit;
        _bit *= 2;
        if (_bit == 256)
            flush();
    }
    void push(int value, int nr_bits)
    {
        int v = 1;
        for (int i = 0; i < nr_bits; i++, v *= 2)
            push((value & v) != 0);
    }
    void flush()
    {
        if (_bit != 1)
        {
            _bit = 1;
            fwrite(&_value, 1, 1, _f);
            _value = 0;
        }
    }
        
private:
    FILE* _f;
    int _bit;
    int _value;
};

class Configuration
{
public:
    Configuration(Slice **slices, int depth, int width)
      : _slices(slices), _depth(depth), _width(width)
      {
          mode(0);
      }
    void mode(int i)
    {
        _i = i;
        _hm = (i & 1) != 0;
        _vm = (i & 2) != 0;
        _dm = (i & 4) != 0;
    }
    bool more() { return _i < ((_depth == _width) ? 8 : 4); }
    void next() { mode(_i+1); }
    bool point(int i, int j) { return !_dm ? _point(i, j) : _point(j, i); }
    bool horzlines(int i, int j) { return !_dm ? _horzlines(i, j) : _vertlines(j, i); }
    bool vertlines(int i, int j) { return !_dm ? _vertlines(i, j) : _horzlines(j, i); }
    bool d1lines(int i, int j)
    { return   !_dm 
             ? (_hm == _vm ? _d1lines(i, j) : _d2lines(i, j))
             : (_hm == _vm ? _d1lines(j, i) : _d2lines(j, i));
    }
    bool d2lines(int i, int j)
    { return   !_dm 
             ? (_hm == _vm ? _d2lines(i, j) : _d1lines(i, j))
             : (_hm == _vm ? _d2lines(j, i) : _d1lines(j, i));
    }
    
    void printData()
    {
        printf("%d, %d\n", _depth, _width);
        for (int i = 0; i < _depth; i++)
        {
            for (int j = 0; j < WIDTH; j++)
                printf(" %d", (_slices[i]->points >> j) & 1);
            printf("  ");
            for (int j = 0; j < 4*WIDTH-3; j++)
                printf(" %d", (_slices[i]->lines >> j) & 1);
            printf("  %5x, %5x\n", _slices[i]->points, _slices[i]->lines);
        }
        printf("\n");
    }
    void print(char *lines[], int offset)
    {
        for (int i = 0; i < 9; i++)
        {
            int len = strlen(lines[i]);
            for (int j = len; j < offset; j++)
                lines[i][j] = ' ';
            for (int j = 0; j < 7; j++)
                lines[i][offset+j] = (j % 2) == 0 && (i % 2) == 0 ? '_' : ' ';
            lines[i][offset + 7] = '\0';
        }
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width; j++)
                lines[i*2][offset+j*2] = point(i, j) ? '*' : '.';
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width-1; j++)
                if (horzlines(i, j))
                    lines[i*2][offset+j*2+1] = '-';
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width; j++)
                if (vertlines(i, j))
                    lines[i*2+1][offset+j*2] = '|';
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
                if (d1lines(i,j))
                    lines[i*2+1][offset+j*2+1] = d2lines(i,j) ? 'X' : '\\';
                else if (d2lines(i,j))
                    lines[i*2+1][offset+j*2+1] = '/';
    }
    int compare(Configuration &rhs)
    {
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width; j++)
            {
                bool lv = point(i, j);
                bool rv = rhs.point(i, j);
                if (lv && !rv) return 1;
                if (rv && !lv) return -1;
            }
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width-1; j++)
            {
                bool lv = horzlines(i, j);
                bool rv = rhs.horzlines(i, j);
                if (lv && !rv) return 1;
                if (rv && !lv) return -1;
            }
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width; j++)
            {
                bool lv = vertlines(i, j);
                bool rv = rhs.vertlines(i, j);
                if (lv && !rv) return 1;
                if (rv && !lv) return -1;
            }
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
            {
                bool lv = d1lines(i, j);
                bool rv = rhs.d1lines(i, j);
                if (lv && !rv) return 1;
                if (rv && !lv) return -1;
            }
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
            {
                bool lv = d2lines(i, j);
                bool rv = rhs.d2lines(i, j);
                if (lv && !rv) return 1;
                if (rv && !lv) return -1;
            }
        return 0;
    }
    bool crossing()
    {
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
                if (d1lines(i, j) && d2lines(i, j))
                    return true;
        return false;
    }
    bool neighconn()
    {
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width-1; j++)
                if (point(i, j) && point(i, j+1) && !horzlines(i, j))
                    return false;
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width; j++)
                if (point(i, j) && point(i+1, j) && !vertlines(i, j))
                    return false;
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
                if (   (point(i, j) && point(i+1, j+1) && !d1lines(i, j))
                    || (point(i, j+1) && point(i+1, j) && !d2lines(i, j)))
                    return false;
        return true;
    }
    bool neighconn_noncross()
    {
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width-1; j++)
                if (point(i, j) && point(i, j+1) && !horzlines(i, j))
                    return false;
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width; j++)
                if (point(i, j) && point(i+1, j) && !vertlines(i, j))
                    return false;
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
                if (   (   (point(i, j) && point(i+1, j+1)) 
                        || (point(i, j+1) && point(i+1, j)))
                    && !d1lines(i, j) && !d2lines(i, j))
                    return false;
        return true;
    }
    void write(OutputBitStream &ostream)
    {
        ostream.push(_depth-1, 3);
        ostream.push(_width-1, 2);
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width; j++)
                ostream.push(point(i, j));        
        for (int i = 0; i < _depth; i++)
            for (int j = 0; j < _width-1; j++)
                if (point(i, j) && point(i, j+1))
                    ostream.push(horzlines(i, j));
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width; j++)
                if (point(i, j) && point(i+1, j))
                     ostream.push(vertlines(i, j));
        for (int i = 0; i < _depth-1; i++)
            for (int j = 0; j < _width-1; j++)
            {
                if (point(i, j) && point(i+1, j+1))
                    ostream.push(d1lines(i, j));
                if (point(i, j+1) && point(i+1, j))
                    ostream.push(d2lines(i, j));
            }        
    }
    
                        
private:
    inline int _hi(int i) { return !_hm ? i : _depth-1-i; }
    inline int _him1(int i) { return !_hm ? i : _depth-2-i; }
    inline int _vj(int j) { return !_vm ? j : _width-1-j; }
    inline int _vjm1(int j) { return !_vm ? j : _width-2-j; }
    inline bool _bit(int v, int i) { return (v & (1 << i)) != 0; }
    inline bool _point(int i, int j) { return _bit(_slices[_hi(i)]->points, _vj(j)); }
    inline bool _horzlines(int i, int j) { return _bit(_slices[_hi(i)]->lines, HLINE(_vjm1(j))); }
    inline bool _vertlines(int i, int j) { return _bit(_slices[1+_him1(i)]->lines, VLINE(_vj(j))); }
    inline bool _d1lines(int i, int j) { return _bit(_slices[1+_him1(i)]->lines, D1LINE(_vjm1(j))); }
    inline bool _d2lines(int i, int j) { return _bit(_slices[1+_him1(i)]->lines, D2LINE(_vjm1(j))); }
    Slice **_slices;
    int _depth;
    int _width;
    int _i;
    bool _hm;
    bool _vm;
    bool _dm;
};

#define NO_RESTRICTION 100
int search_till_nr_points = NO_RESTRICTION;
int search_till_nr_lines = NO_RESTRICTION;

void expand(State* cur_state, Slice *sel_slices[], int depth, int nr_points, int nr_lines, int min, int max, OutputBitStream &ostream)
{
    for (Transition* transition = cur_state->transitions; transition != 0; transition = transition->next)
    {
        for (Slice* slice = transition->slices; slice != 0; slice = slice->next)
        {
            if (slice->points == 0)
            {
                if (!transition->to->first && min == 0)
                {
                    if (depth >= max+1)
                    {
                        bool smaller = false;
                        Configuration first(sel_slices, depth, max+1);
                        Configuration other(sel_slices, depth, max+1);
                        for (other.next(); other.more(); other.next())
                            if (first.compare(other) < 0)
                            {
                                smaller = true;
                                break;
                            }
                        if (!smaller)
                        {
                            counters[nr_points][nr_lines][0]++;
                            
                            first.write(ostream);
                            
                            bool crossing = first.crossing();
                            bool neighconn = first.neighconn();
                            bool neighconn_noncross = !crossing && first.neighconn_noncross();
                            if (!crossing)
                                counters[nr_points][nr_lines][1]++;
                            if (neighconn)
                                counters[nr_points][nr_lines][2]++;
                            if (neighconn_noncross)
                                counters[nr_points][nr_lines][3]++;
                        }
                    }    
                }
            }
            else if (depth < 5)
            {
                int slice_nr_points = 0;
                for (int i = 0; i < WIDTH; i++)
                    if ((slice->points & (1 << i)) != 0)
                        slice_nr_points++;
                int new_nr_points = nr_points + slice_nr_points;
                int slice_nr_lines = 0;
                for (int i = 0; i < MAX_LINES_WITH_WIDTH; i++)
                    if ((slice->lines & (1 << i)) != 0)
                        slice_nr_lines++;
                int new_nr_lines = nr_lines + slice_nr_lines;
                if (new_nr_points > search_till_nr_points || new_nr_lines > search_till_nr_lines)// || new_nr_lines > 5)
                {
                    // skip
                }
                else
                {
                    sel_slices[depth] = slice;
                    expand(transition->to, sel_slices, depth+1,
                           new_nr_points, 
                           new_nr_lines,
                           slice->min < min ? slice->min : min,
                           slice->max > max ? slice->max : max,
                           ostream);
                }
            }
        }
    }
}

void iterateUnique(State* all_con, const char *filename)
{
    OutputBitStream ostream(filename);
    clock_t tick = clock();
    init_counters(4, 5);
    Slice *slices[10];
    expand(all_con, slices, 0, 0, 0, 100, 0, ostream);
    printf("<!-- %lf min -->\n", (clock() - tick)/60.0/CLK_TCK);
}

int main(int argc, char *argv[])
{
    printf("<HTML><HEAD>\n");
    printf("<TITLE>Counting PARR configurations</TITLE>\n");
    printf("</HEAD><BODY BACKGROUND=\"greybgd.jpg\">\n");
    printf("\n");
    printf("<H1>Counting PARR configurations</H1>\n");
    printf("\n");
    printf("This page contains some results with respect to counting\n");
    printf("the number of <A HREF=\"Dpuzzle.html\">PARR's</A> with various properties. This page\n");
    printf("is generated with the program <A HREF=\"PARR_cpp.txt\">PARR.cpp</A>. This program went\n");
    printf("through a number of versions. Only the final version is released and likewise this page\n");
    printf("presents the results of the last version. The page is divided in the methods that were\n");
    printf("employed.\n");
    printf("\n");
    printf("<H2>First method</H2>\n");
    printf("\n");
    printf("The first method I applied was a straight forward counting algorithm. These were\n");
    printf("generated on <A HREF=\"D1508.html#26\">August 26, 2015</A>.\n");
    printf("\n");

    printf("<H3><A NAME=\"all\">All</A></H3>\n");
    printf("\n");
    printf("This table gives the number of all PARR's per number of points (columns)\n");
    printf("and per number connections (rows), including those with crossing diagonal connections,\n");
    printf("where all neighbour points are connected.\n");
    printf("\n");
    all(4, 5, false, "all_o");

    printf("<H3><A NAME=\"all_nc\">No crossings</A></H3>\n");
    printf("\n");
    printf("This table gives the number of all PARR's per number of points (columns)\n");
    printf("and per number of connections (rows) where all neighbouring points are connected,\n");
    printf("except for crossing diagonal connections where either one or the other is included.\n");
    all(4, 5, true, "all_nc_o");

    printf("\n");
    printf("<H2><A NAME=\"sec\">Second method</H2>\n");
    printf("\n");
    printf("The second method I applied involved slicing. These were\n");
    printf("generated on <A HREF=\"D1509.html#3\">September 3, 2015</A>.\n");
    printf("\n");
    
    printf("<H3><A NAME=\"all_con\">Connected</A></H3>\n");
    printf("\n");
    printf("This table gives the number of all PARR's per number of points (columns)\n");
    printf("and per number connections (rows), including those with crossing diagonal connections,\n");
    printf("where all neighbour points are connected and where all points are connected through\n");
    printf("one or more connections.\n");
    printf("\n");
    countWithTransitions(true, false, true, false, true);
    State *all_con = all_states;
    all_states = 0;
    printf("<P>\n<A NAME=\"all_con_o\">\n");
    printf("And when we add all the PARR's where not all neighbour points are connected,\n");
    printf("we get the following table:\n");
    countWithTransitions(false, false, true, false, true);
    State *all_con_o = all_states;
    all_states = 0;
    printf("\n");

    printf("<H3><A NAME=\"all_con_nc\">Connected and no crossings</A></H3>\n");
    printf("\n");
    printf("This table gives the number of all PARR's per number of points (columns)\n");
    printf("and per number of connections (rows) where all neighbouring points are connected,\n");
    printf("except for crossing diagonal connections where either one or the other is included.\n");
    printf("Futhermore, all points are connected through one or more connections.\n");
    countWithTransitions(true, true, true, false, true);
    State *all_con_nc = all_states;
    all_states = 0;
    printf("<P>\n<A NAME=\"all_con_nc_o\">\n");
    printf("And when we add all the PARR's where not all neighbour points are connected,\n");
    printf("we get the following table:\n");
    countWithTransitions(false, true, true, false, false);
    printf("\n");

    printf("<H2>Unique configurations</H2>\n");
    printf("\n");
    printf("The connected configurations still contain many configurations that only\n");
    printf("differ with respect to orientation and translation. I worked on a method\n");
    printf("for finding these. Below the tables for the various kinds of unique\n");
    printf("connected PARR configurations are given. (These were generated on\n");
    printf("<A HREF=\"D1509.html#26\">September 26, 2015</A>.) \n");

    search_till_nr_points = 9;
    iterateUnique(all_con_o, "PARR_con_o_9.data");
    printf("<P>\n");
    printf("<A NAME=\"all_uc_o\">\n");
    printf("All the unique PARR configurations (on %d or fewer points) without restrictions are:\n", search_till_nr_points);
    print_counters_by_lines(4, 5, 0);

    printf("<P>\n");
    printf("<A NAME=\"all_uc_nc_o\">\n");
    printf("The unique PARR configurations with no crossing diagonals are:\n");
    print_counters_by_lines(4, 5, 1);
    
    search_till_nr_points = NO_RESTRICTION;
    iterateUnique(all_con, "PARR_con.data");
    printf("<P>\n");
    printf("<A NAME=\"all_uc\">\n");
    printf("The unique PARR configurations where all neighbour points have a connection, are:\n");
    print_counters_by_lines(4, 5, 0);

    iterateUnique(all_con_nc, "PARR_con_nc.data");
    printf("<P>\n");
    printf("<A NAME=\"all_uc_nc\">\n");
    printf("The unique PARR configurations where all neighbour points have a connection,\n");
    printf("except for crossing diagonal connections where either one or the other is included, are:\n");
    print_counters_by_lines(4, 5, 0);

    printf("<P><HR>\n");
    printf("<ADDRESS>\n");
    printf("<A HREF=\"index.html\">Home</A>\n");
    printf("| <A HREF=\"Dpuzzle.html\">Puzzels</A>\n");
    printf("</ADDRESS>\n");
    printf("\n");
    printf("</BODY></HTML>\n");

    return 0;
}