// Calling sequence: ./subgp2 FILELIST COND3FILE
//
// FILELIST is itself a file.  Every line contains a pair of filenames
// as follows:
//
//   GOODPERMS  TRANSLATIONS
//
// where GOODPERMS is a file containing good permutations (those
// satifying Conditions 1 and 2)
//
// and TRANSLATIONS is a file containing a list of quadruples that
// correspond to ``translation'' permutations by which we should
// conjugate each of the good permutations.
//
// COND3FILE is for the output: permutations satisfying Condition 3 as
// well.

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

#define TRUE 1
#define FALSE 0

#define MAX(a,b) ( (a)>=(b) ? (a) : (b) )
#define MIN(a,b) ( (a)<=(b) ? (a) : (b) )

#define NN 24

#define TransMax 216

FILE *FileList; // File with list of filenames

FILE *PermFile; // File with the good permutations

FILE *TransFile; // File with ``translations'' by which we should
		 // conjugate each good permutation

FILE *Cond3File; // Output file for permutations satisfying
		 // ``Condition~3'' as well as Conditions~1 and~2

int PermCycleCnt[NN][NN]; // How many generated permutations contain 
                          // the cycle (0 C1 C2)?

int Cond3CycleCnt[NN][NN]; // How many generated permutations
                           // containing the cycle (0 C1 C2) satisfy
                           // Condition~3?


int Cond3FailCnt; // How many times does Condition~3 fail ?
int Cond3IxFailCnt[NN]; // How many times does Condition~3 fail, for
  		        // the first time, at each of the possible
  		        // values 0, 1, 2, ..., 23?

int PermCnt,PermFileLineCnt,Cond3Cnt;
	
int ReadTrans(int *NTrans,int Trans[][4]);
void DisplayTrans(int *NTrans,int Trans[][4]);
int ReadPerm(int Perm[]);
void DisplayCntTbl(int Table[][NN],char *TblNam);
void DisplayPerm(int Perm[]);
int Check3Perm(int Perm[]);
int Check4Perm(int Perm[]);
inline int CheckCond3(int Perm[]);
inline int CheckCond3Strong(int Perm[]);
void DisplayCond3(int Perm[]);
void DisplayNotCond3(int Perm[]);

