#define TRUE 1
#define FALSE 0

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

#define eps (1E-10)
#define eps2 (1E-13)

typedef double complex Mtx_t[3][3];

// Identity matrix
Mtx_t Id3=
  {{ 1, 0, 0},
   { 0, 1, 0},
   { 0, 0, 1}};

/* Root of $z^3+3z^2+3=0$ */
#define ZZ (-3.2790187861665935794914426057761899)

/* $\sqrt(-7)$ */
#define sqm7 (2.645751311064590590501615753639*I)

// We assume our unitary group preserves the form $\diag(1,1,-R_0^2)$;
// consequently it acts on the ball of radius $R0$ in $\CC^2$.

// R0 = 1
#define R0sq (1.0)
#define R0 (1.0)

// The form with respect to which the original, given matrices are
// unitary:
Mtx_t FF0=
  {{ 2, 0, 0},
   { 0, 0, 1},
   { 0, 1, 0}};

// We will conjugate the matrices of our group by:
Mtx_t ConjMtx=
  {{1,   0,    0},
   {0, 0.5,  0.5},
   {0, 0.5, -0.5}};
  
Mtx_t ConjMtxInv=
  {{1, 0,  0},
   {0, 1,  1},
   {0, 1, -1}};

// The form with respect to which matrices of the group are
// unitary after being conjugated
Mtx_t FF=
  {{ 1, 0,     0},
   { 0, 1,     0},
   { 0, 0, -R0sq}};

// Point of the complex ball (of radius R0)
typedef struct {
  double s2;  // Squared Euclidean distance from the origin
  double complex z[2]; // Coordinates 
} Pt_t;

#define MaxWordLength 100
#define ReadWordDirective " %100s "
// Group element
typedef struct {
  Pt_t p;  // Image of the origin
  Mtx_t m;  // Normalized matrix
  int wLength; // Word length
  unsigned char w[MaxWordLength+1]; // Word
} Elt_t;

#define MaxNumPts 3000
Pt_t Pts[MaxNumPts];

FILE *PtsFile;
FILE *EltsFile;
FILE *hsnorms;

void PrintPts(Pt_t Pts[],int N);
void PrintPt(Pt_t Pt);
void OutputPts(FILE *file,Pt_t Pts[],int N);
void OutputPt(FILE *file,Pt_t Pt);
void PrintElts(Elt_t Elts[],int N);
void PrintElt(Elt_t Elt);
void OutputElts(FILE *file,Elt_t Elts[],int N);
void OutputElt(FILE *file,Elt_t Elt);
inline void OutputHS2(FILE *file,Elt_t Elt);
int ReadPts(FILE *file,Pt_t Pts[],int MaxNPts);
int ReadElts(FILE *file,Elt_t Elts[],int MaxNElts);
void PrintMtx(Mtx_t m);
void OutputMtx(FILE *file,Mtx_t m);
int EqPt(Pt_t p1,Pt_t p2);
int EqPtEps2(Pt_t p1,Pt_t p2);
int EqMtx(Mtx_t m1,Mtx_t m2);
void MultMtx(Mtx_t result,Mtx_t m1,Mtx_t m2);
void InvMtx(Mtx_t result,Mtx_t m);
inline void NormalizeMtx(Mtx_t mtx);
void ConjugateMtx(Mtx_t result,Mtx_t m);
void AdjMtx(Mtx_t result,Mtx_t mtx);
int CheckUnitary(Mtx_t mtx, Mtx_t form);
void ApplyMtxToElt(Elt_t *newElt,Mtx_t mtx,Elt_t Elt);
inline void ApplyMtxToPt(Pt_t *newPt,Mtx_t mtx,Pt_t Pt);
void ApplyGenToElt(Elt_t *newElt,int GenIx,Elt_t Elt,
		   Mtx_t Gens[],char letters[]);
void MakeGpK();
void PrintMatrices(Mtx_t Matrices[],int N,char *Name);
int InDirichlet(Pt_t P,int NumPts);
double sVal(Pt_t Theta,int NumPts,int *PtIx_p);
inline double complex Scal(Pt_t Pt1,Pt_t Pt2);
inline double RealScal(Pt_t Pt1,Pt_t Pt2);
inline double PtAbsSq(Pt_t Pt);
inline double QQ(Pt_t Pt1,Pt_t Pt2);
inline double QQ1(Pt_t Pt1,Pt_t Pt2);
inline double sDist(Pt_t Pt1,Pt_t Pt2);
inline void UToZero(Pt_t *VNew,Pt_t V,Pt_t U);
inline void ZeroToU(Pt_t *VNew,Pt_t V,Pt_t U);
inline void Midpoint(Pt_t *result, Pt_t U, Pt_t V);
int Order(Mtx_t Mtx,Pt_t *Pt_p);
int ComparePts(const void *Ptr1,const void *Ptr2);
int CompareElts(const void *Ptr1,const void *Ptr2);
int CompareEltsWithLength(const void *Ptr1,const void *Ptr2);