#ifndef NUMERIC_BASE_C
#define NUMERIC_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"

// The identity matrix
Mtx_t Id3=
    {{ 1, 0, 0},
     { 0, 1, 0},
     { 0, 0, 1}};

// The form with respect to which (conjugated) matrices are unitary:
Mtx_t FF1=
    {{ 1, 0,  0},
     { 0, 1,  0},
     { 0, 0, -1}};

Pt_t ZeroPt={0,0};

inline int EqPt(const Pt_t p1,const Pt_t p2,const FLOAT_t eps);
inline void NormalizePt(Pt_t Pt);
inline void AddPts(Pt_t Sum,const Pt_t Pt1,const Pt_t Pt_2);
inline void SubPts(Pt_t Diff,const Pt_t Pt1,const Pt_t Pt_2);
inline void NegPt(Pt_t Neg,const Pt_t Pt);
inline void HalvePt(Pt_t Pt);
inline FLOAT_t EucDist2(const Pt_t Pt1,const Pt_t Pt_2);
inline COMPLEX_t Scal(const Pt_t Pt1,const Pt_t Pt2);
inline FLOAT_t RealScal(const Pt_t Pt1,const Pt_t Pt2);
inline FLOAT_t PtAbsSq(Pt_t const Pt);
inline void NormalizeDiff(Pt_t d,const Pt_t u);

inline FLOAT_t QQ(const Pt_t Pt1,const Pt_t Pt2);
inline FLOAT_t QQ1(const Pt_t Pt1,const Pt_t Pt2);
inline FLOAT_t sDist(const Pt_t Pt1,const Pt_t Pt2);
     
void OutputPt(FILE *ofile,const Pt_t Pt);

inline int EqMtx(Mtx_t m1,Mtx_t m2,const FLOAT_t eps);
inline void MultMtx(Mtx_t result,Mtx_t m1,Mtx_t m2);
inline void InvMtx(Mtx_t result,Mtx_t m);
inline void NormalizeMtx(Mtx_t Mtx);
inline void ConjugateMtx(Mtx_t result,Mtx_t m,Mtx_t c,Mtx_t cInv);
inline void AdjMtx(Mtx_t result,Mtx_t Mtx);
int CheckUnitary(Mtx_t Mtx, Mtx_t form);
int CheckPU(Mtx_t Mtx,Mtx_t form,const FLOAT_t eps);
inline void ApplyMtxToPt(Pt_t newPt,Mtx_t mtx,Pt_t Pt);
inline void OriginToZMtx(Mtx_t C, Mtx_t CInv,const  Pt_t Z);
void OutputMtx(FILE *file, Mtx_t m);

inline FLOAT_t Det3By3Real(FLOAT_t Mtx[3][3]);
inline COMPLEX_t Det3By3(Mtx_t Mtx);
inline void Solve3By3Real(FLOAT_t X[3],FLOAT_t M[3][3],FLOAT_t C[3]);

inline void FixUpElt(Elt_t *elt_p);
void ApplyGenToElt(Elt_t *newElt_p,Elt_t Gen,Elt_t oldElt);
void ApplyEltToEltKK(Elt_t *newElt_p,Elt_t Gen,Elt_t oldElt,
		     Elt_t *KElts,int *KMult,int KSiz);
int CompareElts(const void *Ptr1,const void *Ptr2);
inline int CompareEltsPtsOnly(const void *Ptr1,const void *Ptr2);
int CompareEltsWithWords(const void *Ptr1,const void *Ptr2);
void OutputElt(FILE *ofile,Elt_t elt);
void InputElt(FILE *ifile,Elt_t *elt_p);

void GetEltsString(FILE *ifile,char *NumString,
		   int *NumElts_p,Elt_t **Elts);
void GetElts(FILE *ifile,int *NumElts_p,Elt_t **Elts);
void GetFewElts(FILE *ifile,int *NumFewElts_p,Elt_t **FewElts);
void GetManyElts(FILE *ifile,int *NumManyElts_p,Elt_t **ManyElts);
void GetKKCosets(FILE *ifile,int *NumManyElts_p,Elt_t **ManyElts);
void GetFewKKCosets(FILE *ifile,int *NumManyElts_p,Elt_t **ManyElts);
void GetKElts(FILE *ifile,int *KSiz_p,Elt_t **Elts);
void GetKEltsKK(FILE *ifile,int *KSiz_p,Elt_t **Elts,
		int **KMult_p,int **KInv_p);
void OutputKMult(FILE *ofile,int *KMult,int *KInv,int KSiz);
void SetupPts(Pt_t **Pts_p,int *NumPts_p,
	      Elt_t const *Elts,int const NumElts);
void SetupPtsKK(Pt_t **Pts_p,int *NumPts_p,
		Elt_t *Elts,int const NumElts,
		Elt_t *KElts,int const KSiz);

// Checks 2 points for equality (within $\epsilon$)
inline int EqPt(const Pt_t p1,const Pt_t p2,const FLOAT_t eps) {
  if(CABS(p1[0]-p2[0]) > eps)return FALSE;
  if(CABS(p1[1]-p2[1]) > eps)return FALSE;
  return TRUE;
}