int main(int argc,char *(argv[])) {
  int Perm0[NN],Perm[NN],Perm4[NN],Trans[TransMax][4],NTrans,TransIx,Ix;
  char PermFileName[121],TransFileName[121];

  PermCnt=Cond3Cnt=Cond3FailCnt=0;
  memset(PermCycleCnt,0,sizeof(PermCycleCnt));
  memset(Cond3CycleCnt,0,sizeof(Cond3CycleCnt));
  memset(Cond3IxFailCnt,0,sizeof(Cond3IxFailCnt));

  assert(argc==3);
  FileList=fopen(argv[1],"r");
  assert(FileList != NULL);
  Cond3File=fopen(argv[2],"w");
  assert(Cond3File != NULL);
  
  while(fscanf(FileList,"%*[ \n\t]") != EOF) {
    assert(fscanf(FileList,"%120[^ \t\n] %120[^ \t\n] ",
		  PermFileName,TransFileName) == 2);
    PermFile=fopen(PermFileName,"r");
    assert(PermFile != NULL);
    TransFile=fopen(TransFileName,"r");
    assert(TransFile != NULL);
    
    printf("*** PermFile=%s TransFile=%s ***\n\n",
	   PermFileName,TransFileName);

    assert(ReadTrans(&NTrans,Trans));
    DisplayTrans(&NTrans,Trans);
    assert(NTrans>0);

    PermFileLineCnt=0;
    while(fscanf(PermFile,"%*[ \n\t]") != EOF) {
      assert(ReadPerm(Perm0));
      PermFileLineCnt++;
      //printf("main A: Perm0=");
      //for(Ix=0;Ix<NN;Ix++)printf("%3i",Ix);
      //printf("\n              ");
      //for(Ix=0;Ix<NN;Ix++)printf("%3i",Perm0[Ix]);
      //printf("\n");


      for(TransIx=0;TransIx<NTrans;TransIx++) {
	// Conjugate Perm0 by the map 
	//   j /mapsto j+T(j)
	// where T(j) depends only on j mod 4 and is given by
	//   T(j) = Trans[TransIx][j % 4]
	for(Ix=0;Ix<NN;Ix++)
	  Perm[(Ix+Trans[TransIx][Ix % 4]) % NN]
	    =(Perm0[Ix]+Trans[TransIx][Perm0[Ix] % 4]) % NN;

	// Check that Perm has cycle structure 8(3)
        assert(Check3Perm(Perm));

	// Set up the modification of Perm given by
	//   j /mapsto Perm[j]+4
	for(Ix=0;Ix<NN;Ix++)Perm4[Ix]=(Perm[Ix]+4) % NN;

	// Check that Perm4 has cycle structure 6(4)
	assert(Check4Perm(Perm4));

	PermCnt++;
	PermCycleCnt[Perm[0]][Perm[Perm[0]]]++;

	if(CheckCond3Strong(Perm)) {
	  Cond3Cnt++;
	  Cond3CycleCnt[Perm[0]][Perm[Perm[0]]]++;
	  for(Ix=0;Ix<NN;Ix++)
	    fprintf(Cond3File,"%3i",Perm[Ix]);
	  fprintf(Cond3File,"\n");
	  if(Cond3Cnt % 100 == 0)DisplayCond3(Perm);
	}

	if(PermCnt % 2160000 == 1) {
	  printf("PermCnt =%11i, Cond3Cnt =%9i, Frac = %09.7f\n",
		 PermCnt,Cond3Cnt,(double) Cond3Cnt/(double) PermCnt);

	  printf("Cond3 failures: ");
	  for(Ix=0;Ix<NN/3;Ix++)
	    printf("%09.7f ",(double) Cond3IxFailCnt[Ix]/(double) Cond3FailCnt);
	  printf("\n                ");
	  for(Ix=NN/3;Ix<2*NN/3;Ix++)
	    printf("%09.7f ",(double) Cond3IxFailCnt[Ix]/(double) Cond3FailCnt);
	  printf("\n                ");
	  for(Ix=2*NN/3;Ix<NN;Ix++)
	    printf("%09.7f ",(double) Cond3IxFailCnt[Ix]/(double) Cond3FailCnt);
	  printf("\n");

	  DisplayPerm(Perm);

	  //if(CheckCond3(Perm))
	  //   DisplayCond3(Perm);
	  //else
	  //  DisplayNotCond3(Perm);

	  fflush(stdout);
	  fflush(Cond3File);
	}
	//if(PermCnt >= 21600000)break;
      }
      //if(PermCnt >= 21600000)break;
      //if(PermFileLineCnt == 1000)break;
    }
    fclose(PermFile);
    printf("\n");
    //if(PermCnt >= 21600000)break;
  }
  DisplayCntTbl(PermCycleCnt,"PermCycleCnt");
  DisplayCntTbl(Cond3CycleCnt,"Cond3CycleCnt");
  printf("PermCnt =%11i, Cond3Cnt =%9i, Frac = %09.7f\n",
	 PermCnt,Cond3Cnt,(double) Cond3Cnt/(double) PermCnt);
}

// Reads the translation table
int ReadTrans(int *NTrans,int Trans[][4]) {
  int Ix;

  *NTrans=0;
  while(fscanf(TransFile,"%*[ \n\t]") != EOF) {
    assert(*NTrans < TransMax);
    for(Ix=0;Ix<4;Ix++)
      assert(fscanf(TransFile,"%i ",&Trans[*NTrans][Ix]) == 1);
    (*NTrans)++;
  }
  fclose(TransFile);
  return TRUE;
}

