#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}};

// 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);