// Compare two points of the unit ball of~$\CC^2$: the order to use is
// somewhat arbitrary.  Here it is lexocographic, with the first
// number tested being the distance to~$0$
int ComparePts(const void *Ptr1,const void *Ptr2) {
  const Pt_t *Pt1=(const Pt_t *)Ptr1, *Pt2=(const Pt_t *)Ptr2;
  FLOAT_t AbsSq1,AbsSq2;
  
  // Compare the distances to 0
  AbsSq1=PtAbsSq(*Pt1);
  AbsSq2=PtAbsSq(*Pt2);
  if(AbsSq1 > AbsSq2+epsBig)return +1;
  if(AbsSq2 > AbsSq1+epsBig)return -1;

  // Compare the coordinates of Pt1 and Pt2
  if(CREAL((*Pt1)[0])<CREAL((*Pt2)[0])-epsBigger)return +1;
  if(CREAL((*Pt2)[0])<CREAL((*Pt1)[0])-epsBigger)return -1;
  if(CIMAG((*Pt1)[0])<CIMAG((*Pt2)[0])-epsBigger)return +1;
  if(CIMAG((*Pt2)[0])<CIMAG((*Pt1)[0])-epsBigger)return -1;

  if(CREAL((*Pt1)[1])<CREAL((*Pt2)[1])-epsBigger)return +1;
  if(CREAL((*Pt2)[1])<CREAL((*Pt1)[1])-epsBigger)return -1;
  if(CIMAG((*Pt1)[1])<CIMAG((*Pt2)[1])-epsBigger)return +1;
  if(CIMAG((*Pt2)[1])<CIMAG((*Pt1)[1])-epsBigger)return -1;

  return 0;
}

// Normalize a point in $\CC^2$ so that it has Euclidean length one
inline void NormalizePt(Pt_t Pt) {
  FLOAT_t factor,absz1,absz2;
  int ir,ic;

  absz1=CABS(Pt[0]); absz2=CABS(Pt[1]);
  factor=1.0/SQRT(absz1*absz1+absz2*absz2);
  Pt[0]*=factor; Pt[1]*=factor;
}

// Add two points in $\CC^2$
inline void AddPts(Pt_t Sum,const Pt_t Pt1,const Pt_t Pt2) {
  Sum[0]=Pt1[0]+Pt2[0];
  Sum[1]=Pt1[1]+Pt2[1];
}

// Subtract two points in $\CC^2$
inline void SubPts(Pt_t Diff,const Pt_t Pt1,const Pt_t Pt2) {
  Diff[0]=Pt1[0]-Pt2[0];
  Diff[1]=Pt1[1]-Pt2[1];
}

// Calculate the negative of a point in $\CC^2
inline void NegPt(Pt_t Neg,const Pt_t Pt) {
  Neg[0]=(-Pt[0]);
  Neg[1]=(-Pt[1]);
};

// Halve a point in~$\CC^2$
inline void HalvePt(Pt_t Pt) {
  Pt[0]*=0.5; 
  Pt[1]*=0.5; 
};


// Calculate the Euclidean distance of two points in $\CC^2$
inline FLOAT_t EucDist2(const Pt_t Pt1,const Pt_t Pt2) {
  FLOAT_t a1,a2;

  a1=CABS(Pt1[0]-Pt2[0]);
  a2=CABS(Pt1[1]-Pt2[1]);
  return a1*a1+a2*a2;
}


/* Calculate the scalar product of two points in the complex ball */
inline COMPLEX_t Scal(const Pt_t Pt1,const Pt_t Pt2) {
  return Pt1[0]*CONJ(Pt2[0])+Pt1[1]*CONJ(Pt2[1]);
}

/* Calculate the _real_ scalar product of two points in the complex ball */
inline FLOAT_t RealScal(const Pt_t Pt1,const Pt_t Pt2) {
  return CREAL(Scal(Pt1,Pt2));
}

// Calculate the squared norm of a point
inline FLOAT_t PtAbsSq(Pt_t const Pt) {
  FLOAT_t a1,a2;
  a1=CABS(Pt[0]);
  a2=CABS(Pt[1]);
  return a1*a1+a2*a2;
}

// In $\CC^2$, given a unit vector base point~u and a doubled
// displacement vector~d, this routine normalizes~d so that
//
//   (u + new_d)
//
// is a unit vector in the same direction as
//
//   (u + old_d)
inline void NormalizeDiff(Pt_t d,const Pt_t u) {
  FLOAT_t Abs2SumDiff,Abs2Sum,AbsSum,uFact,dFact;

  Pt_t Test,TestV,dOrig;
  memcpy(dOrig,d,sizeof(Pt_t));

  //printf("Normalize Diff:\n");
  //printf("d="); OutputPt(stdout,d);
  //printf("u="); OutputPt(stdout,u);
	 

  assert(FABS(PtAbsSq(u) - 1.0) < epsSmall);
  Abs2SumDiff=2.0*RealScal(u,d)+PtAbsSq(d);
  Abs2Sum=1.0+Abs2SumDiff;
  AbsSum=SQRT(Abs2Sum);
  uFact=(-Abs2SumDiff)/(AbsSum*(1.0+AbsSum));
  dFact=1.0/AbsSum;

  d[0]=dFact*d[0]+uFact*u[0];
  d[1]=dFact*d[1]+uFact*u[1];
  
  AddPts(Test,u,d);
  AddPts(TestV,u,dOrig);
  //printf("Test ="); OutputPt(stdout,Test);
  //printf("TestV="); OutputPt(stdout,TestV);
  assert(FABS(PtAbsSq(Test))-1.0 < epsSmall);
  assert(CIMAG(CONJ(Test[0])*TestV[0]) < epsSmall);
  assert(CIMAG(CONJ(Test[1])*TestV[1]) < epsSmall);
  assert(CABS(Test[0]*TestV[1]-Test[1]*TestV[0]) < epsSmall);
}