void DisplayTrans(int *NTrans,int Trans[][4]) {
  int TransIx,Ix;

  for(TransIx=0;TransIx<*NTrans;TransIx++) {
    printf("Trans[%03i] = ",TransIx);
    for(Ix=0;Ix<4;Ix++)printf("%3i",Trans[TransIx][Ix]);
    printf("\n");
  }
  printf("\n");
} 


// Reads one permutation
int ReadPerm(int Perm[]) {
  int Ix;
  for(Ix=0;Ix<NN;Ix++)
    assert(fscanf(PermFile,"%i ",&Perm[Ix]) == 1);

  //printf("ReadPerm A: Perm=");
  //for(Ix=0;Ix<NN;Ix++)printf("%3i",Ix);
  //printf("\n                 ");
  //for(Ix=0;Ix<NN;Ix++)printf("%3i",Perm[Ix]);
  //printf("\n");

  return TRUE;
}

void DisplayCntTbl(int Table[][NN],char *TblNam) {
  int Ix1,Ix2;
  
  for(Ix1=0;Ix1<NN;Ix1++) {
    printf("%s[%i,--]:",TblNam,Ix1);
    for(Ix2=0;Ix2<NN;Ix2++) {
      if(Ix2 % 12 == 0)
	printf("\n[%2i,%2i..%2i]",Ix1,Ix2,MIN(Ix2+11,NN-1));
      printf("%8i",Table[Ix1][Ix2]);
    }
    printf("\n");
  }
  printf("\n");
  return;
}
                                                                                                         
// Displays the cycle representation of a permutation 
void DisplayPerm(int Perm[]) {
  int StartIx,Ix,Displayed[NN];

  memset(Displayed,0,sizeof(Displayed));

  StartIx=0;
  while(StartIx<NN) {
    printf("(");
    Ix=StartIx;
    while(TRUE) {
      printf("%3i",Ix);
      Displayed[Ix]=TRUE;
      Ix=Perm[Ix];
      if(Ix == StartIx)break;
    }
    printf(") ");
    for(StartIx++;StartIx<NN && Displayed[StartIx];StartIx++);
  }
  printf("\n");

  return;
}
      
// Check whether a permutations is all 3-cycles (1-cycles NOT allowed)
int Check3Perm(int Perm[]) {
  int Ix,C0,C1,C2,Checked[NN];

  memset(Checked,FALSE,sizeof(Checked));

  //printf("Check3Perm A: Perm=");
  //for(Ix=0;Ix<NN;Ix++)printf("%3i",Ix);
  //printf("\n                   ");
  //for(Ix=0;Ix<NN;Ix++)printf("%3i",Perm[Ix]);
  //printf("\n");

  for(Ix=0;Ix<NN;Ix++)assert(Perm[Ix]>=0 && Perm[Ix]<NN);

  for(Ix=0;Ix<NN;Ix++) {
    if(Checked[Ix])continue;
    C0=Ix;
    C1=Perm[C0];
    if(C0 == C1)return FALSE;
    C2=Perm[C1];
    if(C0 != Perm[C2])return FALSE;
    Checked[C0]=Checked[C1]=Checked[C2]=TRUE;
  }

  return TRUE;
}

// Check whether a permutations is all 4-cycles (1-cycles NOT allowed)
int Check4Perm(int Perm[]) {
  int Ix,C0,C1,C2,C3,Checked[NN];

  memset(Checked,FALSE,sizeof(Checked));

  for(Ix=0;Ix<NN;Ix++) {
    if(Checked[Ix])continue;
    C0=Ix;
    C1=Perm[C0];
    if(C0 == C1)return FALSE;
    C2=Perm[C1];
    if(C0 == C2)return FALSE;
    C3=Perm[C2];
    if(C0 != Perm[C3])return FALSE;
    Checked[C0]=Checked[C1]=Checked[C2]=Checked[C3]=TRUE;
  }

  return TRUE;
}

