// Calling sequence
//   
//  PREFIX=cc/cc-dd YIX0=jj YIX1=kk ./permsearch >> 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 searches for permutations~$B$ of cycle
// structure~$3^{|K|/3}$ which act on the left on a group~$K$ and
// satisfy various relations of the form:
//
//    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$.

// A partially completed permutation, together with further pairs
// to be incorporated.
typedef struct {
  int Cnt; // Number of values filled in so far
  INT_t *B,*BInv;  // B itself and its inverse
} B_t;

void AllocB(B_t *B_p,int KSiz);
void InitB(B_t *B_p,int KSiz);
void FreeB(B_t *B_p);
void CopyB(B_t *ToB_p,B_t *FromB_p,int KSiz);
int CompleteB(B_t *B_p,int KSiz,int NumRels,Rels_t *Rels_p,
	      unsigned char *Excluded,
	      Pair_t NewPair);
inline int xyApply(B_t *B_p,int KSiz,int x, int y,
	   Pair_t *ForB3,Pair_t **ForB3Top_p);
inline void xyNew(B_t *B_p,int KSiz,int x,int y,Pair_t *ForRels_p);
int CheckCompletedB(B_t *B_p,int KSiz,int NumRels,Rels_t *Rels_p);
int Check3Completed(INT_t Perm[],int KSiz);
inline void InitPair(Pair_t *Pair_p,B_t *B_p,int KSiz,
		     int NumRels,Rels_t *Rels_p);
void DoFinishedB(B_t *B_p,int KSiz);
void ShowStack(B_t *BStack,Pair_t *PairStack, int StkIx);
void GetYVals(FILE *ifile,int KSiz,int *NumYVals_p,INT_t yVals[]);
void MarkExcluded(unsigned char *Excluded,int yVal,
		  int KSiz,int *KMult,int *KInv,
		  int NumRels,Rels_t *Rels_p);