// Finds the hyperbolic midpoint between two points of the unit ball
// in~$\CC^2$.
inline void Midpoint(Pt_t result,const  Pt_t U,const  Pt_t V) {
  COMPLEX_t VFact,UFact,den,scal;
  FLOAT_t halfQQ,UVQQ;  

  scal=1.0-Scal(U,V);
  UFact=CABS(scal)*SQRT(1.0-PtAbsSq(V));
  VFact=scal*SQRT(1.0-PtAbsSq(U));
  den=(UFact+VFact);
  UFact/=den; VFact/=den;

  result[0] = UFact*U[0]+VFact*V[0];
  result[1] = UFact*U[1]+VFact*V[1];

  halfQQ=QQ(result,U);
  UVQQ=QQ(U,V);

  assert(FABS(QQ(result,V)-halfQQ) < epsSmall);
  assert(FABS((2.0*halfQQ-1.0)*(2.0*halfQQ-1.0) - UVQQ) < epsSmall);
}


// This calculates, for two points in the complex ball of radius 1,
// $\cosh^2(\hypdist(Pt1,Pt2))$.
inline FLOAT_t QQ(const Pt_t Pt1,const Pt_t Pt2) {
  FLOAT_t a1;

  a1=CABS(1.0-Scal(Pt1,Pt2));
  return (a1*a1)/((1.0-PtAbsSq(Pt1))*(1.0-PtAbsSq(Pt2)));
}

// This calculates $\cosh^2(\hypdist(Pt1,Pt2)$ \emph{times a factor
// depending only on Pt2}.

inline FLOAT_t QQ1(const Pt_t Pt1,const Pt_t Pt2) {
  FLOAT_t a1;

  a1=CABS(1.0-Scal(Pt1,Pt2));
  return (a1*a1)/(1.0-PtAbsSq(Pt1));
}

// This calculates, for two points in the complex ball of radius 1,
// $\tanh(\hypdist(Pt1,Pt2))$.  If one of the points is the origin,
// that value is the Euclidean norm of the other.
inline FLOAT_t sDist(const Pt_t Pt1,const Pt_t Pt2) {
  return SQRT(1.0 - 1.0/QQ(Pt1,Pt2));
}

// Ouputs a point to a file
void OutputPt(FILE *ofile,const Pt_t Pt) {
  int Ix;
  fprintf(ofile,FOShortCnv " ",PtAbsSq(Pt));
  fprintf(ofile,"(" FOSpCnv FOPlCnv "*i," FOSpCnv FOPlCnv "*i)"
	      ,Pt[0],Pt[1]);
  fprintf(ofile,"\n");
}

// Checks 2 matrices for equality
inline int EqMtx(Mtx_t m1,Mtx_t m2,const FLOAT_t eps) {
  int ir,ic;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      if(CABS(m1[ir][ic]-m2[ir][ic]) > eps)return FALSE;
  return TRUE;
}

/* Multiply 2 matrices */
inline void MultMtx(Mtx_t result,Mtx_t m1,Mtx_t m2) {
  int ir,i,ic;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) {
      result[ir][ic]=0;
      for(i=0;i<3;i++)
	result[ir][ic]+=m1[ir][i]*m2[i][ic];
    }
  return;
} 

// Invert a matrix
inline void InvMtx(Mtx_t result,Mtx_t m) {
  COMPLEX_t det,invdet;
  int ir,ic,ix;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      result[ir][ic]=
         m[(ic+1)%3][(ir+1)%3]*m[(ic+2)%3][(ir+2)%3]
	-m[(ic+1)%3][(ir+2)%3]*m[(ic+2)%3][(ir+1)%3];
  det=0;
  for(ix=0;ix<3;ix++)det+=m[1][ix]*result[ix][1];
  invdet=1.0/det;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      result[ir][ic]*=invdet;
  return;
}

/* Normalize a matrix, mutiplying it by a scalar. */
inline void NormalizeMtx(Mtx_t Mtx) {
  COMPLEX_t f;
  int ir,ic;

  f=1.0/Mtx[2][2];
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) 
      Mtx[ir][ic]*=f;
  return;
}

/* Calculate the conjugate matrix cmc^-1 */
inline void ConjugateMtx(Mtx_t result,Mtx_t m,Mtx_t c,Mtx_t cInv) {
  Mtx_t tmp;
  MultMtx(tmp,c,m);
  MultMtx(result,tmp,cInv);
  return;
} 

/* Calculate the (standard) adjoint of a matrix */
inline void AdjMtx(Mtx_t result,Mtx_t Mtx) {
  int ir,ic;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      result[ir][ic]=CONJ(Mtx[ic][ir]);
  return;
} 

/* Checks for unitarity with respect to a sesquilinear form */
inline int CheckUnitary(Mtx_t Mtx,Mtx_t form) {
  int ir,ic,OK;
  Mtx_t adj,prod1,prod2;
  AdjMtx(adj,Mtx);
  MultMtx(prod1,adj,form);
  MultMtx(prod2,prod1,Mtx);
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) 
      if(CABS(prod2[ir][ic]-form[ir][ic])>epsBig)return FALSE;
  return TRUE;
}

