/* Coded by Antti Karttunen.
   Some parts are quite generic, some parts are
   extremely kludgous.

   Compile with command like:
cc -I /dos/glow/traparot.dir/gd1_3/ -L /dos/glow/traparot.dir/gd1_3/ -g -o transham transham.c  -lgd -lm 

   You can start this with the arguments:
   ./transham 4 1234 1 2 1 2 1 > tulokset.out
   to get the 18 solutions beginning with the quick six through S3.

  CHANGES

    21.Feb.2000   karttu    Cleaned a bit. Added -p option for
                            finding palindromic Hamiltonian cycles.

 */

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

typedef unsigned long int Permutation;


int n_elems;
int n_vertices;
int elem_binwidth; /* In bits. */
int forbidden_elem = 0;
Permutation single_elem_mask; /* = (1<<(elem_binwidth))-1 */
Permutation double_elem_mask; /* = (1<<(2*elem_binwidth))-1 */
Permutation start_perm, end_perm;

int palindromes_flag = 0;
int every_other_transposition_must_be = 0;
int tot_solutions=0;

int *result_transp_array;
unsigned char *mark_array;

int bin_width(n)
{
    int j;
    for(j=0; n; j++, n>>=1) ;
    return(j);
}

int factorial(n)
{
    if(n < 2) { return(1); }
    else { return(n*factorial(n-1)); }
}

int adj_transpose(Permutation p,int transp_ind)
{
    int shift_count = (((n_elems-transp_ind)-1)*elem_binwidth);
    Permutation adj_pair = ((p >> shift_count) & double_elem_mask);

    adj_pair = (((adj_pair << elem_binwidth)
                 | (adj_pair >> elem_binwidth))
                & double_elem_mask);

    return((p & ~(double_elem_mask << shift_count))
            |
           (adj_pair << shift_count));
}

/* Is it possible to compute this easily
   from the other end? (for j from 1 to n do ...)
 */


int PermRevLexRank(Permutation p)
{
  int n,i,j,r = 0;
  Permutation a,b;

  for(j = n_elems; j; j--)
   {
      int b, a = (p&single_elem_mask);
      r += ((j-a)*(factorial(j-1)));
      p >>= elem_binwidth;
      for(i=0; (b=((p >> i)&single_elem_mask)); i += elem_binwidth)
       {
          if(b > a)
           {
             p = ((p & ~(single_elem_mask << i)) | ((b-1) << i));
           }
       }
   }
  return(r);
}

/* Algorithm 2.17: TrotterJonhsonRank(p) from
 * page 60 of Combinatorial Algorithms, Generation, Enumeration and Search,
 * Donald L. Kreher and Douglas R. Stinson. ISBN 0-8493-3988-X  CRC Press, 1998
 * See URL: http://www.math.mtu.edu/~kreher/cages.html
 */

int TrotterJohnsonRank(Permutation p)
{
    int i,j,k,r=0;
    Permutation b;

    for(j=2; j <= n_elems; j++)
     {
       k = 1;
       for(i=((n_elems-1)*elem_binwidth);
            ((b=((p >> i)&single_elem_mask)) != j);
           i -= elem_binwidth)
        {
          if(b < j) { k++; }
        }
       if(r&1) { r = ((j*r) + k - 1); }
       else    { r = ((j*r) + j - k); }
     }

    return(r);
}


void output_perm(FILE *outfp, Permutation p)
{
    if(!p) { return; }
    output_perm(outfp, (p >> elem_binwidth));
    putc(('0' + (p & single_elem_mask)),outfp);
}

/* **************************************************************** */

/* **************************************************************** */

#ifdef NOGIFS

struct IMPARA
{
    int dummy;
};

struct IMPARA *im_parameters = NULL;

int IM_draw_connection(Permutation prev_perm,Permutation next_perm,
                       int transposition,
                       struct IMPARA *im_pars)
{
}

void IM_output_it(struct IMPARA *im_pars, char *filename)
{
}

#else

/************************************************************************/
/*                                                                      */
/*               Image Related Functions Using GD 1.3                   */
/*                                                                      */
/************************************************************************/

#include "gd.h"
#include "gdfontt.h"

#define SIXTY_DEGREES 1.0471975511965976 /* In radians with gcl (acos 0.5) */

struct IMPARA
{
    gdImagePtr im; /* Image to which we draw on. */
    int linecolor; /* The connecting lines are drawn with this color. */
    int permcolor; /* The permutations 1234, 3241, etc. with this color. */
    int transcolor; /* The transposition numbers 1,2,3 or 4 with this. */
    int radius;
    int x_offset;
    int y_offset;
    int im_width;
    int im_height;
    int ptext_x_displ;
    int ptext_y_displ;
    char *filename_prefix;
    gdFontPtr font; /* set e.g. to gdFontSmall */
};

