#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>


/* Multiplication by a constant is much faster than not, so pick a
   power of 2.  */
#define MAX_N           8

/* The actual size of the square matrix to be manipulated.  */
static int N;

/* File names and descriptor for results output.  */
static FILE *outfile;
#define OUTFILENAME     "results.txt"


/* The state of the entire computation.  We'll be copying this around as we 
   prune potential submatricies, and constrain future permutations.  */

struct column_state
{
  /* The actual data of the column.  */
  unsigned char row[MAX_N];

  /* The number of non-constrained rows of the column.  Since we always
     constrain the diagonal, this number is never larger than N-1.  */
  unsigned char n_unconstrained;

  /* Counters used by the non-recursive permutation generator.  */
  unsigned char perm_count[MAX_N];

  /* Before beginning permutation on this column, ROW_MAP[N] is nonzero
     iff it has been constrained to a particular value.  After beginning
     permutation, the first N_UNCONSTRAINED entries of ROW_MAP point to
     the actual rows that are unconstrained.  This allows the permutation
     generator to quickly skip over the rows that are constrained.  */
  unsigned char row_map[MAX_N];
};

struct matrix_state
{
  struct column_state col[MAX_N];
};


/* Initialize M such that all columns are set to the canonical first
   permutation.  Make use of quandle axiom (i) and fix the diagonal
   such that it is not modified by the permutation generator.  */

static void
initialize_matrix (struct matrix_state *m)
{
  int i, j, n = N;

  memset (m, 0, sizeof (*m));

  for (i = 0; i < n; ++i)
    {
      for (j = 0; j < n; ++j)
        m->col[i].row[j] = j;

      m->col[i].n_unconstrained = n - 1;
      m->col[i].row_map[i] = 1;
    }
}


/* Add a constraint to a column for which we have yet to generate
   permutations.  Return true iff the new constraint does not conflict
   with any existing constraint.  */

static bool
constrain_permutation (struct column_state *col, int row, unsigned char val)
{
  if (col->row_map[row])
    return col->row[row] == val;
  else
    {
      col->row_map[row] = 1;
      col->n_unconstrained--;

      /* If the existing value in that row isn't the value we want, we
         need to find the value we want in another row, and swap out for
         the value that's currently in the row.  This step is required
         so that we preserve axiom (ii), that each column is a permutation
         and that each value appears exactly once.  */
      if (col->row[row] != val)
        {
          int i;
          unsigned char t;

          for (i = 0; ; ++i)
            if (col->row[i] == val)
              break;

          /* Verify that row I hasn't already been constrained.  */
          if (col->row_map[i])
            return false;

          t = col->row[row];
          col->row[row] = col->row[i];
          col->row[i] = t;
        }

      return true;
    }
}


/* Turn the boolean constrained/unconstrained markers in COL->ROW_MAP into
   a mapping to the N unconstrained rows.  */

static void
compute_row_map (struct column_state *col)
{
  int i, j;

  for (i = j = 0; i < N; ++i)
    if (col->row_map[i] == 0)
      col->row_map[j++] = i;

  for (; j < MAX_N; ++j)
    col->row_map[j] = -1;
}


/* Move to the next permutation using Heap's method.
   Sedgewick, Permutation Generation Methods, Computing Surveys, 1977.

   Return false iff there are no permutations remaining.

   Special case for our purposes: Skip elements that have been marked
   as constrained .  This means doing N_UNCONSTRAINED permutations around
   the remaining elements.  This allows us to automatically satisfy quandle
   axiom (i), and partially satisfy quandle axiom (iii).  */

static bool
next_permutation (struct column_state *col)
{
  int n;

  for (n = 0; n < col->n_unconstrained; n++)
    {
      unsigned char count = col->perm_count[n];

      if (count < n)
        {
          unsigned long x, y;
          unsigned char t;

          x = col->row_map[n & 1 ? count : 0];
          y = col->row_map[n];

          t = col->row[x];
          col->row[x] = col->row[y];
          col->row[y] = t;

          col->perm_count[n] = count + 1;
          return true;
        }
      else
        col->perm_count[n] = 0;
    }

  return false;
}