// Checks for projective unitarity with respect to a sesquilinear
// form
inline int CheckPU(Mtx_t Mtx,Mtx_t form,FLOAT_t eps) {
  int ir,ic,OK;
  Mtx_t adj,prod1,prod2;
  FLOAT_t fact;
  AdjMtx(adj,Mtx);
  MultMtx(prod1,adj,form);
  MultMtx(prod2,prod1,Mtx);
  assert(CABS(form[2][2]) > epsBig);
  fact=prod2[2][2]/form[2][2];
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) 
      if(CABS(prod2[ir][ic]-fact*form[ir][ic])>eps)return FALSE;
  return TRUE;
}

// Apply a matrix, on the left, as a LFT, to a point of the unit
// ball in~$\CC^2$.
inline void ApplyMtxToPt(Pt_t newPt,Mtx_t mtx,Pt_t Pt) {
  COMPLEX_t f;

  f=1.0/(mtx[2][0]*Pt[0]+mtx[2][1]*Pt[1]+mtx[2][2]);
  newPt[0]=f*(mtx[0][0]*Pt[0]+mtx[0][1]*Pt[1]+mtx[0][2]);
  newPt[1]=f*(mtx[1][0]*Pt[0]+mtx[1][1]*Pt[1]+mtx[1][2]);
}

// Given a point~$z$ in the unit ball of~$\CC^2$, find a matrix,
// projectively unitary with respect to $\diag(1,1,-1)$ which maps~$0$
// to~$z$ as a LFT.  Find also the (projective) inverse of this
// matrix.
inline void OriginToZMtx(Mtx_t C, Mtx_t CInv,const Pt_t Z) {
  FLOAT_t AbsZ2,Root;

  AbsZ2=PtAbsSq(Z);
  Root=SQRT(1.0-AbsZ2);

  C[0][0] = Z[0]*CONJ(Z[0])/(1.0+Root) + Root;
  C[0][1] = Z[0]*CONJ(Z[1])/(1.0+Root);
  C[1][0] = Z[1]*CONJ(Z[0])/(1.0+Root);
  C[1][1] = Z[1]*CONJ(Z[1])/(1.0+Root) + Root;

  C[0][2]=Z[0]; C[2][0]=CONJ(Z[0]);
  C[1][2]=Z[1]; C[2][1]=CONJ(Z[1]);

  C[2][2]=CInv[2][2]=1.0;

  memcpy(CInv,C,sizeof(Mtx_t));
  CInv[0][2]=(-CInv[0][2]); CInv[2][0]=(-CInv[2][0]); 
  CInv[1][2]=(-CInv[1][2]); CInv[2][1]=(-CInv[2][1]); 

  //{
  //  Mtx_t Prod;
  //  Pt_t Pt;
  //  assert(CheckPU(C,FF1,epsSmall));
  //  assert(CheckPU(CInv,FF1,epsSmall));
  //
  //  ApplyMtxToPt(Pt,C,ZeroPt);
  //  assert(EqPt(Z,Pt,epsSmall));
  //  ApplyMtxToPt(Pt,CInv,Z);
  //  assert(EqPt(ZeroPt,Pt,epsSmall));
  //  MultMtx(Prod,C,CInv);
  //  NormalizeMtx(Prod);
  //  assert(EqMtx(Prod,Id3,epsBig));
  //}
}

/* Output a matrix to a file. */
void OutputMtx(FILE *ofile, Mtx_t m) {
  FLOAT_t x,y;
  int ir,ic;
  for(ir=0;ir<3;ir++) {
    for(ic=0;ic<3;ic++) {
      x=CREAL(m[ir][ic]); y=CIMAG(m[ir][ic]);
      if(FABS(x)<epsSmall) x=0.0; if(FABS(y)<epsSmall) y=0.0;
      fprintf(ofile,"  " FOSpCnv FOPlCnv "*i ",x,y);
    }
    fprintf(ofile,"\n");
  }
  return;
}

// Calculates the determinant of a real 3 by 3 matrix
inline FLOAT_t Det3By3Real(FLOAT_t Mtx[3][3]){
  return 
     Mtx[0][0]*(Mtx[1][1]*Mtx[2][2]-Mtx[2][1]*Mtx[1][2])
    +Mtx[1][0]*(Mtx[2][1]*Mtx[0][2]-Mtx[0][1]*Mtx[2][2])
    +Mtx[2][0]*(Mtx[0][1]*Mtx[1][2]-Mtx[1][1]*Mtx[0][2]);
}

// Calculates the determinant of a complex 3 by 3 matrix
inline COMPLEX_t Det3By3(Mtx_t Mtx){
  return 
     Mtx[0][0]*(Mtx[1][1]*Mtx[2][2]-Mtx[2][1]*Mtx[1][2])
    +Mtx[1][0]*(Mtx[2][1]*Mtx[0][2]-Mtx[0][1]*Mtx[2][2])
    +Mtx[2][0]*(Mtx[0][1]*Mtx[1][2]-Mtx[1][1]*Mtx[0][2]);
}

// Solves the real 3 by 3 linear system: MX=C.
inline void Solve3By3Real(FLOAT_t X[3],FLOAT_t M[3][3],FLOAT_t C[3]){
  FLOAT_t det;
  det=Det3By3Real(M);
  X[0]=1.0/det*(
     C[0]*(M[1][1]*M[2][2]-M[2][1]*M[1][2])
    +C[1]*(M[2][1]*M[0][2]-M[0][1]*M[2][2])
    +C[2]*(M[0][1]*M[1][2]-M[1][1]*M[0][2]));
  X[1]=1.0/det*(
     C[0]*(M[1][2]*M[2][0]-M[2][2]*M[1][0])
    +C[1]*(M[2][2]*M[0][0]-M[0][2]*M[2][0])
    +C[2]*(M[0][2]*M[1][0]-M[1][2]*M[0][0]));
  X[2]=1.0/det*(
     C[0]*(M[1][0]*M[2][1]-M[2][0]*M[1][1])
    +C[1]*(M[2][0]*M[0][1]-M[0][0]*M[2][1])
    +C[2]*(M[0][0]*M[1][1]-M[1][0]*M[0][1]));
}