struct IMPARA *im_parameters = NULL;

void IM_initialize_image(struct IMPARA *im_pars)
{
    int white, black;

    /* Allocate the image: 64 pixels across by 64 pixels tall */
    im_pars->im = gdImageCreate(im_pars->im_width, im_pars->im_height);

    /* Allocate the color white (red, green and blue all maximum).
       Since this is the first color in a new image, it will
       be the background color.
     */
    white = gdImageColorAllocate(im_pars->im, 255, 255, 255);

    /* Allocate the color black (red, green and blue all minimum).
    black = gdImageColorAllocate(im_pars->im, 0, 0, 0);

    /*    im_pars->linecolor = black; */
    im_pars->linecolor = gdImageColorAllocate(im_pars->im,0,0,255); /* Blue. */
    im_pars->transcolor = gdImageColorAllocate(im_pars->im,255,0,0); /* Red. */
    im_pars->permcolor = gdImageColorAllocate(im_pars->im,0,255,0); /* Green */


}



struct IMPARA *IM_allocate_image(int radius,int x_offset,int y_offset, int im_width, int im_height,int ptext_x_displ,int ptext_y_displ,char *filename_prefix)
{
    struct IMPARA *im_pars;

    im_pars = ((struct IMPARA *) calloc(1,sizeof(struct IMPARA)));

    if(NULL == im_pars)
     {
       fprintf(stderr,"\nFailed: could not allocate struct IMPARA of size %u\n",
                sizeof(struct IMPARA));
       exit(1);
     }

    im_pars->im_width  = im_width;
    im_pars->im_height = im_height;
    im_pars->radius    = radius;
    im_pars->x_offset  = x_offset;
    im_pars->y_offset  = y_offset;
    im_pars->ptext_x_displ = ptext_x_displ;
    im_pars->ptext_y_displ = ptext_y_displ;
    im_pars->font      = gdFontTiny;
    im_pars->filename_prefix = filename_prefix;
    IM_initialize_image(im_pars);

    return(im_pars);
}





/* This works nice only with n_elems = 4, with 24 permutations. */
int get_x_and_y(int *px, int *py, Permutation p, int radius)
{
    int tj_rank = TrotterJohnsonRank(p);
    int distance, sector;

    distance = tj_rank % 8;
    if(tj_rank & 4) { distance = 7 - distance; }

    sector = tj_rank >> 2; /* tj_rank / 4, 0/4 - 23/4 -> 0 - 5 */

    *px = cos(+1*SIXTY_DEGREES*sector)*radius*(distance+1);
    *py = sin(+1*SIXTY_DEGREES*sector)*radius*(distance+1);
    return(tj_rank);
}

int IM_draw_connection(Permutation prev_perm,Permutation next_perm,
                       int transposition,
                       struct IMPARA *im_pars)
{
    int px1,py1,tj1rank;
    int px2,py2,tj2rank;
    char permbuf[81];

    tj1rank = get_x_and_y(&px1,&py1,prev_perm,im_pars->radius);
    tj2rank = get_x_and_y(&px2,&py2,next_perm,im_pars->radius);

/*
    fprintf(stderr,"Drawing line from %u,%u to %u,%u with color %d\n",
                     (px1+im_pars->x_offset), (py1+im_pars->y_offset),
                     (px2+im_pars->x_offset), (py2+im_pars->y_offset),
                     im_pars->linecolor);
    fflush(stderr);
 */

    gdImageLine(im_pars->im, (px1+im_pars->x_offset), (py1+im_pars->y_offset),
                             (px2+im_pars->x_offset), (py2+im_pars->y_offset),
                             im_pars->linecolor);

    if(01234 == prev_perm)
     {
       gdImageArc(im_pars->im,(px1+im_pars->x_offset),(py1+im_pars->y_offset),
                  4,4,0,360,im_pars->linecolor);

       sprintf(permbuf,"%04lo",prev_perm);
       gdImageString(im_pars->im, im_pars->font,
                     (px1+im_pars->x_offset+im_pars->ptext_x_displ),
                     (py1+im_pars->y_offset+im_pars->ptext_y_displ),
                     permbuf, im_pars->permcolor);
     }


    gdImageArc(im_pars->im, (px2+im_pars->x_offset), (py2+im_pars->y_offset),
               4,4,0,360,im_pars->linecolor);

    sprintf(permbuf,"%04lo",next_perm);
    gdImageString(im_pars->im, im_pars->font,
                 (px2+im_pars->x_offset+im_pars->ptext_x_displ),
                 (py2+im_pars->y_offset+im_pars->ptext_y_displ),
                  permbuf, im_pars->permcolor);

    gdImageChar(im_pars->im, im_pars->font,
                  (((px1+px2)/2)+im_pars->x_offset+3),
                  (((py1+py2)/2)+im_pars->y_offset),
                  (('a'-1)+transposition), /* Was: ('0'+transposition) */
                im_pars->transcolor);

}

