/*  zeta-orbits.c
**  
**  Frank Sottile    18 March 1998
**  
**  This file finds all irreducible, full-support permutations  in the 
**  symmetric group, up to the equivalence generated by cyclic shift, 
**  taking inverse,  and conjugaton by the longest element.
**     
**  It first finds representatives of this equivalence relation, but 
**  does not check for reducibility.  In this first part, it steps through
**  the symmetric group of order n in lexicographic order, beginning with the 
**  least derangement.   Then, if a permutation zeta satisfies
**    1)  zeta is a derangement,                     (fixpoint)
**    2)  k := #{ a | a<zeta(a) } \leq n/2,  and     (kbig)
**    3)  l := |zeta| \geq n-1                       (length)
**  it considers pi further.  We impose condition 2), as one of zeta and 
**  zeta^-1 satisfy this, and they are equivalent, and condition 3) is 
**  necessary for irreducibilty.
**  
**  Given such a zeta, it constructs all permutations equivalent to zeta
**  which also satisfy these conditions (there are at most 2n of these).
**  These are the cyclic shifts (orbit of conjugation by the cycle (12...n)
**  of  zeta and  w_0 zeta^-1 w_0).   Then it compares this set with all
**  representative previously encountered with the same k and l.  If a match 
**  is not found, it adds zeta to the list of representatives, computes the 
**  next lexicographic permutation, and repeats this step.
**
**  The second part of the program is dedicated to discarding the reducible 
**  permutations from this list.
**
**
**  VARIABLES
**  zeta        This is the current permutation in one-line notation
**          zeta[0]   = k := #{a | a< zeta(a) }
**          zeta[n+1] = length of zeta.
**  orbit       The equivalence orbit of zeta
**  perms       This 4-dimensional array stores equivalence classes 
**                of permutations sorted by k, length
**  frequency   This keeps track of numbers for each k, length
**  cycle[1+n +n/2];     /*  a permutation parsed into cycles 
**       spot:       0    |   1....n   | n+1.....n+n/2
**               # cycles | cycles  in |  lengths of 
**                        | lex order  |  cycles
**
**     example: zeta=(153)(472)(68) <----> 315324768332  
**
*/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define   n    7         /*  The index of the symmetric group */
#define   PermSize  n+2      /*  length of the permutation data */
#define   Numks     1+n/2     /*  number of possible k's */
#define   Nranks    2+n*n/4-n   /*  range of ranks */


int nextperm(int pi[PermSize]);  /*  finds the lexocographically next 
				permutation, and returns the last ascent. */ 
int fixpoint(int pi[PermSize]);  /*  checks for derangement.  */
int kbig(int pi[PermSize]);     /*  checks if k> n/2. */
int length(int pi[PermSize]);   /*  computes the length of zeta. */
int make_orbit(int zeta[PermSize],  int orbit[2*n+1][n+1]);
                         /*  Creates the equivalence orbit   */
int Pfrequency(int frequency[Numks][Nranks]);       
                            /*  Prints out a table of frequencies  */
int make_cycle(int pi[PermSize], int cycle[1+n+n/2]);
                            /* parses a permutation into cycles */
int Pcycle(int pi[PermSize], int cycle[1+n+n/2]);                    
                             /* Prints a permutation in cycle notation  */
int check_irred(int cycle[1+n+n/2]);
                            /* checks if irreducible */
