#ifndef PERMSEARCH_BASE_C
#define PERMSEARCH_BASE_C

#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"

#define INT_t short int

#define MaxKSiz 864

// Structure containing the permutations corresponding to the various
// pairs $(k_1,k_2)$ and $(k_3,k_4)$.  These are for relations of the
// form
// 
//   B k1 = k2 B     and     B k3 = k4 B^-1
typedef struct {
  // Each of these has the structure of an array of _NumRels_
  // permutations, each of which is an array of _PermSiz_ values.
  INT_t *k1,*k2,*k3,*k4;
} Rels_t;

// A pair of values to be used in constructing a permutation.  This
// corresponds to a choice of the form B(x)=y.
typedef struct {
  INT_t x,y; // Pair of values: B(x) = y
} Pair_t;

void GetKRels(FILE *ifile,int *NumRels_p,Rels_t *Rels_p,
	      int KSiz,Elt_t *KElts,int *KMult);
int NamToPerm(INT_t *Perm,char *kNam,
	      int KSiz, Elt_t *KElts, int *KMult);
void CheckKRels(int NumRels, int PermSiz,Rels_t *Rels_p);
void GetYVals(FILE *ifile,int KSiz,int *NumYVals,INT_t yVals[]);

// 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 GetKRels(FILE *ifile,int *NumRels_p,Rels_t *Rels_p,
	      int KSiz,Elt_t *KElts,int *KMult) {

  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)*KSiz*sizeof(INT_t));
  assert((Rels_p->k1) != NULL);
  (Rels_p->k2)=malloc((*NumRels_p)*KSiz*sizeof(INT_t));
  assert((Rels_p->k2) != NULL);
  (Rels_p->k3)=malloc((*NumRels_p)*KSiz*sizeof(INT_t));
  assert((Rels_p->k3) != NULL);
  (Rels_p->k4)=malloc((*NumRels_p)*KSiz*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(NamToPerm(k1Perm,k1Nam,KSiz,KElts,KMult));
    k1Perm+=KSiz;
    assert(NamToPerm(k2Perm,k2Nam,KSiz,KElts,KMult));
    k2Perm+=KSiz;
  }

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

  CheckKRels(*NumRels_p,KSiz,Rels_p);
}

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

  Elt_t *KElt_p;
  int KIx,Ix,*KMult_p;

  KElt_p=KElts;
  KMult_p=KMult;
  for(KIx=0;KIx<KSiz;KIx++) {
    if(strcmp(KElt_p->Word,kNam) == 0) {
      for(Ix=0;Ix<KSiz;Ix++) 
	Perm[Ix]=KMult_p[Ix];
      return TRUE;
    }
    KElt_p++;
    KMult_p+=KSiz;
  }
  return FALSE;
}