void IM_output_it(struct IMPARA *im_pars, char *filename)
{
     /* Open a file for writing. "wb" means "write binary", important
        under MSDOS, harmless under Unix. */
     FILE *outfp = fopen(filename, "wb");
     if(NULL == outfp)
      {
        fprintf(stderr,"\nCould not open a file \"%s\" for writing!\n",
                        filename);
        exit(1);
      }

     /* Output the image to the disk file. */
     gdImageGif(im_pars->im, outfp);    

     /* Close the file. */
     fclose(outfp);

     /* Destroy the image in memory. */
     gdImageDestroy(im_pars->im);

     IM_initialize_image(im_pars); /* And allocate new image. */
}

#endif


/* **************************************************************** */

/* **************************************************************** */


int traverse(int level,int prev_trans,Permutation p,Permutation ep)
{
    int t;
    Permutation q,r,s,eq;

    if((level == n_vertices) /* Found one hamiltonian path... */
        || (palindromes_flag && (level > (n_vertices>>1)))
      )
     {
       int i;
       output_perm(stdout,start_perm); printf(" -> ");
       for(i=1; i < n_vertices; i++)
        {
          printf("%u",result_transp_array[i]);
          /* if(i < (n_vertices-1)) */ { putchar(' '); }
        }
       if(palindromes_flag) { printf("(central: )"); }
       else { printf("-> "); }
       output_perm(stdout,p); printf("\n");

       for(i=1, r = start_perm; i < n_vertices; i++)
        {
          s = adj_transpose(r,result_transp_array[i]);
          printf("%u (%u) ", PermRevLexRank(r),  TrotterJohnsonRank(r));
          if(im_parameters)
           {
             IM_draw_connection(r,s,result_transp_array[i],im_parameters);
           }
          r = s;
        }
       printf("%u (%u) \n", PermRevLexRank(r), TrotterJohnsonRank(r));

       fflush(stdout);
       tot_solutions++;
#ifndef NOGIFS
       if(im_parameters)
        {
          char filename[2055];
          sprintf(filename,"%s%d.gif",im_parameters->filename_prefix,tot_solutions);
          IM_output_it(im_parameters,filename);
        }
#endif
     }

    mark_array[p] = 1;
    if(palindromes_flag) { mark_array[ep] = 1; }

    for(t=((every_other_transposition_must_be && (level&1))
                  ? every_other_transposition_must_be : 1);
          t < ((every_other_transposition_must_be && (level&1))
                  ? (every_other_transposition_must_be+1) : n_elems);
        t++)
     {
       if(t != prev_trans) /* Do not try the same transposition twice in row */
        {
          q = adj_transpose(p,t);
          if(palindromes_flag)
           {
             eq = adj_transpose(ep,t);
           }
          if((q & single_elem_mask) != forbidden_elem) /* Forbidden permutations. (all in one specific coset, which fixes the last element of permutation as forbidden_elem) */
           {

/* The first condition is for the joining transposition of symmetric halves. */
             if((palindromes_flag && (p == eq) && (q == ep)
                  && (level == (n_vertices>>1)))
                ||
                 (!mark_array[q] /* Not encountered yet. */
                   && (!palindromes_flag || !mark_array[eq]))
               )
              {
                if(palindromes_flag)
                 {
                   result_transp_array[n_vertices-level] = t;
                 }
                result_transp_array[level] = t;
                traverse(level+1, t, q, eq);
              }
           }
        }
     }

    mark_array[p] = 0;
    if(palindromes_flag) { mark_array[ep] = 0; }
}

#ifdef NOGIFS

int draw_whole_perm4_graph()
{
}

#else

int draw_whole_perm4_graph()
{
    int t;
    Permutation p = 01234, next_p, q;
    int min_successor = 1;

    while(min_successor != -1)
     {
       int suc;
       min_successor = -1;
       for(t=1; t < n_elems; t++)
        {
          q = adj_transpose(p,t);
          if((suc = TrotterJohnsonRank(q)) > TrotterJohnsonRank(p))
           {
             if((-1 == min_successor) || (suc < min_successor))
              { min_successor = suc; next_p = q; }
             if(im_parameters)
              {
                IM_draw_connection(p,q,t,im_parameters);
              }
           }
        }
       p = next_p;
     }

    if(im_parameters)
     {
       char filename[2055];
       sprintf(filename,"%s%d.gif",im_parameters->filename_prefix,tot_solutions);
       IM_output_it(im_parameters,filename);
     }
}

#endif