/* Validate quandle axiom (iii).  */

static bool
q3test (struct matrix_state *m)
{
  int i, j, k, n = N;

  for (i = 0; i < n; ++i)
    for (j = 0; j < n; ++j)
      for (k = 0; k < n; ++k)
        if (m->col[k].row[m->col[j].row[i]]
            != m->col[m->col[k].row[j]].row[m->col[k].row[i]])
          return false;

  return true;
}


/* Constrain a partially generated matrix via quandle axiom (iii).
   Examine columns up to COL and see if there is any information we
   can make use of.  Return false if we find that axiom (iii) has
   already been violated.  */

static bool
q3test_partial (struct matrix_state *m, int col)
{
  int i, j, k, n = N;

  for (k = 0; k < col; ++k)
    for (i = 0; i < n; ++i)
      {
        unsigned char Mki = m->col[k].row[i];

        for (j = 0; j < col; ++j)
          {
            unsigned char Mkj = m->col[k].row[j];
            unsigned char Mji = m->col[j].row[i];
            unsigned char MkMji = m->col[k].row[Mji];

            if (Mkj < col)
              {
                if (MkMji != m->col[Mkj].row[Mki])
                  return false;
              }
            else if (!constrain_permutation (&m->col[Mkj], Mki, MkMji))
              return false;
          }
      }

  return true;
}


/* Output a matrix in Maple format.  */

static void
output_matrix (struct matrix_state *m)
{
  char buf[256], *p = buf;
  int i, j, n = N;
  ssize_t buf_len;

  /* Construct the entire output in an internal buffer.  */

  strcpy (p, "matrix([");
  p = strchr (p, '\0');

  for (i = 0; i < n; ++i)
    {
      *p++ = '[';
      for (j = 0; j < n; ++j)
        {
          if (MAX_N < 9)
            *p++ = '0' + m->col[j].row[i] + 1;
          else
            p += sprintf (p, "%d", m->col[j].row[i] + 1);
          if (j < n - 1)
            *p++ = ',';
        }
      *p++ = ']';
      if (i < n - 1)
        *p++ = ',';
    }
  strcpy (p, "]),\n");
  buf_len = strlen (buf);

  /* The buffer size chosen should be easily large enough.  */
  if (buf_len > sizeof(buf))
    abort ();

  fwrite (buf, 1, buf_len, outfile);
}


/* Process all matrix permutations beginning with column COL.  Columns 
   0 through COL-1 have been initialized during program startup, and 
   correspond to the permutation state we read from the log file.  */

static void
quandleslist (struct matrix_state *m, int col)
{
  compute_row_map (&m->col[col]);

  if (col == N - 1)
    {
      do
        {
          if (q3test (m))
            output_matrix (m);
        }
      while (next_permutation (&m->col[col]));
    }
  else
    {
      do
        {
          struct matrix_state sub = *m;
          if (q3test_partial (&sub, col + 1))
            quandleslist (&sub, col + 1);
        }
      while (next_permutation (&m->col[col]));
    }
}


int
main(int ac, char **av)
{
  struct matrix_state init;

  if (ac != 2)
    {
      fprintf (stderr, "Usage %s <n>\n"
               "  N   is the size of the (square) matrix space to search\n",
               av[0]);
      return 1;
    }

  N = atoi (av[1]);

  if (N <= 0 || N > MAX_N)
    {
      fprintf (stderr, "N should be between 1 and %d.\n", MAX_N);
      return 1;
    }

  outfile = fopen (OUTFILENAME, "w");

  initialize_matrix (&init);
  quandleslist (&init, 0);

  fclose (outfile);
  return 0;
}