// Calling sequence
//   
//  PREFIX=cc/cc-dd MM=mm 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 order~3 with no
// fixed points which act on the left on a set consisting of MM
// copies of the 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 each copy of~$K$ by
// left-multiplication.
//
// The set on which $B$~acts, namely MM~copies of the group~$K$ enters
// into the conditions only through its structure as a left-$K$-set.
// Consequently any permutation which preserves its structure as a
// left $K$-set will conjugate any valid~$B$ to another such. The
// group of such permutations is the wreath-product of the symmetric
// group on MM letters, permuting the MM copies of~$K$, and
// $K$~itself, acting by right multiplication on any single copy
// of~$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 GetKRelsMM(FILE *ifile,int *NumRels_p,Rels_t *Rels_p,
		int MM,int KSiz,Elt_t *KElts,int *KMult);

int NamToPermMM(INT_t *Perm,char *kNam,
		int MM,int KSiz, Elt_t *KElts, int *KMult);

void AllocB(B_t *B_p,int PermSiz);
void InitB(B_t *B_p,int PermSiz);
void FreeB(B_t *B_p);
void CopyB(B_t *ToB_p,B_t *FromB_p,int PermSiz);
int CompleteB(B_t *B_p,int PermSiz,int NumRels,Rels_t *Rels_p,
	      unsigned char *Excluded,
	      Pair_t NewPair);
inline int xyApply(B_t *B_p,int PermSiz,int x, int y,
	   Pair_t *ForB3,Pair_t **ForB3Top_p);
inline void xyNew(B_t *B_p,int x,int y,Pair_t *ForRels_p);
int CheckCompletedB(B_t *B_p,int PermSiz,int NumRels,Rels_t *Rels_p);
int Check3Completed(INT_t Perm[],int PermSiz);
inline void InitPair(Pair_t *Pair_p,B_t *B_p,int PermSiz,
		     int NumRels,Rels_t *Rels_p);
void DoFinishedB(B_t *B_p,int PermSiz);
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 MM,int KSiz,int *KMult,int *KInv,
		  int NumRels,Rels_t *Rels_p);

int FinishedBCnt=0;
FILE *logFile;