int FinishedBCnt=0;
FILE *logFile;

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

  // Stack index
  int StkIx;
  // Stack of partial permutations
  B_t BStruct[MaxKSiz];
  // Parallel stack of new facts of the form $B(x)=y$ to incorporate
  // into the partial permutation
  Pair_t NewPair[MaxKSiz];

  // List of y-values to be used in imposing $B(0)=y$ or excluding
  // $B(x)=yx$.
  int NumYVals,yValIx0,yValIx1,yValIx,yVal;
  INT_t yVals[MaxKSiz];

  // Table of excluded pairs $(x,y)$: we don't want $B(x)=y$
  unsigned char *Excluded;

  char *yIx0_p,*yIx1_p;
  int OKFg,Ix;

  // Input files
  FILE *KFile, *KRelFile, *yValFile;

  char *Prefix_p,*NN_p,FileNam[MAX_FILENAME_LEN]="../";

  // Get the prefix
  Prefix_p=getenv("PREFIX");
  assert(Prefix_p != NULL);

  // Get the indexes of the first and last y-values, which will occur
  // the initial pair --- (x=0,yVals[Ix]), which means B(0)=yVals[Ix]
  yIx0_p=getenv("YIX0");
  assert(yIx0_p != NULL);
  assert(sscanf(yIx0_p," %i ",&yValIx0) == 1);
  assert(yValIx0 >=0);
  yIx1_p=getenv("YIX1");
  assert(yIx1_p != NULL);
  assert(sscanf(yIx1_p," %i ",&yValIx1) == 1);
  assert(yValIx1>=yValIx0);

  // Open the files
  {
    size_t len=strlen(Prefix_p),
      leny0=strlen(yIx0_p), leny1=strlen(yIx1_p);
    assert(3+len+1+leny0+1+leny1+16<MAX_FILENAME_LEN);
    // Output files
    strcpy(&FileNam[3],Prefix_p);
    strcpy(&FileNam[3]+len,"-");
    strcpy(&FileNam[3]+len+1,yIx0_p);
    strcpy(&FileNam[3]+len+1+leny0,"-");
    strcpy(&FileNam[3]+len+1+leny0+1,yIx1_p);
    strcpy(&FileNam[3]+len+1+leny0+1+leny1,"-permsearch.log");
    logFile=fopen(FileNam,"w");
    assert(logFile != NULL);
    fprintf(logFile,"%s\n",FileNam);
    // Input files
    strcpy(&FileNam[3]+len,"-K");
    fprintf(logFile,"%s\n",FileNam);
    KFile=fopen(FileNam,"r");
    assert(KFile != NULL);

    strcpy(&FileNam[3]+len,"-KRel");
    fprintf(logFile,"%s\n",FileNam);
    KRelFile=fopen(FileNam,"r");
    assert(KRelFile != NULL);

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

    fflush(logFile);
  }

  // Get the elements of~$K$
  GetKEltsKK(KFile,&KSiz,&KElts,&KMult,&KInv);
  if(strcmp(KElts[0].Word,"") == 0)
    strcpy(KElts[0].Word,"1");
  fprintf(logFile,"GetKEltsKK, KSiz=%i\n",KSiz);
  assert(KSiz < SHRT_MAX);
  // Get the pairs~$(k,k')$ which constrain~$B$
  GetKRels(KRelFile,&NumRels,&Rels,KSiz,KElts,KMult);
  fprintf(logFile,"GetKRels, NumRels=%i\n",NumRels);
  // Get the list of yVals
  GetYVals(yValFile,KSiz,&NumYVals,yVals);
  fprintf(logFile,"GetYVals, NumYVals=%i\n",NumYVals);
  assert(yValIx1<NumYVals);

  Excluded=malloc(KSiz*KSiz*sizeof(unsigned char));
  assert(Excluded != NULL);

  // Mark all the B-structures as not having space allocated
  for(StkIx=0;StkIx<MaxKSiz;StkIx++)
    BStruct[StkIx].B=(INT_t *)NULL;

  for(yValIx=yValIx0;yValIx<=yValIx1;yValIx++) {

    // Set up the table of excluded values
    memset(Excluded,0,KSiz*KSiz*sizeof(unsigned char));
    
    yVal=yVals[yValIx];
    fprintf(logFile,"\nExcluded:\n");
    for(Ix=0;Ix<yValIx;Ix++) {
      fprintf(logFile,"B[x] = k_{%i} x\n",yVals[Ix]);
      MarkExcluded(Excluded,yVals[Ix],KSiz,KMult,KInv,NumRels,&Rels);
    }
    fprintf(logFile,"\nDesired:\nB[0] = k_{%i}\n\n",yVal);

    // Set up the base element on the stack
    StkIx=0;
    InitB(&BStruct[StkIx],KSiz);
    NewPair[StkIx].x=0; NewPair[StkIx].y=yVal;

    while(TRUE) {
      if(StkIx==1) 
	ShowStack(BStruct,NewPair,StkIx);
      if(StkIx < 0)break;
      for(;NewPair[StkIx].y < KSiz;NewPair[StkIx].y++) {
	if((BStruct[StkIx].BInv)[NewPair[StkIx].y] == -1 &&
	   NewPair[StkIx].y != NewPair[StkIx].x) break;
      }
      // Have we considered all possibilities for adding to
      // this partial B?
      if(NewPair[StkIx].y >= KSiz || 
	 (StkIx==0 && NewPair[StkIx].y > yVal)) {
	StkIx--;
	continue;
      }
      StkIx++;
      assert(StkIx < KSiz);
      // Copy the old BStruct to the next step on the stack

      CopyB(&BStruct[StkIx],&BStruct[StkIx-1],KSiz);
      // Add the fact from Pair, and fill in as much as possible based on
      // that and the constraints.
      OKFg=CompleteB(&BStruct[StkIx],KSiz,NumRels,&Rels,Excluded,
		     NewPair[StkIx-1]);
      // Next time, we'll try the next value for~$y$ in $B(x)=y$.
      NewPair[StkIx-1].y++;
      // Was it possible to meet the constraints?
      if(OKFg)  {
	// Have we finished?
	if(BStruct[StkIx].Cnt == KSiz) {
	  DoFinishedB(&BStruct[StkIx],KSiz);
	  StkIx--;
	}
	// If not, set up the Pair which needs to be varied 
	// to further determine the partial B
	else {
	  InitPair(&NewPair[StkIx],&BStruct[StkIx],KSiz,
		   NumRels,&Rels);
	}
      }
      else {
	StkIx--;
      }
    }
    fprintf(logFile,"\nyValIx=%i yVal=%i --- FinishedBCnt=%i\n",
	    yValIx,yVal,FinishedBCnt);
  }
  fprintf(logFile,"\nyValIx0=%s yValIx1=%s --- FinishedBCnt=%i DONE\n",
	  yIx0_p,yIx1_p,FinishedBCnt);
}