int zeta[PermSize];
main()
{
 int i,ii, j,jj, asc, k, l, Ntest, in_there;
 int   orbit[2*n+1][n+1],  cycle[1+n+n/2],  pi[PermSize],
   perms[Numks][Nranks][60][n+1], frequency[Numks][Nranks],
   final_count[Numks][Nranks]; 

 /*  Initialize frequency   */ 
 for ( k=1; k <= n/2; ++k)
   for ( l=0; l <= Nranks-1; ++l){
     frequency[k][l]=0; 
     final_count[k][l]=0; 
   }

 /*  Initialize the permutation    */   
 zeta[n]=n;
 for ( i=1; i <= n/2; ++i) {
   zeta[2*i-1]=2*i;
   zeta[2*i]=2*i-1;
 }
asc=1;
  while(asc>0){

    if ( fixpoint(zeta) && kbig(zeta) && length(zeta)){

     make_orbit(zeta,orbit);

      k=zeta[0];
      l=zeta[n+1]+1-n;
      Ntest=frequency[k][l];
      in_there=0;
      if(Ntest>0){
      jj=Ntest-1;
      while( !in_there && (jj< 2*n*Ntest)){
	ii= (jj%(Ntest)) +1;
	j = 1+ jj/Ntest;
	
	in_there = (perms[k][l][ii][1]==orbit[j][1]);
	i=2;
	while (in_there&& i<=n){
	  in_there = in_there && 
	             (perms[k][l][ii][i]==orbit[j][i]);
	  ++i;
	  }
	++jj;
      }
      
      if (2*zeta[0]==n){
	
	for (i=1; i<=n; ++i)
	  pi[zeta[i]]=i;

	make_orbit(pi,orbit);
	k=zeta[0];
	l=zeta[n+1]+1-n;
	Ntest=frequency[k][l];
	jj=0;
	  while( !in_there && (jj< 2*n*Ntest)){
	    ii= (jj%(Ntest)) +1;
	    j = 1+ jj/Ntest;
	    
	    in_there = (perms[k][l][ii][1]==orbit[j][1]);
	    i=2;
	    while (in_there&& i<=n){
	      in_there = in_there && 
		(perms[k][l][ii][i]==orbit[j][i]);
	      ++i;
	    }
	    ++jj;
	  }
	}
      }

      if (!in_there) {
	frequency[k][l]=frequency[k][l]+1;
	for ( i=1; i <= n; ++i) 
	  perms[k][l][frequency[k][l]][i]=zeta[i];
      }
    }

    asc=nextperm(zeta);
   
    /*    if (asc==1) {
      Pfrequency(frequency);
      printf(" asc=1  zeta=");
	for (i=1; i<=n; ++i)
	  printf("%d",zeta[i]);
      printf("\n");
      }*/

  }

  /*   Print frequency table*/
  Pfrequency(frequency);

  for (k=1; k<=n/2; ++k){
    for (l=0; l<=Nranks-1; ++l){
      for (ii=1; ii<=frequency[k][l]; ++ii){ 
 
	  for (i=1; i<=n; ++i)
	    zeta[i]=perms[k][l][ii][i];

	  make_cycle(zeta,cycle);

	  if ((cycle[0]==1) || (check_irred(cycle))) 
	    final_count[k][l]=final_count[k][l]+1;

      }
    }
  }

  
  /*   Print final frequency table*/
  Pfrequency(final_count);


  printf("  %d\n",clock());

	      }
/*check_irred
**              checks if the permutation is irreducible 
*/
int check_irred(int cycle[1+n+n/2])
{
  int a, b, btest, i,ii, j, jj, crossing, place, Nmove, not_done;
  int Blocks[1+n+n/2], temp[1+n+n/2];

  not_done=1;

  for (i=0; i<= n+n/2; ++i)
    Blocks[i]=cycle[i];


  btest=2;
  place=Blocks[n+1];

  while (not_done){
    
    crossing=0;
    
    /* Checks to see if the first and btest-th block are crossing */
    for (i=1; i < Blocks[n+1]; ++i){
      for (j=i+1; j<=Blocks[n+1]; ++j){
	a=(Blocks[i]+Blocks[j]- abs(Blocks[i]-Blocks[j]));
	b=(Blocks[i]+Blocks[j]+ abs(Blocks[i]-Blocks[j]));
	
	for (ii=place+1; ii<place+Blocks[n+btest]; ++ii){
	  for (jj=ii+1; jj<=place+Blocks[n+btest]; ++jj){
	    crossing = crossing ||  
	      ( a < Blocks[ii]+Blocks[jj] - abs(Blocks[ii]-Blocks[jj]) ) &&
	      ( Blocks[ii]+Blocks[jj] - abs(Blocks[ii]-Blocks[jj]) <b) &&
	      ( b < Blocks[ii]+Blocks[jj] + abs(Blocks[ii]-Blocks[jj]));
	  }}
      }
    }
    
    if (!(crossing)){
      if (btest<Blocks[0]){
	place=place+Blocks[n+btest];
	++btest;}
      else
	return crossing;
    }
    
    if (crossing){
      if ((btest==2)&&(2==Blocks[0]))
	return crossing;
      else{
	/*         Merge 1st and btest-th Blocks  */
	for (i=0; i<= n+n/2; ++i)
	  temp[i]=Blocks[i];
	
	Blocks[0]=Blocks[0]-1;
	Blocks[n+1]=Blocks[n+1]+Blocks[n+btest];
	for (j=btest; j<n/2; ++j)
	  Blocks[n+j]=Blocks[n+j+1];
	
	if (btest>2){
	  Nmove=0;
	  for (j=2; j<btest; ++j)
	    Nmove=Nmove+temp[n+j];
	  
	  for (j=1; j<=temp[n+btest]; ++j)
	    Blocks[temp[n+1]+j]=temp[temp[n+1]+Nmove+j];
	  
	  for (j=1; j<=Nmove; ++j)
	    Blocks[Blocks[n+1]+j]=temp[temp[n+1]+j];
	  
	}
	btest=2;
	place=Blocks[n+1];
	
	crossing=0;
	
      }
    }
  }
  
}



