// Calling sequence
//   
//  ./C11-2-perm1 < C11-2-BPermFile
//

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>

#include "../numeric/SU21-numeric.h"
#include "../numeric/SU21-numeric-base.c"
#include "../permsearch/SU21-permsearch-base.c"

// This program deals with permuatations~$B$ which act on the group~$K$
// of cardinality~288 and satisfy 10~relations of the forms:
//
//    B k_1 = k_2 B  and B k_3 = k_4 B^-1
//
// for various given pairs of elements~$(k_1,k_2)$ and $(k_3,k_4)$
// of~$K$, where the~$k_i$ act on~$K$ by left-multiplication.
//
// The set on which $B$~acts, which is a copy of the group~$K$, enters
// into the conditions only through its structure as a left-$K$-set.
// Consequently any permutation of~$K$ which preserves its structure
// as a left $K$-set will conjugate any valid~$B$ to another such.  It
// is obvious (and very well-known) that a permutation preserves the
// structure of~$K$ as a left-$K$-set if and only if it is
// right-multiplication by some~$k\in K$.

// Elements~0 through~143 of~$K$ are of the form:
//
//   K1^l K2 K1^m K2
//
// while elements~144 through~287 are of the form:
//
//   K1^l K2 K1^m

// This program checks, for each permutation, in how many ways 
// conjugation by a right translation will give us a conjugated
// permutation satisfying:
//
//   (R_k B R_k^{-1}) (1) = B(1)
//


#define KSIZ 288
#define MaxConjCnt 120

#define MaxNumSigs 2000
#define NumSpecial 6
#define NumBCnt 6

void ShowSigTbl(int NumSigs,INT_t SigTbl[MaxNumSigs][NumSpecial],
		int SigCnt[MaxNumSigs]);
void yVal0Finish(int yVal0Ix,int yVal0,int BCnt[],int TotBCnt[],
		 int *SigCnt,int *TotSigCnt,
		 int NumSigs,INT_t SigTbl[MaxNumSigs][NumSpecial]);
inline int TrivialPower(int KSiz,int n,
			INT_t P1[],INT_t P2[],INT_t P3[],INT_t P4[]);
inline int nCycle(int KSiz,int n,
		  INT_t P1[],INT_t P2[],INT_t P3[],INT_t P4[]);