// Allocates space for a B-structure
void AllocB(B_t *B_p,int KSiz) {
    (B_p->B)=malloc(KSiz*sizeof(INT_t));
    assert(B_p->B != NULL);
    (B_p->BInv)=malloc(KSiz*sizeof(INT_t));
    assert(B_p->BInv != NULL);
}

// Frees space used by a B-structure
void FreeB(B_t *B_p) {
    free(B_p->B);
    free(B_p->BInv);
}


// Initializes a B-structure
void InitB(B_t *B_p,int KSiz) {
  // If necessary, allocates space
  if(B_p->B == NULL) AllocB(B_p,KSiz);

  (B_p->Cnt)=0;
  memset(B_p->B,-1,KSiz*sizeof(INT_t));
  memset(B_p->BInv,-1,KSiz*sizeof(INT_t));
}

// Copies the data from one B-structure to another
void CopyB(B_t *ToB_p,B_t *FromB_p,int KSiz) {
  // If necessary, allocates space
  if(ToB_p->B == NULL) AllocB(ToB_p,KSiz);

  (ToB_p->Cnt)=(FromB_p->Cnt);
  memcpy(ToB_p->B,FromB_p->B,KSiz*sizeof(INT_t));
  memcpy(ToB_p->BInv,FromB_p->BInv,KSiz*sizeof(INT_t));
}

// Completes a _partial_ permutation B, based on one new datum:
// B(xNew)=yNew.  Returns TRUE if this is possible; FALSE if it is
// not.
int CompleteB(B_t *B_p,int KSiz,int NumRels,Rels_t *Rels_p,
	      unsigned char *Excluded,
	      Pair_t NewPair) {
  int RelIx;
  register int x,y,z,x2,y2;
  INT_t *k1Perm,*k2Perm,*k3Perm,*k4Perm;
  // Pair to which we need to apply the relations of the form
  //
  //   B k1 = k2 B   and B k3 = k4 B^-1
  //
  // This pair will not yet have been applied to B and BInv
  Pair_t ForRels;
  // Pairs to which we need to apply the fact that B consists entirely
  // of 3-cycles.  These pairs will already have been applied to B and
  // BInv.
  Pair_t ForB3[KSiz],*ForB3Top;

  ForB3Top=&ForB3[0];
  memcpy(&ForRels,&NewPair,sizeof(Pair_t));

  while(TRUE) {
    // If there are any pairs $(x,y)$ such that $B(x)=y$ which have
    // not yet been applied to B and BInv, apply them, together with
    // all further pairs deducible from the relations of the form:
    //
    //   B k1 = k2 B  and  B k3 = k4^-1 B
    if(ForRels.x != -1) {
      x=ForRels.x;
      y=ForRels.y;
      // Mark ForRels as empty
      ForRels.x=(-1);
      // It's possible that the fact B(x)=y has already been applied;
      // if not, apply it.
      if((B_p->B)[x] != y) {
	if(Excluded[KSiz*x+y]) return FALSE;

	k1Perm=(Rels_p->k1); k2Perm=(Rels_p->k2);
	k3Perm=(Rels_p->k3); k4Perm=(Rels_p->k4);
	for(RelIx=0;RelIx<NumRels;RelIx++) {
	  if(! xyApply(B_p,KSiz,k1Perm[x],k2Perm[y],ForB3,&ForB3Top))
	    return FALSE;
	  if(! xyApply(B_p,KSiz,k3Perm[y],k4Perm[x],ForB3,&ForB3Top))
	    return FALSE;
	  k1Perm+=KSiz; k2Perm+=KSiz;
	  k3Perm+=KSiz; k4Perm+=KSiz;
	}
      }
    }
    // If pairs are available which have not yet been checked
    // against the requirement that B consist only of 3-cycles,
    // see what can be deduced from them:
    else if(ForB3Top > ForB3) {
      ForB3Top--;
      x=ForB3Top->x;
      y=ForB3Top->y;
      assert((B_p->B)[x]==y);
      assert((B_p->BInv)[y]==x);
      // Proceed forwards ---  B: x->y->z->x
      z=(B_p->B)[y];
      if(z != -1) {
	if(z == x)return FALSE;
	x2=(B_p->B)[z];
	if(x2 != -1) {
	  if(x2 != x)return FALSE;
	}
	else
	  xyNew(B_p,KSiz,z,x,&ForRels);
      }
      else {
	// Proceed backwards --- BInv: y->x->z->y
	z=(B_p->BInv)[x];
	if(z != -1) {
	  if(z==y)return FALSE;
	  y2=(B_p->BInv)[z];
	  if(y2 != -1) {
	    if(y2 != y)return FALSE;
	  }
	  else
	    xyNew(B_p,KSiz,y,z,&ForRels);
	}
      }
    }
    // Nothing is left to do.
    else {
      assert(CheckCompletedB(B_p,KSiz,NumRels,Rels_p));
      return TRUE;
    }
  }
}