// Fixes up a partially prepared element:
//   (.) normalizes the matrix
//   (.) calculates the image point, $g(0)$
//   (.) calculates the squared distance $d(g(0),0)^2$

inline void FixUpElt(Elt_t *elt_p){
  NormalizeMtx(elt_p->Mtx);
  elt_p->Pt[0]=elt_p->Mtx[0][2];
  elt_p->Pt[1]=elt_p->Mtx[1][2];
  elt_p->d2=PtAbsSq(elt_p->Pt);
}

/* Apply a generator, on the left to a group element. */
void ApplyGenToElt(Elt_t *NewElt_p,Elt_t Gen,Elt_t OldElt) {
  size_t LenGenWord,LenEltWord;

  // Concatenate the two words
  LenGenWord=strlen(Gen.Word);
  LenEltWord=strlen(OldElt.Word);
  assert(LenGenWord+LenEltWord<=MaxWordLength);
  if(strcmp(Gen.Word,"1") == 0) {
    strcpy(NewElt_p->Word,OldElt.Word);
  }
  else {
    strcpy(NewElt_p->Word,Gen.Word);
    if(strcmp(OldElt.Word,"1") != 0)
      strcpy((NewElt_p->Word)+LenGenWord,OldElt.Word);
  }
  // Add up the two lengths
  NewElt_p->Len = Gen.Len+OldElt.Len;
  // Multiply the two matrices
  MultMtx(NewElt_p->Mtx,Gen.Mtx,OldElt.Mtx);
  NormalizeMtx(NewElt_p->Mtx);
  // Set up the point and the squared distance
  NewElt_p->Pt[0]=NewElt_p->Mtx[0][2];
  NewElt_p->Pt[1]=NewElt_p->Mtx[1][2];
  NewElt_p->d2=PtAbsSq(NewElt_p->Pt);
  // Set up the flags
  NewElt_p->KIx=-1;
  NewElt_p->KStart=Gen.KIx;
}

// Apply a group element, on the left to a second group element.  When
// constructing the word, take account of the possibility that the
// first element is from K and that the second element starts with
// another element of K.
void ApplyEltToEltKK(Elt_t *NewElt_p,Elt_t OldElt1,Elt_t OldElt2,
		     Elt_t *KElts,int *KMult,int KSiz){
  size_t LenOldElt1Word,LenOldKWord,LenNewKWord;
  int NewKIx;

  // Multiply the two matrices
  MultMtx(NewElt_p->Mtx,OldElt1.Mtx,OldElt2.Mtx);
  // Normalize the matrix, set up the point, set up the squared
  // distance
  FixUpElt(NewElt_p);

  //printf("ApplyEltToEltKK: %s * %s",OldElt1.Word,OldElt2.Word);

  assert(OldElt1.KStart >=0 && OldElt1.KStart < KSiz);
  assert(OldElt2.KStart >=0 && OldElt2.KStart < KSiz);

  // Is OldElt1 in K?
  if(OldElt1.KIx != -1) {
    assert(OldElt1.KIx == OldElt1.KStart);
    NewKIx=KMult[KSiz*OldElt1.KIx+OldElt2.KStart];
    assert(NewKIx >= 0 && NewKIx < KSiz);
    // Set up length
    NewElt_p->Len=
      OldElt2.Len-KElts[OldElt2.KStart].Len+KElts[NewKIx].Len;
    // Set up word
    LenOldKWord=strlen(KElts[OldElt2.KStart].Word);
    LenNewKWord=strlen(KElts[NewKIx].Word);
    if(strlen(OldElt2.Word)-LenOldKWord+LenNewKWord
       <= MaxWordLength) {
      strcpy(NewElt_p->Word,KElts[NewKIx].Word);
      strcpy((NewElt_p->Word)+LenNewKWord,(OldElt2.Word)+LenOldKWord);
    }
    else {
      strcpy(NewElt_p->Word,KElts[NewKIx].Word);
      strcpy((NewElt_p->Word)+LenNewKWord,"...");
    }
    // Set up KIx and KStart
    NewElt_p->KStart=NewKIx;
    if(OldElt2.KIx == -1)
      NewElt_p->KIx=(-1);
    else
      NewElt_p->KIx=NewKIx;
  }
  else {
    // Set up length
    NewElt_p->Len=OldElt1.Len+OldElt2.Len;
    // Set up word
    LenOldElt1Word=strlen(OldElt1.Word);
    if(LenOldElt1Word+strlen(OldElt2.Word) <= MaxWordLength) {
      strcpy(NewElt_p->Word,OldElt1.Word);
      strcpy((NewElt_p->Word)+LenOldElt1Word,OldElt2.Word);
    }
    else {
      LenNewKWord=strlen(KElts[OldElt1.KStart].Word);
      strcpy(NewElt_p->Word,KElts[OldElt1.KStart].Word);
      strcpy((NewElt_p->Word)+LenNewKWord,"...");
    }
    // Set up KIx and KStart
    NewElt_p->KStart=OldElt1.KStart;
    NewElt_p->KIx=(-1);
  }
  //printf(" = %s ; KIx=%i KStart=%i\n",
  // NewElt_p->Word,NewElt_p->KIx,NewElt_p->KStart);
}