int main(int argc,char *(argv[])) {
  // Information about the group $K$
  int KSizFromFile,*KMult,*KInv;
  Elt_t *KElts;

  INT_t KMult2[KSIZ][KSIZ],KInv2[KSIZ];

  // Table which shows for each pair, $(x,y)$, for which it hold
  // that $B(x)=y$, what pair $(1,y')$ this can be conjugated to by
  // right-multiplication, i.e.
  //
  //   (R_x^{-1} B R_x) (1) = (R_x^{1} B) (x) = R_x^{-1} (y) = y x^{-1}
  //
  // This table gives its answer in the form of an index to the yVal
  // table.
  INT_t xyTbl[KSIZ][KSIZ],xyVal;
  // Similar table, showing whether yx^-1 is one of the ``special''
  // values used to calculate the signature.
  char xySpecialTbl[KSIZ][KSIZ],xySpecial;

  // Input permutation B, and its conjugates by right translations
  INT_t B[KSIZ],BConj[MaxConjCnt][KSIZ],BTest[KSIZ];

  // Input files
  FILE *KFile,*yValFile;

  char FileNam[MAX_FILENAME_LEN]="../";
  char Prefix[16]="C11/C11-2-shift",*Ix0_p;

  // Special values for BTest[0]
  int SpecialList[NumSpecial]={144,145,148,149,150,151};
  int Special,SpecialIx;
  // The ``signature'' of a permutation B is a tuple of counts, one
  // for each element of SpecialList.  If Special is one of the values
  // in SpecialList, the associated count is the number of elements x0
  // of K so that conjugating B by right-translation by x0 leads to a 
  // permutation BTest satisfying BTest[0]=Special.

  // Table of signatures
  INT_t SigTbl[MaxNumSigs][NumSpecial];
  int NumSigs=0,SigIx;
  int SigCnt[MaxNumSigs],TotSigCnt[MaxNumSigs];
  // Signature  of B
  INT_t BSig[NumSpecial];

  // List of y-values which were used for imposing $B(0)=y$ or excluding
  // $B(x)=yx$.
  int NumYVals,yValIx,yVal0Ix,yVal0IxOld,yVal0,yVal0Old,yVal;
  INT_t yVals[KSIZ];

  int x,y,x0,x0Inv;
  int MinFg,Comparison,FoundFg,ConjIx,ConjCnt;
  int BCnt[NumBCnt];
  int TotBCnt[NumBCnt];

  // Names and permutations for left-multiplication by certain
  // elements of~$K$
#define NumKNames 9
  char KNames[NumKNames][10]=
    {"1","K1","K2","K2K1K2","K1^4K2","K1^5K2","K1^6K2","K1^3K2",
     "K2K1^4K2"};
  INT_t KPerms[NumKNames][KSIZ];
#define Perm1      KPerms[0]
#define PermK1     KPerms[1]
#define PermK2     KPerms[2]
#define PermK2K1K2 KPerms[3]
#define PermK1u4K2 KPerms[4]
#define PermK1u5K2 KPerms[5]
#define PermK1u6K2 KPerms[6]
#define PermK1u3K2 KPerms[7]
#define PermK2K1u4K2 KPerms[8]

  // Open the files
  {
    size_t len=strlen(Prefix);
    assert(3+len+6<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix);
    // Input files
    strcpy(&FileNam[3]+len,"-K");
    printf("%s\n",FileNam);
    KFile=fopen(FileNam,"r");
    assert(KFile != NULL);

    strcpy(&FileNam[3]+len,"-yVal");
    printf("%s\n",FileNam);
    yValFile=fopen(FileNam,"r");
    assert(yValFile != NULL);
    fflush(stdout);
  }

  // Get the elements of~$K$
  GetKEltsKK(KFile,&KSizFromFile,&KElts,&KMult,&KInv);
  assert(KSizFromFile == KSIZ);
  assert(strcmp(KElts[0].Word,"") == 0);
  strcpy(KElts[0].Word,"1");
  assert(strcmp(KElts[144].Word,"K2") == 0);
  printf("GetKEltsKK, KSiz=%i\n",KSizFromFile);
  
  // Construct permutations corresponding to left-multiplication by
  // certain elements of~$K$
  {
    int Ix;
    for(Ix=0;Ix<NumKNames;Ix++)
      assert(NamToPerm(KPerms[Ix],KNames[Ix],KSIZ,KElts,KMult));
  }

  // Get the list of yVals
  GetYVals(yValFile,KSIZ,&NumYVals,yVals);
  printf("GetYVals, NumYVals=%i\n",NumYVals);

  // Copy the multiplication and inverse tables for~$K$ into a more
  // convenient format
  for(x=0;x<KSIZ;x++) {
    KInv2[x]=KInv[x];
    for(y=0;y<KSIZ;y++) 
      KMult2[x][y]=KMult[KSIZ*x+y];
  }
  free(KInv);
  free(KMult);

  memset(xyTbl,-1,sizeof(xyTbl));
  for(yValIx=0;yValIx<NumYVals;yValIx++) {
    yVal=yVals[yValIx];
    for(x=0;x<KSIZ;x++) {
      // y = k_{yVal} x,   yx^{-1}=k_{yVal}
      y=KMult2[yVal][x];
      xyTbl[x][y]=yValIx;
    }
  }

  memset(xySpecialTbl,-1,sizeof(xySpecialTbl));
  for(SpecialIx=0;SpecialIx<NumSpecial;SpecialIx++) {
    Special=SpecialList[SpecialIx];
    for(x=0;x<KSIZ;x++) {
      // y = k_{Special} x,   yx^{-1}=k_{Special}
      y=KMult2[Special][x];
      xySpecialTbl[x][y]=SpecialIx;
    }
  }

  // BCnt[0] is the number of permutations B read from stdin.
 
  // BCnt[1] and BCnt[2] are increased only if the input permutation B
  // is smallest among all conjugations, BTest, by
  // right-multiplications which satisfy BTest[0]=B[0].  BCnt[2]
  // counts the number of input permutations satisfying this
  // condition; BCnt[1] counts, for each input permutation satisfying
  // this condition, the number of different conjugations BTest by
  // right-multiplications such that BTest[0]=B[0].

  // BCnt[3] through BCnt[5] refer to permutations which satisfy
  // further relations.

  // The Total values refer to the entire run; the plain values refer
  // to a single value of B[0].

  memset(BCnt,0,sizeof(BCnt));
  memset(TotBCnt,0,sizeof(TotBCnt));

  // SigCnt counts the number of B's, excluding those which fail the
  // conjugate test, which have various signatures.
  memset(SigCnt,0,sizeof(SigCnt));
  memset(TotSigCnt,0,sizeof(TotSigCnt));

  yVal0IxOld=(-1);

  while(scanf("%*[ \n\t]") != EOF) {
    if(BCnt[0] % 10000 == 0 && BCnt[0] != 0) {
      //  if(BCnt[0] % 1 == 0) {
      printf("BCnt[0]=%i BCnt[1]=%i BCnt[2]=%i\n"
	     "  BCnt[3]=%i BCnt[4]=%i BCnt[5]=%i\n",
	     BCnt[0],BCnt[1],BCnt[2],BCnt[3],BCnt[4],BCnt[5]);
      //if(BCnt[0] % 500000 == 0)
      //ShowSigTbl(NumSigs,SigTbl,SigCnt);
      fflush(stdout);
    }

    // Read a permutation B
    x=0;
    assert(scanf("B= %i ",&y) == 1);
    B[x]=y;
    for(x=1;x<KSIZ;x++) {
      assert(scanf(" %i ",&y) == 1);
      B[x]=y;
    }      
    // Check B(0)
    yVal0Ix=xyTbl[0][B[0]];
    assert(yVal0Ix != -1);
    yVal0=yVals[yVal0Ix];
    assert(yVal0 == B[0]);

    // Has yVal0 changed?
    if(yVal0Ix != yVal0IxOld && yVal0IxOld != -1) {
      assert(yVal0Ix > yVal0IxOld);
      yVal0Finish(yVal0IxOld,yVal0Old,BCnt,TotBCnt,
		  SigCnt,TotSigCnt,NumSigs,SigTbl);
    }

    BCnt[0]++;
    yVal0IxOld=yVal0Ix;
    yVal0Old=yVal0;

    memset(BSig,0,sizeof(BSig));
    xySpecial=xySpecialTbl[0][yVal0];
    if(xySpecial != -1) BSig[xySpecial]++;

    // Consider possible conjugations, BTest, by right-translations
    // which satisfy BTest[0]=B[0]

    // Is this B the minimum among all such conjugates?
    MinFg=TRUE;
    // How many such conjugates are there?
    memcpy(BConj[0],B,sizeof(B));
    ConjCnt=1;
    
    for(x0=1;x0<KSIZ;x0++) {
      // Does conjugation by right-translation by x0 give us
      // a new permutation, BTest, where BTest[0] lies among the
      // special values?
      xySpecial=xySpecialTbl[x0][B[x0]];
      if(xySpecial != -1) BSig[xySpecial]++;

      // Does conjugation by right-tranlation by x0 give us a
      // new permutation, BTest, where BTest[0] lies in yVals?
      yValIx=xyTbl[x0][B[x0]];
      if(yValIx != -1) {
	// We assume that during the construction of this permutation B
	// (1) the value B[0]=yVal0=yVals[yVal0Ix] was imposed,
	// and (2) the possibility of conjugating so that
	// B[0]=yVals[yValIx] for yValIx<yVal0Ix was excluded.
	assert(yValIx>=yVal0Ix);
	// Does conjugation by right-translation by x0 give us
	// a new permutation, BTest, with BTest[0]=B[0]?
	if(yValIx == yVal0Ix) {
	  // Calculate BTest
	  x0Inv=KInv2[x0];
	  for(x=0;x<KSIZ;x++) {
	    register int y,xNew,yNew;
	    y=B[x];
	    xNew=KMult2[x][x0Inv];
	    yNew=KMult2[y][x0Inv];
	    BTest[xNew]=yNew;
	  }
	  assert(BTest[0] == yVal0);
	  // If BTest is less than B, we will skip this B
	  Comparison=memcmp(BTest,B,sizeof(B));
	  if(Comparison < 0) {MinFg=FALSE; break;}
	  // If BTest is the same as B, then we don't have a new
	  // conjugate.
	  if(Comparison == 0) continue;
	  // Otherwise, compare BTest to other conjugates we have found
	  FoundFg=FALSE;
	  for(ConjIx=1;ConjIx<ConjCnt;ConjIx++)
	    if(memcmp(BTest,BConj[ConjIx],sizeof(BTest)) == 0)
	      {FoundFg=TRUE; break;}
	  if(FoundFg) continue;
	  // BTest is a new conjugate
	  assert(ConjCnt<MaxConjCnt);
	  memcpy(BConj[ConjCnt],BTest,sizeof(BTest));
	  ConjCnt++;
	}
      }
    }
    // If this B isn't the minimum among conjugations, BTest,  by
    // right-translation satisfying BTest[0]=Special0, don't deal with it
    // at all.
    if(!MinFg)continue;

    // Try to find the signature of B in the signature table
    FoundFg=FALSE;
    for(SigIx=0;SigIx<NumSigs;SigIx++)
      if(memcmp(BSig,SigTbl[SigIx],sizeof(BSig))==0) {FoundFg=TRUE; break;}
    // Add it to the table, if necessary
    if(!FoundFg) {
      assert(NumSigs < MaxNumSigs);
      SigIx=NumSigs;
      memcpy(SigTbl[NumSigs],BSig,sizeof(BSig));
      NumSigs++;
    }
    // Count permutations of each signature
    SigCnt[SigIx]++;

    BCnt[1]+=ConjCnt;
    BCnt[2]++;

    assert(nCycle(KSIZ,3,Perm1,Perm1,Perm1,B));
    assert(nCycle(KSIZ,24,Perm1,Perm1,PermK1,B));
    assert(TrivialPower(KSIZ,3,Perm1,Perm1,PermK2K1u4K2,B));
    assert(TrivialPower(KSIZ,3,Perm1,PermK2K1u4K2,PermK2K1u4K2,B));

    if(nCycle(KSIZ,3,Perm1,Perm1,PermK2K1u4K2,B)) {
      BCnt[3]++;
      if(nCycle(KSIZ,3,Perm1,PermK2K1u4K2,PermK2K1u4K2,B)) {
	BCnt[4]++;
	// IsPUOne((K1^4*K2*A*K1^5*K2*A)^3);
	if(TrivialPower(KSIZ,3,PermK1u4K2,B,PermK1u5K2,B))
	  BCnt[5]++;
      }
    }
  }
  yVal0Finish(yVal0IxOld,yVal0Old,BCnt,TotBCnt,
	      SigCnt,TotSigCnt,NumSigs,SigTbl);

  printf("\nTotBCnt[0]=%i TotBCnt[1]=%i TotBCnt[2]=%i\n"
	 "  TotBCnt[3]=%i TotBCnt[4]=%i TotBCnt[5]=%i\n\n",
	 TotBCnt[0],TotBCnt[1],TotBCnt[2],
	 TotBCnt[3],TotBCnt[4],TotBCnt[5]);
  ShowSigTbl(NumSigs,SigTbl,TotSigCnt);
}