// Applies the fact B(x)=y to B and BInv. Returns FALSE if this is
// impossible OR if it has already been done.  Otherwise returns TRUE.
inline int xyApply(B_t *B_p,int KSiz,int x, int y,
		   Pair_t *ForB3,Pair_t **ForB3Top_p) {
  // See if the application is possible AND new
  if((B_p->B)[x] != -1)return FALSE;
  if((B_p->BInv)[y] != -1)return FALSE;
  // Apply to B and BInv
  (B_p->B)[x]=y;
  (B_p->BInv)[y]=x;
  (B_p->Cnt)++;

  // Put this pair on the stack to check against the condition that
  // B consists of only 3-cycles
  assert((*ForB3Top_p) - ForB3 < KSiz);
  (*ForB3Top_p)->x = x;   (*ForB3Top_p)->y = y;
  (*ForB3Top_p)++;
  return TRUE;
}

// Inserts a new fact --- namely B(x)=y --- into *ForRels_p
inline void xyNew(B_t *B_p,int KSiz,int x,int y,Pair_t *ForRels_p) {
  // Check that *ForRels_p is empty
  assert(ForRels_p->x == -1);
  // and add the new fact
  (ForRels_p->x) = x; (ForRels_p->y) = y;
}

// Checks that a _partially_ completed B-structure is consistent
int CheckCompletedB(B_t *B_p,int KSiz,int NumRels,Rels_t *Rels_p) {
  INT_t *k1Perm,*k2Perm,*k3Perm,*k4Perm;
  int RelIx,CntB,CntBInv;
  register int x,y;

  // Check that B and BInv are inverse _partial_ permutations:
  CntB=CntBInv=0;
  for(x=0;x<KSiz;x++) {
    y=(B_p->B)[x];
    if(y>=0 && y<KSiz) {
      CntB++;
      assert((B_p->BInv)[y]==x);
    }
    else
      assert(y==(-1));

    y=(B_p->BInv)[x];
    if(y>=0 && y<KSiz)
      CntBInv++;
    else
      assert(y==(-1));
  }
  assert(CntB == B_p->Cnt);
  assert(CntBInv == B_p->Cnt);

  // Check that all the relations of the form:
  //
  //   B k1 = k2 B  and  B k3 = k4^-1 B
  // 
  // have been applied
  k1Perm=(Rels_p->k1); k2Perm=(Rels_p->k2);
  k3Perm=(Rels_p->k3); k4Perm=(Rels_p->k4);
  for(RelIx=0;RelIx<NumRels;RelIx++) {
    for(x=0;x<KSiz;x++) {
      y=(B_p->B)[x];
      if(y != -1) {
	assert((B_p->B)[k1Perm[x]] == k2Perm[y]);
	assert((B_p->B)[k3Perm[y]] == k4Perm[x]);
	// Here, we are checking that the subgroup generated
	// by B and the k1, k2, k3  and k4 acts freely:
	assert(k1Perm[x] != x || RelIx==0);
	assert(k3Perm[y] != x);
      }
    }
    k1Perm+=KSiz; k2Perm+=KSiz;
    k3Perm+=KSiz; k4Perm+=KSiz;
  }
  // Check that all deductions possible from the fact that
  // B is composed of 3-cycles have already been used:
  assert(Check3Completed(B_p->B,KSiz));

  return TRUE;
}