// Compare two group elements: the order to use is somewhat arbitrary.
// Here it is lexocographic, with the first number tested being
// the distance of $g(0)$ from $0$
//
// The word and the ``length'' are NOT taken into account.

int CompareElts(const void *Ptr1,const void *Ptr2) {
  const Elt_t *Elt1=(const Elt_t *)Ptr1, *Elt2=(const Elt_t *)Ptr2;
  int Ix,Row,Col,tmp;

  //printf("\nCompareElts: Elt1=");
  //OutputElt(stdout,*Elt1);
  //printf("Elt2=");
  //OutputElt(stdout,*Elt2);

  // Compare the squared Euclidean distances of $g_1(0)$ to~$0$ and
  // $g_2(0)$ to~$0$
  if(Elt1->d2 > Elt2->d2+epsBig)return +1;
  if(Elt2->d2 > Elt1->d2+epsBig)return -1;

  // Compare the coordinates of $g_1(0)$ and $g_2(0)$
  if(CREAL(Elt1->Pt[0])<CREAL(Elt2->Pt[0])-epsBigger)return +1;
  if(CREAL(Elt2->Pt[0])<CREAL(Elt1->Pt[0])-epsBigger)return -1;
  if(CIMAG(Elt1->Pt[0])<CIMAG(Elt2->Pt[0])-epsBigger)return +1;
  if(CIMAG(Elt2->Pt[0])<CIMAG(Elt1->Pt[0])-epsBigger)return -1;

  if(CREAL(Elt1->Pt[1])<CREAL(Elt2->Pt[1])-epsBigger)return +1;
  if(CREAL(Elt2->Pt[1])<CREAL(Elt1->Pt[1])-epsBigger)return -1;
  if(CIMAG(Elt1->Pt[1])<CIMAG(Elt2->Pt[1])-epsBigger)return +1;
  if(CIMAG(Elt2->Pt[1])<CIMAG(Elt1->Pt[1])-epsBigger)return -1;

  // Compare the coordinates of the matrices of $g_1$ and $g_2$
  for(Row=0;Row<3;Row++)
    for(Col=0;Col<3;Col++) {
    if(CREAL((Elt1->Mtx)[Row][Col])<CREAL((Elt2->Mtx)[Row][Col])-epsBigger)
      return +1;
    if(CREAL((Elt2->Mtx)[Row][Col])<CREAL((Elt1->Mtx)[Row][Col])-epsBigger)
      return -1;
    if(CIMAG((Elt1->Mtx)[Row][Col])<CIMAG((Elt2->Mtx)[Row][Col])-epsBigger)
      return +1;
    if(CIMAG((Elt2->Mtx)[Row][Col])<CIMAG((Elt1->Mtx)[Row][Col])-epsBigger)
      return -1;
    }

  return 0;
}

// Compare two group elements, taking account only of the result of
// the element when applied as a LFT to the origin of $B(\CC^2).
//
// The word and the ``length'' are NOT taken into account.

inline int CompareEltsPtsOnly(const void *Ptr1,const void *Ptr2) {
  const Elt_t *Elt1=(const Elt_t *)Ptr1, *Elt2=(const Elt_t *)Ptr2;
  int Ix,Row,Col,tmp;
 
  return ComparePts(Elt1->Pt,Elt2->Pt);
}

// Compare two group elements: the order to use is somewhat arbitrary.
// Here it is lexocographic, with the first number tested being the
// distance from $g(0)$ to $0$.
//  If the two elements are equal as matrices, which means that they
// represent the same element of the group, one then compares the
// ``lengths'' of their expressions as words and finally the words
// themselves.

int CompareEltsWithWords(const void *Ptr1,const void *Ptr2) {
  const Elt_t *Elt1=(const Elt_t *)Ptr1, *Elt2=(const Elt_t *)Ptr2;
  int tmp;
  
  tmp=CompareElts(Elt1,Elt2);
  if(tmp != 0)return tmp;
  // Compare the word lengths
  if(Elt1->Len > Elt2->Len)return +1;
  if(Elt2->Len > Elt1->Len)return -1;

  return strncmp(Elt1->Word,Elt2->Word,MaxWordLength);
}

// Outputs an element to a file
void OutputElt(FILE *file,Elt_t elt) {
  if(strlen(elt.Word) == 0) {
    assert(elt.KIx == 0);
    fprintf(file,"1 ");
  }
  else
    fprintf(file,"%s ",elt.Word);
  fprintf(file,"Len=%i KStart=%i KIx=%i\n",
	  elt.Len,elt.KStart,elt.KIx);
  fprintf(file,FOShortCnv " (" FOSpCnv FOPlCnv "*i," FOSpCnv FOPlCnv "*i)\n",
	  elt.d2,elt.Pt[0],elt.Pt[1]);
  OutputMtx(file,elt.Mtx);
}