main(int argc, char **argv)
{
    int perm_binwidth;
    int start_transpos = 1;
    int start_level = 0;
    int t_i,arg_i;
    char *p;
    Permutation per,eper;
    char *progname = *argv;

    while((argv[1]) && ('-' == (*(argv[1]))))
     {
       char c;

       switch(c = *(1+argv[1]))
        {
          case 'p': /* Find palindromic cycles. */
           {
             palindromes_flag++;
             break;
           }
          case '1': case '2': case '3': case '4': case '5':
          case '6': case '7': case '8': case '9':
           {
             every_other_transposition_must_be = (c - '0');
             break;
           }
          default:
           {
             fprintf(stderr,"Unknown option: %c\n",c);
             goto usage;
           }
        }
       argc--;
       argv++;
     }

    if(argc < 3)
     {
usage:
       fprintf(stderr,
"Usage: %s [-p] [-1] n starting_permutation [transpos1 [transpos2 [transpos3 ...]]]\n",
                progname);
       exit(1);
     }


    n_elems = atoi(*(argv+1));
    n_vertices = factorial(n_elems);

    elem_binwidth = bin_width(n_elems);
    perm_binwidth = n_elems*elem_binwidth;
    single_elem_mask = (1<<(elem_binwidth))-1;
    double_elem_mask = (1<<(2*elem_binwidth))-1;

    for(p = *(argv+2), start_perm = 0; *p; p++)
     {
       if(!isdigit(*p) || ('0' == *p) || (*p > ('0' + n_elems)))
        {
          fprintf(stderr,
             "%s: Only digits 1 - %u allowed in starting permutation!\n",
                     progname, n_elems);
          exit(1);
        }

       start_perm <<= elem_binwidth;
       start_perm += (*p - '0');
     }

    printf("n=%u, v=%u, ew=%u, pw=%u, se_mask=%lu, de_mask=%lu, start_perm=%lo\n",
              n_elems, n_vertices,
              elem_binwidth, perm_binwidth,
              single_elem_mask, double_elem_mask, start_perm);
    printf("Starting permutation: ");
    output_perm(stdout,start_perm);
    printf("\n");

#ifndef NOGIFS
    if((*((progname)+strlen(progname)-1)) == 'd')
     {
       int radius    = 60;
       int x_offset  = 250;
       int y_offset  = 250;
       int im_width  = 500+20;
       int im_height = 500;
       int ptext_x_displ = 5;
       int ptext_y_displ = 2;
       char *filename_prefix = "sol";
       im_parameters = IM_allocate_image(radius,x_offset,y_offset,im_width,im_height,ptext_x_displ,ptext_y_displ, filename_prefix);
     }
#endif


/* This was just for testing, commented out:
************************************************************
     {
       int t;
       Permutation q = start_perm;
       printf("After series of transpositions (1 2) - (%u %u): ",
               (n_elems-1),n_elems);
       for(t=1; t < n_elems; t++)
        {
          q = adj_transpose(q,t);
          output_perm(stdout,q);
          if(t < (n_elems-1)) { putchar(' '); }
        }
       printf("\n");
     }
************************************************************
 */

    if(NULL==(result_transp_array = ((int *) calloc(n_vertices,sizeof(int)))))
     {
       fprintf(stderr,"%s: Could not allocate result_transp_array of %u vertices\n",
                progname, n_vertices);
       exit(1);
     }

    if(NULL == (mark_array
      = ((unsigned char *) calloc((1<<perm_binwidth),sizeof(unsigned char)))))
     {
       fprintf(stderr,"%s: Could not allocate mark_array of %u bytes\n",
                progname, (1<<perm_binwidth));
       exit(1);
     }

    if(im_parameters) { draw_whole_perm4_graph(); }

    if(palindromes_flag)
     {
       eper = end_perm = adj_transpose(start_perm,(n_elems-1));
     }

    /* Start with command like: ./transham 4 1324 joo > tulokset.out */
    for(per = start_perm, t_i = 1, arg_i = t_i+2; argv[arg_i]; arg_i++, t_i++)
     {
       int t;
       t = atoi(argv[arg_i]);
       if((t < 1) || (t > (n_elems-1)))
        {
          fprintf(stderr,
             "%s: Only transpositions 1 - %u (i.e. (1 2) - (%u %u)) allowed!\n",
                     progname, (n_elems-1),(n_elems-1),n_elems);
          exit(1);
        }
       result_transp_array[t_i] = t;
       mark_array[per] = 1;
       per = adj_transpose(per,t);
       output_perm(stdout,per);

       if(palindromes_flag)
        {
          result_transp_array[n_vertices-t_i] = t;
          mark_array[eper] = 1;
          eper = adj_transpose(eper,t);
          putchar(' '); putchar('('); output_perm(stdout,eper); putchar(')');
        }

       if(argv[arg_i+1]) { putchar(' '); } else { putchar('\n'); }
     }

    traverse((argc-2),0,per,eper);


    printf("Total solutions=%u\n", tot_solutions);
}