// Checks a collection of pairs $(k_1,k_2)$, which give conditions of
// the form:
//
//    B k_1 = k_2 B  and B k_1 = k_2 B^-1
// 
// for consistency and completeness.
void CheckKRels(int NumRels, int PermSiz,Rels_t *Rels_p) {

  INT_t *k1C_p,*k2C_p,*k1D_p,*k2D_p,*k1_p,*k2_p;
  INT_t *k3C_p,*k4C_p,*k3D_p,*k4D_p,*k3_p,*k4_p;
  INT_t Perm1[PermSiz],Perm2[PermSiz];
  int RelIxC,RelIxD,RelIx,Ix;
  int OKFg;

  k1C_p=(Rels_p->k1); k2C_p=(Rels_p->k2);
  k3C_p=(Rels_p->k3); k4C_p=(Rels_p->k4);
  for(RelIxC=0;RelIxC<NumRels;RelIxC++) {
    // Check that the first of the pairs $(k_1,k_2)$ corresponds to
    // the trivial relation,  B 1 = 1 B.
    if(RelIxC==0)
      for(Ix=0;Ix<PermSiz;Ix++) {
	assert(k1C_p[Ix]==Ix);
	assert(k2C_p[Ix]==Ix);
      }
    k1D_p=(Rels_p->k1); k2D_p=(Rels_p->k2);
    k3D_p=(Rels_p->k3); k4D_p=(Rels_p->k4);
    for(RelIxD=0;RelIxD<NumRels;RelIxD++) {
      if(RelIxC != RelIxD) {
	assert(memcmp(k1C_p,k1D_p,PermSiz*sizeof(INT_t)) != 0);
	assert(memcmp(k2C_p,k2D_p,PermSiz*sizeof(INT_t)) != 0);
	assert(memcmp(k3C_p,k3D_p,PermSiz*sizeof(INT_t)) != 0);
	assert(memcmp(k4C_p,k4D_p,PermSiz*sizeof(INT_t)) != 0);
      }

      // Compose a relation of the form: B k1C = k2C B
      // with another of the form:       B k1D = k2D B
      for(Ix=0;Ix<PermSiz;Ix++) {
	Perm1[Ix]=k1D_p[k1C_p[Ix]];
	Perm2[Ix]=k2D_p[k2C_p[Ix]];
      }
      OKFg=FALSE;
      k1_p=(Rels_p->k1); k2_p=(Rels_p->k2);
      for(RelIx=0;RelIx<NumRels;RelIx++) {
	if(memcmp(k1_p,Perm1,PermSiz*sizeof(INT_t)) == 0 &&
	   memcmp(k2_p,Perm2,PermSiz*sizeof(INT_t)) == 0) {
	  OKFg=TRUE;
	  break;
	}
	k1_p+=PermSiz; k2_p+=PermSiz;
      }
      assert(OKFg);

      // Compose a relation of the form: B k1C = k2C B
      // with another of the form:       B k3D = k4D B^-1
      for(Ix=0;Ix<PermSiz;Ix++) {
	Perm1[Ix]=k3D_p[k2C_p[Ix]];
	Perm2[Ix]=k4D_p[k1C_p[Ix]];
      }
      OKFg=FALSE;
      k3_p=(Rels_p->k3);  k4_p=(Rels_p->k4);
      for(RelIx=0;RelIx<NumRels;RelIx++) {
	if(memcmp(k3_p,Perm1,PermSiz*sizeof(INT_t)) == 0 &&
	   memcmp(k4_p,Perm2,PermSiz*sizeof(INT_t)) == 0) {
	  OKFg=TRUE;
	  break;
	}
	k3_p+=PermSiz; k4_p+=PermSiz;
      }
      assert(OKFg);

      // Compose a relation of the form: B k3C = k4C B^-1
      // with another of the form:       B k1D = k2D B
      for(Ix=0;Ix<PermSiz;Ix++) {
	Perm1[Ix]=k1D_p[k3C_p[Ix]];
	Perm2[Ix]=k2D_p[k4C_p[Ix]];
      }
      OKFg=FALSE;
      k3_p=(Rels_p->k3);  k4_p=(Rels_p->k4);
      for(RelIx=0;RelIx<NumRels;RelIx++) {
	if(memcmp(k3_p,Perm1,PermSiz*sizeof(INT_t)) == 0 &&
	   memcmp(k4_p,Perm2,PermSiz*sizeof(INT_t)) == 0) {
	  OKFg=TRUE;
	  break;
	}
	k3_p+=PermSiz; k4_p+=PermSiz;
      }
      assert(OKFg);

      // Compose a relation of the form: B k3C = k4C B^-1
      // with another of the form:       B k3D = k4D B^-1
      for(Ix=0;Ix<PermSiz;Ix++) {
	Perm1[Ix]=k3D_p[k4C_p[Ix]];
	Perm2[Ix]=k4D_p[k3C_p[Ix]];
      }
      OKFg=FALSE;
      k1_p=(Rels_p->k1);  k2_p=(Rels_p->k2);
      for(RelIx=0;RelIx<NumRels;RelIx++) {
	if(memcmp(k1_p,Perm1,PermSiz*sizeof(INT_t)) == 0 &&
	   memcmp(k2_p,Perm2,PermSiz*sizeof(INT_t)) == 0) {
	  OKFg=TRUE;
	  break;
	}
	k1_p+=PermSiz; k2_p+=PermSiz;
      }
      assert(OKFg);

      k1D_p+=PermSiz; k2D_p+=PermSiz; k3D_p+=PermSiz; k4D_p+=PermSiz; 
    }
    k1C_p+=PermSiz; k2C_p+=PermSiz; k3C_p+=PermSiz; k4C_p+=PermSiz;
  }  
}

// Gets the list of yVals
void GetYVals(FILE *ifile,int KSiz,int *NumYVals,INT_t yVals[]) {
  int yValIx,y;
  
  (*NumYVals)=0;
  while(fscanf(ifile,"%*[ \t\n]") != EOF) {
    assert(*NumYVals < MaxKSiz);
    assert(fscanf(ifile," %i ",&y) == 1);
    assert(y>=0);
    assert(y<=KSiz);
    yVals[*NumYVals]=y;
    (*NumYVals)++;
  }
}
#endif