// Inputs an element from a file
void InputElt(FILE *ifile,Elt_t *elt_p) {
  FLOAT_t x,y;
  int rIx,cIx;
  
  assert(fscanf(ifile,ReadWordCnv,&elt_p->Word) == 1);
  assert(fscanf(ifile,"Len=%i ",&elt_p->Len) == 1);
  assert(fscanf(ifile,"KStart=%i ",&elt_p->KStart) == 1);
  assert(fscanf(ifile,"KIx=%i ",&elt_p->KIx) == 1);

  if(strcmp(elt_p->Word,"1") == 0) {
    assert(elt_p->KIx == 0);
    strcpy(elt_p->Word,"");
  }

  assert(fscanf(ifile,FICnv,&elt_p->d2) == 1);
  assert(fscanf(ifile,"(" FICnv FICnv "*i, ",&x,&y) == 2);
  elt_p->Pt[0]=x+I*y;
  assert(fscanf(ifile,FICnv FICnv "*i) ",&x,&y) == 2);
  elt_p->Pt[1]=x+I*y;
  elt_p->d2=PtAbsSq(elt_p->Pt);
  
  for(rIx=0;rIx<3;rIx++)
    for(cIx=0;cIx<3;cIx++) {
      assert(fscanf(ifile,FICnv FICnv "*i ",&x,&y) == 2);
      elt_p->Mtx[rIx][cIx]=x+I*y;
    }
}

// This routine gets elements from a file.  It checks that the first
// line reads: ``NumString''=NumElts, and that the file actually
// contains.  exactly NumElts elements.

void GetEltsString(FILE *ifile,char *NumString,
		   int *NumElts_p,Elt_t **Elts) {
  char Format[55];
  int LenNumString,EltIx;
  LenNumString=strlen(NumString);
  assert(LenNumString<=50);
  strcpy(Format,NumString);
  strcpy(Format+LenNumString,"=%i ");

  assert(fscanf(ifile,Format,NumElts_p) == 1);
  *Elts=malloc((*NumElts_p)*sizeof(Elt_t));

  for(EltIx=0;EltIx<*NumElts_p;EltIx++) {
    InputElt(ifile,&(*Elts)[EltIx]);
  }
  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
}


//void GetElts(FILE *ifile,int *NumElts_p,Elt_t **Elts) {
//  int EltIx;
//
//  assert(fscanf(ifile,"NumElts=%i ",NumElts_p) == 1);
//  *Elts=malloc((*NumElts_p)*sizeof(Elt_t));
//
//  for(EltIx=0;EltIx<*NumElts_p;EltIx++) {
//    InputElt(ifile,&(*Elts)[EltIx]);
//  }
//  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
//}

//void GetFewElts(FILE *ifile,int *NumFewElts_p,Elt_t **FewElts) {
//  int FewEltIx;
//
//  assert(fscanf(ifile,"NumFewElts=%i ",NumFewElts_p) == 1);
//  *FewElts=malloc((*NumFewElts_p)*sizeof(Elt_t));
//
//  for(FewEltIx=0;FewEltIx<*NumFewElts_p;FewEltIx++) {
//    //printf("GetFewElts: FewEltIx=%i\n",FewEltIx);
//    InputElt(ifile,&(*FewElts)[FewEltIx]);
//  }
//  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
//}

//void GetManyElts(FILE *ifile,int *NumManyElts_p,Elt_t **ManyElts) {
//  int ManyEltIx;
//
//  assert(fscanf(ifile,"NumManyElts=%i ",NumManyElts_p) == 1);
//  *ManyElts=malloc((*NumManyElts_p)*sizeof(Elt_t));
//
//  for(ManyEltIx=0;ManyEltIx<*NumManyElts_p;ManyEltIx++) {
//    InputElt(ifile,&(*ManyElts)[ManyEltIx]);
//  }
//  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
//}

//void GetKKCosets(FILE *ifile,int *NumElts_p,Elt_t **Elts) {
//  int EltIx;
//
//  assert(fscanf(ifile,"NumKKCosets=%i ",NumElts_p) == 1);
//  *Elts=malloc((*NumElts_p)*sizeof(Elt_t));
//
//  for(EltIx=0;EltIx<*NumElts_p;EltIx++) {
//    InputElt(ifile,&(*Elts)[EltIx]);
//  }
//  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
//}

//void GetFewKKCosets(FILE *ifile,int *NumFewElts_p,Elt_t **FewElts) {
//  int FewEltIx;
//
//  assert(fscanf(ifile,"NumFewKKCosets=%i ",NumFewElts_p) == 1);
//  *FewElts=malloc((*NumFewElts_p)*sizeof(Elt_t));
//
//  for(FewEltIx=0;FewEltIx<*NumFewElts_p;FewEltIx++) {
//    //printf("GetFewElts: FewEltIx=%i\n",FewEltIx);
//    InputElt(ifile,&(*FewElts)[FewEltIx]);
//  }
//  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
//}

void GetKElts(FILE *ifile,int *KSiz_p,Elt_t **Elts) {
  int EltIx;

  assert(fscanf(ifile,"KSiz=%i ",KSiz_p) == 1);
  *Elts=malloc((*KSiz_p)*sizeof(Elt_t));

  for(EltIx=0;EltIx<*KSiz_p;EltIx++) {
    InputElt(ifile,&(*Elts)[EltIx]);
  }
  assert(fscanf(ifile,"%*[ \n\t]") == EOF);
}