int main(int argc,char *(argv[])) {
  // Number of copies of~$K$ to be used
  int MM;
  char *MM_p;
  // Information about the group $K$
  int KSiz,PermSiz,*KMult,*KInv;
  Elt_t *KElts;
  // Information about the relations
  int NumRels;
  Rels_t Rels;

#define MaxStkSiz 15
  // Stack index
  int StkIx;
  // Stack of partial permutations
  B_t BStruct[MaxStkSiz];
  // Parallel stack of new facts of the form $B(x)=y$ to incorporate
  // into the partial permutation
  Pair_t NewPair[MaxStkSiz];
  // Parallel stack: last value to use for~$y$
  INT_t yLast[MaxStkSiz];

  // 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 value of~$m$
  MM_p=getenv("MM");
  assert(MM_p != NULL);
  assert(sscanf(MM_p," %i ",&MM) == 1);
  assert(MM>=1);
  assert(MM<=6);  // No reason for this upper limit except that it
                  // holds for all the examples we have in mind.

  // 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);
  PermSiz=MM*KSiz;
  assert(PermSiz < SHRT_MAX);
  // Get the pairs~$(k_1,k_2)$ and $(k_3,k_4)$ which constrain~$B$
  GetKRelsMM(KRelFile,&NumRels,&Rels,MM,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(PermSiz*PermSiz*sizeof(unsigned char));
  assert(Excluded != NULL);

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

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

    // Set up the table of excluded values
    memset(Excluded,0,PermSiz*PermSiz*sizeof(unsigned char));
    
    yVal=yVals[yValIx];
    // yVal=KSiz should come last
    assert(yVal<KSiz || (yVal==KSiz && yValIx==yValIx1));
    // ... and is not to be used unless MM>1
    if(yVal==KSiz && MM==1)break;

    fprintf(logFile,"\nExcluded:\n");
    for(Ix=0;Ix<yValIx;Ix++) {
      fprintf(logFile,"B[x] = k_{%i} x\n",yVals[Ix]);
      MarkExcluded(Excluded,yVals[Ix],MM,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],PermSiz);
    NewPair[StkIx].x=0; NewPair[StkIx].y=yVal;
    yLast[StkIx]=yVal;

    while(TRUE) {
      if(StkIx==1) 
	ShowStack(BStruct,NewPair,StkIx);
      if(StkIx < 0)break;
      for(;NewPair[StkIx].y <= yLast[StkIx];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 > yLast[StkIx]) {
	StkIx--;
	continue;
      }
      StkIx++;
      assert(StkIx < MaxStkSiz);

      // Copy the old BStruct to the next step on the stack
      CopyB(&BStruct[StkIx],&BStruct[StkIx-1],PermSiz);
      // Add the fact from Pair, and fill in as much as possible based on
      // that and the constraints.
      OKFg=CompleteB(&BStruct[StkIx],PermSiz,NumRels,&Rels,Excluded,
		     NewPair[StkIx-1]);
      // Next time, we'll try the next value for~$y$ in $B(x)=y$.
      INT_t yOld=NewPair[StkIx-1].y;
      NewPair[StkIx-1].y++;
      // Was it possible to meet the constraints?
      if(OKFg)  {
	// Have we finished?
	if(BStruct[StkIx].Cnt == PermSiz) {
	  DoFinishedB(&BStruct[StkIx],PermSiz);
	  StkIx--;
	}
	// If not, set up the Pair which needs to be varied 
	// to further determine the partial B
	else {
	  InitPair(&NewPair[StkIx],&BStruct[StkIx],PermSiz,
		   NumRels,&Rels);
	  // fprintf(logFile,"MMsearch Q1: NewPair.x=%i NewPair.y=%i\n",
	  //	 NewPair[StkIx].x,NewPair[StkIx].y);

	  // and the limit value for $y$
	  yLast[StkIx]=yLast[StkIx-1];
	  if(yLast[StkIx] != PermSiz-1) {
	    if(yLast[StkIx] < KSiz) 
	      yLast[StkIx]=KSiz;
	    else
	      if(yOld >= yLast[StkIx] || NewPair[StkIx].x >= yLast[StkIx])
		yLast[StkIx]+=KSiz;
	    if(yLast[StkIx] >= PermSiz) yLast[StkIx]=PermSiz-1;
	  }
	  // fprintf(logFile,"MMsearch Q2: yLast=%i\n",yLast[StkIx]);
	}
      }
      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 PermSiz) {
    (B_p->B)=malloc(PermSiz*sizeof(INT_t));
    assert(B_p->B != NULL);
    (B_p->BInv)=malloc(PermSiz*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 PermSiz) {
  // If necessary, allocates space
  if(B_p->B == NULL) AllocB(B_p,PermSiz);

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

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

  (ToB_p->Cnt)=(FromB_p->Cnt);
  memcpy(ToB_p->B,FromB_p->B,PermSiz*sizeof(INT_t));
  memcpy(ToB_p->BInv,FromB_p->BInv,PermSiz*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 PermSiz,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[PermSiz],*ForB3Top;

  // fprintf(logFile,"CompleteB A: NewPair.x=%i NewPair.y=%i\n",
  // NewPair.x,NewPair.y);
  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) {
	// fprintf(logFile,"CompleteB E: PermSiz=%i x=%i y=%i\n",PermSiz,x,y);
	if(Excluded[PermSiz*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,PermSiz,k1Perm[x],k2Perm[y],ForB3,&ForB3Top))
	    return FALSE;
	  if(! xyApply(B_p,PermSiz,k3Perm[y],k4Perm[x],ForB3,&ForB3Top))
	    return FALSE;
	  k1Perm+=PermSiz; k2Perm+=PermSiz;
	  k3Perm+=PermSiz; k4Perm+=PermSiz;
	}
      }
    }
    // 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,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,y,z,&ForRels);
	}
      }
    }
    // Nothing is left to do.
    else {
      assert(CheckCompletedB(B_p,PermSiz,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 PermSiz,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 < PermSiz);
  (*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 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 PermSiz,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<PermSiz;x++) {
    y=(B_p->B)[x];
    if(y>=0 && y<PermSiz) {
      CntB++;
      assert((B_p->BInv)[y]==x);
    }
    else
      assert(y==(-1));

    y=(B_p->BInv)[x];
    if(y>=0 && y<PermSiz)
      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<PermSiz;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+=PermSiz; k2Perm+=PermSiz;
    k3Perm+=PermSiz; k4Perm+=PermSiz;
  }
  // Check that all deductions possible from the fact that
  // B is composed of 3-cycles have already been used:
  assert(Check3Completed(B_p->B,PermSiz));

  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 PermSiz) {
  INT_t PermX[PermSiz+1],Iterates[4][PermSiz+1];
  int IterIx;
  register int Ix;


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

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

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

  // For a chain C0 -> C1 -> C2 -> C3
  for(Ix=0;Ix<PermSiz;Ix++) {
    assert(Iterates[0][Ix] != Iterates[1][Ix]); // C0 != C1
    assert(Iterates[0][Ix] != Iterates[2][Ix]); // C0 != C2
    if(Iterates[2][Ix] != PermSiz)
      // 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 PermSiz,
		       int NumRels,Rels_t *Rels_p) {
    int x,xSav,DedCnt,MaxDedCnt,RelIx;
    INT_t *k1Perm,*k4Perm;

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

    for(x=0;x<PermSiz;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+=PermSiz; k4Perm+=PermSiz;
	}
	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 PermSiz) {
    int Ix,Cnt;
    FinishedBCnt++; 
    // fprintf(logFile,"***FinishedBCnt=%i\n",FinishedBCnt);

    // printf("B=\n");
    printf("[\n");
    Cnt=0;
    for(Ix=0;Ix<PermSiz;Ix++) {
      printf("%3i",B_p->B[Ix]);
      if(Ix<PermSiz-1)printf(",");
      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 MM,int KSiz,int *KMult,int *KInv,
		    int NumRels,Rels_t *Rels_p) {
    int KIx,RelIx,x,yVal12,yVal34,PermSiz=MM*KSiz;
    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++) {
	int y12,y34,MMShift;
	y12=KMult[KSiz*yVal12+x];
	y34=KMult[KSiz*yVal34+x];
	for(MMShift=0;MMShift<PermSiz;MMShift+=KSiz) {
	  Excluded[PermSiz*(x+MMShift)+(y12+MMShift)]=TRUE;
	  Excluded[PermSiz*(x+MMShift)+(y34+MMShift)]=TRUE;
	}
      }
      k1Perm+=PermSiz; k2Perm+=PermSiz;
      k3Perm+=PermSiz; k4Perm+=PermSiz;
    }
    //fprintf(logFile,"\n");
  }

  // Sets up the pairs $(k_1,k_2)$ and $(k_3,k_4)$ for the various
  // relations of the form
  //
  //    B k_1 = k_2 B    and    B k_3 = k_4 B^-1
  // 
  void GetKRelsMM(FILE *ifile,int *NumRels_p,Rels_t *Rels_p,
		  int MM,int KSiz,Elt_t *KElts,int *KMult) {

    int PermSiz=MM*KSiz;
    int RelIx;
    INT_t *k1Perm,*k2Perm,*k3Perm,*k4Perm;
    char k1Nam[MaxKEltNamLength],k2Nam[MaxKEltNamLength];
    char k3Nam[MaxKEltNamLength],k4Nam[MaxKEltNamLength];

    assert(fscanf(ifile,"NumRels=%i ",NumRels_p) == 1);
    (Rels_p->k1)=malloc((*NumRels_p)*PermSiz*sizeof(INT_t));
    assert((Rels_p->k1) != NULL);
    (Rels_p->k2)=malloc((*NumRels_p)*PermSiz*sizeof(INT_t));
    assert((Rels_p->k2) != NULL);
    (Rels_p->k3)=malloc((*NumRels_p)*PermSiz*sizeof(INT_t));
    assert((Rels_p->k3) != NULL);
    (Rels_p->k4)=malloc((*NumRels_p)*PermSiz*sizeof(INT_t));
    assert((Rels_p->k4) != NULL);

    k1Perm=(Rels_p->k1); k2Perm=(Rels_p->k2);
    for(RelIx=0;RelIx<*NumRels_p;RelIx++) {
      assert(fscanf(ifile,"+ " ReadKEltNamCnv ReadKEltNamCnv,
		    k1Nam,k2Nam) == 2);
      assert(NamToPermMM(k1Perm,k1Nam,MM,KSiz,KElts,KMult));
      k1Perm+=PermSiz;
      assert(NamToPermMM(k2Perm,k2Nam,MM,KSiz,KElts,KMult));
      k2Perm+=PermSiz;
    }

    k3Perm=(Rels_p->k3); k4Perm=(Rels_p->k4);
    for(RelIx=0;RelIx<*NumRels_p;RelIx++) {
      assert(fscanf(ifile,"- " ReadKEltNamCnv ReadKEltNamCnv,
		    k3Nam,k4Nam) == 2);
      assert(NamToPermMM(k3Perm,k3Nam,MM,KSiz,KElts,KMult));
      k3Perm+=PermSiz;
      assert(NamToPermMM(k4Perm,k4Nam,MM,KSiz,KElts,KMult));
      k4Perm+=PermSiz;
    }

    CheckKRels(*NumRels_p,PermSiz,Rels_p);
  }

  // Converts the name of an element of~$K$ into the permutation
  // corresponding to left-multiplication by that element on MM copies
  // of~$K$
  int NamToPermMM(INT_t *Perm,char *kNam,
		  int MM,int KSiz, Elt_t *KElts, int *KMult) {

    Elt_t *KElt_p;
    int KIx,Ix,*KMult_p,PermSiz=MM*KSiz,MMShift;

    // fprintf(logFile,"NamToPermMM A: kNam=%s KSiz=%i\n",kNam,KSiz);
    KElt_p=KElts;
    KMult_p=KMult;
    for(KIx=0;KIx<KSiz;KIx++) {
      if(strcmp(KElt_p->Word,kNam) == 0) {
	// fprintf(logFile,"NamToPermMM B: KElt_p->Word=%s\n",KElt_p->Word);
	for(Ix=0;Ix<KSiz;Ix++) 
	  for(MMShift=0;MMShift<PermSiz;MMShift+=KSiz)
	    Perm[Ix+MMShift]=KMult_p[Ix]+MMShift;
	return TRUE;
      }
      KElt_p++;
      KMult_p+=KSiz;
    }
    return FALSE;
  }