/*Pcycle                    ** Prints a permutation in cycle notation 
*/
int Pcycle(int pi[PermSize], int cycle[1+n+n/2])                   
{
  int place, i, ii;
  place=1;
  for (i=1; i<=cycle[0]; ++i){
    printf(" (");
    for (ii=1; ii<=cycle[n+i]; ++ii){
      printf(" %d",cycle[place]);
      ++place;
    }
    printf(" ) ");
  }
/*  printf(" \n");*/
}

/*make_cycle
** This parses a permutation into cycles
*/
int make_cycle(int pi[PermSize], int cycle[1+n+n/2])
{
  int i, ncycles, begin, length, counter, test;

/*  Initialize cycle */
  for ( i=0; i <= n+n/2; ++i)
   cycle[i]=0;


  ncycles=1;
  begin=1;
  length=1;
  counter=1;

  while(counter<n){
    cycle[counter]=begin;

    while (!(pi[cycle[counter]]==begin)){
      cycle[counter+1]=pi[cycle[counter]];
      ++counter;
      ++length;
    }
    cycle[n+ncycles]=length;
  if (counter<n){
      test=2;
      ++ncycles;
      ++begin;
     while ( test <= counter ) {
	if (begin == cycle[test] ){
	  ++begin;
	  test=2;
	}
	else
	  ++test;
      }
  }
  ++counter;
  length=1;

  }

  cycle[0]=ncycles;
}


/*Pfrequency
**  This prints out a table of frequencies
*/
int Pfrequency(int frequency[Numks][Nranks])
{
  int k,l;
  printf("\n");
  printf("Table of Frequencies\n");
  printf("\\ k");
  for (k=1; k<=n/2; ++k)
    printf("  %d  ",k);
  printf("\n");
  printf("l\\_________________\n");
  
  for (l=0; l<= Nranks-1; ++l){
    printf("%d | ",l+n-1);
    for (k=1; k<=n/2; ++k)
      printf("  %d ",frequency[k][l]);
    printf("\n");
  }
  printf("\n");
  
}


/*make_orbit(zeta,orbit)
**   This makes the  equivalence orbit
*/

int make_orbit(int zeta[PermSize],  int orbit[2*n+1][n+1])
{
  int i,j;
  for ( i=1; i <= n; ++i) {
    orbit[1][i]=zeta[i];
    orbit[n+1][n+1-zeta[i]]=n+1-i;
  }
  
  for ( j=2; j <= n; ++j) {
    for ( i=1; i <= n; ++i) {
      orbit[j][(i%n)+1]=(orbit[j-1][i]%n)+1; 
      orbit[n+j][(i%n)+1]=(orbit[n+j-1][i]%n)+1;
    }
  }
}


/*length(k,zeta)
**  This computes the length of zeta
*/

int length(int pi[PermSize])
{
  int i,rank,a,b;
  rank=0;
  for ( a=2; a<=n; ++a ){
    for ( b=1; b<a; ++b ){    
      if ( a<pi[a] ) {
        if ( b<pi[b] && pi[a]<pi[b]){
          rank = rank -1;
        }
        if ( b>pi[b]){
          rank = rank -1;
          if (pi[a]>pi[b]){
            rank = rank + 1;
	  }
        }
      }
      if ( a>pi[a] ) {
        if ( b<pi[b] && pi[a]<pi[b]){
          rank = rank +1;
        }
        if ( b>pi[b] && pi[a]<pi[b]){
          rank = rank -1;
        }
      }
    }
  }
pi[n+1]=rank;
return (rank >= n-1);
}

/*kbig(k,zeta)
**  This checks to see if k<=n/2*/

int kbig(int pi[PermSize])
{
  int i,k;
  k=0;
  for ( i=1; i <= n; ++i) 
    if ( i<pi[i]) k=k+1;
   pi[0]=k;
  return(k<= n/2);
}


/*fixpoint(zeta)
**  This checks to see if the permutation zeta has a fixed point
*/
int fixpoint(int pi[PermSize])
{
  int i,fp;
  fp=1;
  for ( i = 2; i<=n; ++i)
    fp = (fp && !(pi[i]==i));
  return fp;
}


/*nextperm(pi):
**  This finds the next permutation in the lexicographic order and it
**  returns the value of the last ascent of pi, to ensure the program 
**  halts.
*/
int nextperm(int pi[PermSize])
{
int asc, i, temp, bigger;

asc = 0;
 
for ( i = 1; i<n; ++i)
  if ( pi[i] < pi[i+1] )
    asc = i;

if ( asc > 0 ) {
  temp = pi[asc];

  for ( i = asc+1; i <= n; ++i){
    if ( pi[i] > temp )
      bigger = i;
  }
  pi[asc]=pi[bigger];
  pi[bigger]=temp;
  for ( i=1; i <= (n-asc)/2; ++i){
    temp=pi[asc+i];
    pi[asc+i]=pi[n+1-i];
    pi[n+1-i]=temp;
  }
}

return asc;

}