// Checks ``Condition 3'' for a permutation.  This condition is that
// if
//
//   \gamma_A(j)=11+Perm[j]
//   \gamma_B(j)=17+Perm[j]
//
// then for every $j$ there exist indices $p$, $q$, and $r$ so that
//
//   \gamma_p(\gamma_p(\gamma_r(j)))=j
//
inline int CheckCond3(int Perm[]) {
  int j,jA,jB,jAA,jAB,jBA,jBB,jAAA,jAAB,jABA,jABB,jBAA,jBAB,jBBA,jBBB;

#define gammaA(x) (Perm[x]<13 ? Perm[x]+11 : Perm[x]-13)
#define gammaB(x) (Perm[x]<7 ? Perm[x]+17 : Perm[x]-7)

  for(j=0;j<24;j++) {
    jA=gammaA(j);
    jB=gammaB(j);

    jAA=gammaA(jA);
    jAB=gammaB(jA);

    jBA=gammaA(jB);
    jBB=gammaB(jB);

    jAAA=gammaA(jAA);
    jAAB=gammaB(jAA);
    jABA=gammaA(jAB);
    jABB=gammaB(jAB);

    jBAA=gammaA(jBA);
    jBAB=gammaB(jBA);
    jBBA=gammaA(jBB);
    jBBB=gammaB(jBB);

    if(j != jAAA && j != jAAB && j != jABA && j != jABB &&
       j != jBAA && j != jBAB && j != jBBA && j != jBBB) {
      Cond3FailCnt++;
      Cond3IxFailCnt[j]++;
      return FALSE;
	}
  }
  return TRUE;
}


// Checks the STRONG version of ``Condition 3'' for a permutation.
// This condition is that if
//
//   \gamma_A(j)=11+Perm[j]
//   \gamma_B(j)=17+Perm[j]
//
// then for every $j$ there exist indices $p$, $q$, $p'$, and $q'$ so that
//
//         \gamma_p (\gamma_q (\gamma_0))) = j
//   \gamma_{p'} (\gamma_{q'} (\gamma_1))) = j
//
inline int CheckCond3Strong(int Perm[]) {
  int j,jA,jB,jAA,jAB,jBA,jBB,jAAA,jAAB,jABA,jABB,jBAA,jBAB,jBBA,jBBB;

  for(j=0;j<24;j++) {
    jA=gammaA(j);
    jB=gammaB(j);

    jAA=gammaA(jA);
    jAB=gammaB(jA);

    jBA=gammaA(jB);
    jBB=gammaB(jB);

    jAAA=gammaA(jAA);
    jAAB=gammaB(jAA);
    jABA=gammaA(jAB);
    jABB=gammaB(jAB);

    jBAA=gammaA(jBA);
    jBAB=gammaB(jBA);
    jBBA=gammaA(jBB);
    jBBB=gammaB(jBB);

    if((j != jAAA && j != jAAB && j != jABA && j != jABB) ||
       (j != jBAA && j != jBAB && j != jBBA && j != jBBB)) {
      Cond3FailCnt++;
      Cond3IxFailCnt[j]++;
      return FALSE;
	}
  }
  return TRUE;
}


// Displays the TRUTH of Condition 3 for a permutation

void DisplayCond3(int Perm[]) {
  int Choice[3],Summand[3],Chain[4],StepIx,Ix,GoodCnt;
  int Options[2]={11,-7};
  char Flag[2];

  printf("\nThe TRUTH of Condition 3 for Perm:\n");
  DisplayPerm(Perm);
  printf("\n");
  printf("   Ix ->  ");
  for(Ix=0;Ix<NN;Ix++)printf("  %02i",Ix);
  printf("\n");
  printf(" Perm[Ix] ");
  for(Ix=0;Ix<NN;Ix++)printf("  %02i",Perm[Ix]);
  printf("\n\n");

  for(Ix=0;Ix<NN;Ix++) {
    GoodCnt=0;
    for(Choice[0]=0;Choice[0]<2;Choice[0]++)
      for(Choice[1]=0;Choice[1]<2;Choice[1]++)
	for(Choice[2]=0;Choice[2]<2;Choice[2]++) {
	  Chain[0]=Ix;
	  for(StepIx=0;StepIx<3;StepIx++) {
	    Summand[StepIx]=Options[Choice[StepIx]];
	    Chain[StepIx+1]=(Perm[Chain[StepIx]]+Summand[StepIx]+24) % 24;
	  }
	  if(Chain[3] != Chain[0])continue;
	  GoodCnt++;
	  if(Chain[0] == MIN(Chain[0],MIN(Chain[1],Chain[2])))
	    strncpy(Flag,"*",2);
	  else
	    strncpy(Flag," ",2);
	  for(StepIx=0;StepIx<3;StepIx++)
	    printf("%02i%1s-> %02i + %2i=",
		   Chain[StepIx],Flag,Perm[Chain[StepIx]],Summand[StepIx]);
	  printf("%02i%1s",Chain[3],Flag);
	  if(GoodCnt>1)printf("  +++ %i +++",GoodCnt);
	  printf("\n");
	}
    assert(GoodCnt>0);
  }
  printf("\n");
  return;
}