void ShowSigTbl(int NumSigs,INT_t SigTbl[MaxNumSigs][NumSpecial],
		int SigCnt[MaxNumSigs]) {
  int SigIx,SpecialIx;

  printf("NumSigs=%i\n\n",NumSigs);
  printf("SigTbl:\n");

  for(SigIx=0;SigIx<NumSigs;SigIx++) {
    if(SigCnt[SigIx] != 0) {
      printf("%8i  ",SigCnt[SigIx]);
      for(SpecialIx=0;SpecialIx<NumSpecial;SpecialIx++)
	printf("%3i",SigTbl[SigIx][SpecialIx]);
      printf("\n");
    }
  }
  printf("\n");
}

// Finishes off for a particular value of yVal0 (=B[0])
//
// Displays information, updates the total (overall) counts, and rezeroes
// the counts relevant to a particular value of B[0]

void yVal0Finish(int yVal0Ix,int yVal0,int BCnt[],int TotBCnt[],
		 int *SigCnt,int *TotSigCnt,
		 int NumSigs,INT_t SigTbl[MaxNumSigs][NumSpecial]) {

  int SigIx,BCntIx;

  printf("\nyVal0Ix=%i yVal0=%i\n",yVal0Ix,yVal0);
  printf("BCnt[0]=%i BCnt[1]=%i BCnt[2]=%i\n"
	 "  BCnt[3]=%i BCnt[4]=%i BCnt[5]=%i\n\n",
	 BCnt[0],BCnt[1],BCnt[2],BCnt[3],BCnt[4],BCnt[5]);

  if(BCnt[0]==0) {
    for(BCntIx=1;BCntIx<NumBCnt;BCntIx++)
      assert(BCnt[BCntIx]==0);
  }
  else {
    ShowSigTbl(NumSigs,SigTbl,SigCnt);
    for(BCntIx=0;BCntIx<NumBCnt;BCntIx++)
      (TotBCnt[BCntIx])+=(BCnt[BCntIx]);
    memset(BCnt,0,NumBCnt*sizeof(int));

    for(SigIx=0;SigIx<NumSigs;SigIx++) {
      TotSigCnt[SigIx]+=SigCnt[SigIx];
      SigCnt[SigIx]=0;
    }
  }
  fflush(stdout);
}

// Checks whether the product, Q,  of four given permutations,
// satisfies Q^n = 1
// 
inline int TrivialPower(int KSiz,int n,
			INT_t P1[],INT_t P2[],INT_t P3[],INT_t P4[]) {

  int x,y,Ix;
  for(x=0;x<KSiz;x++) {
    y=x;
    for(Ix=0;Ix<n;Ix++)
      y=P1[P2[P3[P4[y]]]];
    if(y != x) return FALSE;
  }
  return TRUE;
}

// Checks whether the product, Q,  of four given permutations,
// satisfies Q^n = 1, and is made up entirely of n-cycles (no fixed
// points allowed).
// 
inline int nCycle(int KSiz,int n,
		  INT_t P1[],INT_t P2[],INT_t P3[],INT_t P4[]) {
  int x,y,Ix;

  for(x=0;x<KSiz;x++) {
    y=x;
    for(Ix=0;Ix<n-1;Ix++) {
      y=P1[P2[P3[P4[y]]]];
      if(y == x) return FALSE;
    }
    y=P1[P2[P3[P4[y]]]];
    if(y != x) return FALSE;
  }
  return TRUE;
}