// This routine loads the K-elements and sets up the multiplication
// table and the inverse table for~K
void GetKEltsKK(FILE *ifile,int *KSiz_p,Elt_t **Elts,
		int **KMult_p,int **KInv_p) {
  int EltIx,EltIx2,EltIxProd,EltIxInv,FoundFg;
  Mtx_t Prod,Inv;

  assert(fscanf(ifile,"KSiz=%i ",KSiz_p) == 1);
  *Elts=malloc((*KSiz_p)*sizeof(Elt_t));
  assert(Elts != NULL);
  (*KMult_p)=malloc((*KSiz_p)*(*KSiz_p)*sizeof(int));
  assert((*KMult_p) != NULL);
  (*KInv_p)=malloc((*KSiz_p)*sizeof(int));
  assert((*KInv_p) != NULL);

  for(EltIx=0;EltIx<*KSiz_p;EltIx++) {
    InputElt(ifile,&(*Elts)[EltIx]);
  }
  assert(fscanf(ifile,"%*[ \n\t]") == EOF);

  for(EltIx=0;EltIx<*KSiz_p;EltIx++)
    for(EltIx2=0;EltIx2<*KSiz_p;EltIx2++) {
      MultMtx(Prod,(*Elts)[EltIx].Mtx,(*Elts)[EltIx2].Mtx);
      NormalizeMtx(Prod);
      FoundFg=FALSE;
      for(EltIxProd=0;EltIxProd<*KSiz_p;EltIxProd++)
	if(EqMtx(Prod,(*Elts)[EltIxProd].Mtx,epsSmall)) {
	  (*KMult_p)[EltIx*(*KSiz_p)+EltIx2]=EltIxProd;
	  FoundFg=TRUE;
	  break;
	}
      assert(FoundFg);
    }
  for(EltIx=0;EltIx<*KSiz_p;EltIx++) {
    InvMtx(Inv,(*Elts)[EltIx].Mtx);
    NormalizeMtx(Inv);
    FoundFg=FALSE;
    for(EltIxInv=0;EltIxInv<*KSiz_p;EltIxInv++)
      if(EqMtx(Inv,(*Elts)[EltIxInv].Mtx,epsSmall)) {
	(*KInv_p)[EltIx]=EltIxInv;
	FoundFg=TRUE;
	break;
      }
    assert(FoundFg);
  }
  // OutputKMult(stdout,(*KMult_p),(*KInv_p),*KSiz_p);
}

// Output a K multiplication table
void OutputKMult(FILE *ofile,int *KMult,int *KInv,int KSiz) {
  int Ix1,Ix2;
  fprintf(ofile,"     INV    ");
  for(Ix2=0;Ix2<KSiz;Ix2++)
    fprintf(ofile,"%4i",Ix2);
  fprintf(ofile,"\n\n");
    
  for(Ix1=0;Ix1<KSiz;Ix1++) {
    fprintf(ofile,"%4i%4i    ",Ix1,KInv[Ix1]);
    for(Ix2=0;Ix2<KSiz;Ix2++)
      fprintf(ofile,"%4i",KMult[Ix1*KSiz+Ix2]);
    fprintf(ofile,"\n");
  }
  fprintf(ofile,"\n");
}

// Sets up the list of points in the orbit of the origin on the basis
// of a list of elements g such that g(0) is to be on the list.
// 
// The image points ought already to have been calculated.  This
// routine merely sorts them, eliminates duplicates, and stores them
// in a new place.

void SetupPts(Pt_t **Pts_p,int *NumPts_p,
	      Elt_t const *Elts,int const NumElts) {
  int PtIx;

  // Allocate storage
  *Pts_p=malloc(NumElts*sizeof(Pt_t));
  assert(*Pts_p != NULL);
  // Copy in the points
  for(PtIx=0;PtIx<NumElts;PtIx++)
    memcpy(&(*Pts_p)[PtIx],&(Elts[PtIx].Pt),sizeof(Pt_t));
  // Sort the list
  qsort((*Pts_p),(size_t) NumElts,sizeof(Pt_t),ComparePts);
  // Eliminate Duplicates
  *NumPts_p=1;
  for(PtIx=1;PtIx<NumElts;PtIx++)
    if(! EqPt((*Pts_p)[PtIx],(*Pts_p)[PtIx-1],epsBig)) {
      memcpy(&(*Pts_p)[*NumPts_p],&(*Pts_p)[PtIx],sizeof(Pt_t));
      (*NumPts_p)++;
    }
}

// Sets up the list of points in the orbit of the origin on the basis
// of a list of double-K-cosets KgK such that the orbits Kg(0) ought
// to be included on the list.
// 
// This routine calculates the K-orbits of the images g(0), sorts
// them, eliminates duplicates, and stores them in a new place.

void SetupPtsKK(Pt_t **Pts_p,int *NumPts_p,
		Elt_t *Elts,int const NumElts,
		Elt_t *KElts,int const KSiz) {
  int CosetIx,KIx,PtIx;
  Pt_t *Pt_p;

  // Allocate storage
  *Pts_p=malloc(NumElts*KSiz*sizeof(Pt_t));
  assert(*Pts_p != NULL);
  // Calculate the points
  Pt_p=*Pts_p;
  for(CosetIx=0;CosetIx<NumElts;CosetIx++)
    for(KIx=0;KIx<KSiz;KIx++) {
      ApplyMtxToPt(*Pt_p,KElts[KIx].Mtx,Elts[CosetIx].Pt);
      Pt_p++;
    }
  // Sort the list
  qsort((*Pts_p),(size_t) (NumElts*KSiz),sizeof(Pt_t),ComparePts);
  // Eliminate Duplicates
  *NumPts_p=1;
  for(PtIx=1;PtIx<NumElts*KSiz;PtIx++)
    if(! EqPt((*Pts_p)[PtIx],(*Pts_p)[PtIx-1],epsBig)) {
      memcpy(&(*Pts_p)[*NumPts_p],&(*Pts_p)[PtIx],sizeof(Pt_t));
      (*NumPts_p)++;
    }
}

#endif