// Checks whether a _partial_ permutation Perm might be finished so as
// to be made up exclusively of 3-cycles

// Moreover, checks whether all entries deducible from the fact that
// upon completion Perm is to consist exclusively of 3-cycles have
// already been filled in.
int Check3Completed(INT_t Perm[],int KSiz) {
  INT_t PermX[KSiz+1],Iterates[4][KSiz+1];
  int IterIx;
  register int Ix;

  for(Ix=0;Ix<KSiz;Ix++)
    if(Perm[Ix] == -1)
      PermX[Ix]=KSiz;
    else if (Perm[Ix]>=0 && Perm[Ix]<KSiz)
      PermX[Ix]=Perm[Ix];
    else
      assert(FALSE);
  PermX[KSiz]=KSiz;

  for(Ix=0;Ix<KSiz;Ix++)
    Iterates[0][Ix]=Ix;

  for(IterIx=1;IterIx<4;IterIx++)
    for(Ix=0;Ix<=KSiz;Ix++)
      Iterates[IterIx][Ix]=PermX[Iterates[IterIx-1][Ix]];

  // For a chain C0 -> C1 -> C2 -> C3
  for(Ix=0;Ix<KSiz;Ix++) {
    assert(Iterates[0][Ix] != Iterates[1][Ix]); // C0 != C1
    assert(Iterates[0][Ix] != Iterates[2][Ix]); // C0 != C2
    if(Iterates[2][Ix] != KSiz)
      // Here C1 and C2 are filled in
      assert(Iterates[0][Ix] == Iterates[3][Ix]); // C0 == C3
  }
  return TRUE;
}

// Given a _partial_ B, set up the initial value for a pair $(x,y)$
// corresponding to a relation $B(x)=y$ which should be applied to B
// to make it more complete.  Naturally, x should be a value for which
// B is not yet defined.
//
// Moreover, we give preference to an x for whose value will determine
// many others.  If $B(x)=y$ then also
//
//   B(k1 x) = k2 B(x) = k2 y    and    B(k3 y) = k4 B^-1(y) = k4 x
//
// Consequently, if the value of B^-1(k1 x) is known, any value of
// y will determine a 3-cycle:
//
//   B:  B^-1(k1 x) -> k1 x -> k2 y -> B^-1(k1 x)
// 
// and similarly, if the value of B(k4 x) is known, any value of y will
// determine a 3-cycle:
//
//   B:  k3 y -> k4 x -> B(k4 x) -> k3 y
//
// We want as many of these deductions to be available as possible.
inline void InitPair(Pair_t *Pair_p,B_t *B_p,int KSiz,
		     int NumRels,Rels_t *Rels_p) {
  int x,xSav,DedCnt,MaxDedCnt,RelIx;
  INT_t *k1Perm,*k4Perm;

  assert(B_p->Cnt != KSiz);
  
  Pair_p->y=0;
  MaxDedCnt=(-1);

  for(x=0;x<KSiz;x++)
    if((B_p->B)[x] == -1) {
      DedCnt=0;
      k1Perm=(Rels_p->k1); k4Perm=(Rels_p->k4);
      for(RelIx=0;RelIx<NumRels;RelIx++) {
	if((B_p->BInv)[k1Perm[x]] != -1)DedCnt++;
	if((B_p->B)[k4Perm[x]] != -1)DedCnt++;
	k1Perm+=KSiz; k4Perm+=KSiz;
      }
      if(DedCnt>MaxDedCnt) {
	Pair_p->x = x;
	MaxDedCnt=DedCnt;
      }
    }
  //fprintf(logFile,"InitPair X: Pair=(%i,%i) DedCnt=%i\n",
  //	 Pair_p->x,Pair_p->y,MaxDedCnt);
}