// Displays the FAILURE of Condition 3 for a permutation

void DisplayNotCond3(int Perm[]) {
  int Choice[3],Summand[3],Chain[4],StepIx,Ix,OK,ShowFg,
    OldSummand[3],OldChain[4];
  int Options[2]={11,-7};

  printf("\nThe FAILURE of Condition 3 for Perm:\n");
  DisplayPerm(Perm);
  printf("\n");
  printf("   Ix ->  ");
  for(Ix=0;Ix<NN;Ix++)printf("  %02i",Ix);
  printf("\n");
  printf(" Perm[Ix] ");
  for(Ix=0;Ix<NN;Ix++)printf("  %02i",Perm[Ix]);
  printf("\n\n");

  for(Ix=0;Ix<NN;Ix++) {
    OK=FALSE;
    for(Choice[0]=0;Choice[0]<2;Choice[0]++) {
      for(Choice[1]=0;Choice[1]<2;Choice[1]++) {
	for(Choice[2]=0;Choice[2]<2;Choice[2]++) {
	  Chain[0]=Ix;
	  for(StepIx=0;StepIx<3;StepIx++) {
	    Summand[StepIx]=Options[Choice[StepIx]];
	    Chain[StepIx+1]=(Perm[Chain[StepIx]]+Summand[StepIx]+24) % 24;
	  }
	  if(Chain[3] == Chain[0]) {
	    OK=TRUE;
	    break;
	  }
	}
	if(OK)break;
      }
      if(OK)break;
    }
    if(!OK)break;
  }

  assert(!OK);

  // OldSummand and OldChain start out with impossible values
  memset(OldSummand,NN,sizeof(OldSummand));
  memset(OldChain,-1,sizeof(OldChain));

  for(Choice[0]=0;Choice[0]<2;Choice[0]++) 
    for(Choice[1]=0;Choice[1]<2;Choice[1]++) 
	for(Choice[2]=0;Choice[2]<2;Choice[2]++) {
	  Chain[0]=Ix;
	  for(StepIx=0;StepIx<3;StepIx++) {
	    Summand[StepIx]=Options[Choice[StepIx]];
	    Chain[StepIx+1]=(Perm[Chain[StepIx]]+Summand[StepIx]+24) % 24;
	  }
	  ShowFg=FALSE;
	  for(StepIx=0;StepIx<3;StepIx++) {
	    if(OldChain[StepIx] != Chain[StepIx])ShowFg=TRUE;
	    if(ShowFg)
	      printf("%02i -> %02i",
		   Chain[StepIx],Perm[Chain[StepIx]]);
	    else
	      printf("        ");
	    if(OldSummand[StepIx] != Summand[StepIx])ShowFg=TRUE;
	    if(ShowFg)
	      printf(" + %2i =",Summand[StepIx]);
	    else
	      printf("       ");
	  }
	  printf("%02i",Chain[3]);
	  printf("\n");
	  memcpy(OldSummand,Summand,sizeof(OldSummand));
	  memcpy(OldChain,Chain,sizeof(OldChain));
	}
  printf("\n");
  return;
}