// Treat a finished B
void DoFinishedB(B_t *B_p,int KSiz) {
  int Ix,Cnt;
  FinishedBCnt++; 
  // fprintf(logFile,"***FinishedBCnt=%i\n",FinishedBCnt);

  printf("B=\n");
  Cnt=0;
  for(Ix=0;Ix<KSiz;Ix++) {
    printf("%4i",B_p->B[Ix]);
    Cnt++;
    if(Cnt % 20 == 0)printf("\n");
  }
  if(Cnt % 20 != 0)printf("\n");
  printf("\n");
  //fflush(stdout);
}

// Shows some information on the stack of partially completed B's
void ShowStack(B_t *BStack,Pair_t *PairStack, int StkIx) {
  B_t *B_p;
  Pair_t *Pair_p;
  int Ix;

  B_p=BStack;
  Pair_p=PairStack;

  for(Ix=0;Ix<=StkIx;Ix++) {
    fprintf(logFile,"Cnt=%i (%i,%i) ",B_p->Cnt,Pair_p->x,Pair_p->y-1);
    B_p++;
    Pair_p++;
  }
  fprintf(logFile,"\n");
  fflush(logFile);
}


// Marks certain pairs $(x,y)$ as excluded --- we don't want B(x)=y
// when y=k_{yVal} x, i.e. we don't want B(x)=k_{yVal} x
void MarkExcluded(unsigned char *Excluded,int yVal,
		  int KSiz,int *KMult,int *KInv,
		  int NumRels,Rels_t *Rels_p) {
  int KIx,RelIx,x,yVal12,yVal34,y;
  INT_t *k1Perm,*k2Perm,*k3Perm,*k4Perm;

  // If we're excluding:
  //
  //   B(x) = k_{yVal} x   i.e.  B^-1(y) = k_{yVal}^-1 x
  //
  // and if we have relations of the form:
  //
  //   B k1 = k2 B    and    B k3 = k4 B^-1 
  //
  // then we must also exclude:
  //
  //   k_{yVal} k1 x = B k1 x = k2 B x       and
  //   k_{yVal}^-1 k4 x = B^-1 k4 x = k3 B x 
  //
  // i.e.
  //
  //   B x = k2^-1 k_{yVal} k1 x
  //   B x = k3^-1 k_{yVal}^-1 k4 x

  //fprintf(logFile,"MarkExcluded, yVal=%i\n",yVal);

  k1Perm=Rels_p->k1; k2Perm=Rels_p->k2;
  k3Perm=Rels_p->k3; k4Perm=Rels_p->k4;  
  for(RelIx=0;RelIx<NumRels;RelIx++) {
    yVal12=KMult[KSiz*KMult[KSiz*KInv[k2Perm[0]]
			    +yVal]
		 +k1Perm[0]];
    yVal34=KMult[KSiz*KMult[KSiz*KInv[k3Perm[0]]
			    +KInv[yVal]]
		 +k4Perm[0]];
    //fprintf(logFile,"%4i%4i",yVal12,yVal34);
    for(x=0;x<KSiz;x++) {
      y=KMult[KSiz*yVal12+x];
      Excluded[KSiz*x+y]=TRUE;
      y=KMult[KSiz*yVal34+x];
      Excluded[KSiz*x+y]=TRUE;
    }
    k1Perm+=KSiz; k2Perm+=KSiz;
    k3Perm+=KSiz; k4Perm+=KSiz;
  }
  //fprintf(logFile,"\n");
}