#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>

#define TRUE 1
#define FALSE 0

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

/* Square root of 6 */
#define sq6 (2.44948974278317809819)
#define r sq6
/* 3rd root of unity */
#define z3 (-0.5+0.86602540378443864676*I)
#define om z3

/* $\sqrt(6)-2$ */
#define R0sq (0.44948974278317809819)
/* R0 =$\sqrt(\sqrt(6)-2))$
   We assume our unitary group preserves the form $\diag(1,1,-R0)$;
   consequently it acts on the ball of radius $R0$ in $\CC^2$. */
#define R0 (0.67043996210188582268)

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

typedef double complex Mtx_t[3][3];

// The form with respect to which our matrices are unitary
Mtx_t FF=
  {{ 1, 0,     0},
   { 0, 1,     0},
   { 0, 0, -R0sq}};

// Generators for a big finite subgroup $K$
Mtx_t KU=
  {{ r*(om-1)/6, r*(om-1)/6, 0}, 
   {-r*(om-1)/6, r*(om-1)/6, 0}, 
   {          0,          0, 1}};

Mtx_t KV=
  {{1,  0, 0}, 
   {0, -1, 0}, 
   {0,  0, 1}};

// The group $K$ itself
#define GpKSiz 48
Mtx_t GpK[GpKSiz];

// Other generators of our group $\Gamma$
#define NoGens 10
#define MaxNoGens 10
Mtx_t Gens[MaxNoGens]={
  // This first one is _AA_
  {{      1/3.0*r*om+1/3.0*(2*r+3), 0,            -1/3.0*r*om+1/3.0*r},
   {                             0, 1,                              0}, 
   {1/3.0*(-2*r-6)*om+1/3.0*(-r-3), 0, 1/3.0*(-r-3)*om+1/3.0*(-2*r-3)}},

  {{       -1/6.0*r*om+1/6.0*(r+6), -1/6.0*r*om+1/6.0*r,                            -om}, 
   {       -1/6.0*r*om+1/6.0*(r+6),  1/6.0*r*om-1/6.0*r,                            -om}, 
   {1/3.0*(-2*r-6)*om+1/3.0*(-r-3),                   0, 1/3.0*(-r-3)*om+1/3.0*(-2*r-3)}},

  {{        1/6.0*(r+6)*om-1/6.0*r, -1/6.0*r*om+1/6.0*r,                             om}, 
   {       1/6.0*(-r-6)*om+1/6.0*r, -1/6.0*r*om+1/6.0*r,                            -om}, 
   {1/3.0*(-2*r-6)*om+1/3.0*(-r-3),                   0, 1/3.0*(-2*r-3)*om+1/3.0*(-r-3)}},

  {{ 1/2.0*(r+2), -1/3.0*r*om-1/6.0*r,          -1/3.0*r*om+1/3.0*r}, 
   {1/2.0*(-r-2), -1/3.0*r*om-1/6.0*r,           1/3.0*r*om-1/3.0*r}, 
   {    (r+2)*om,                   0, 1/3.0*(2*r+6)*om+1/3.0*(r+3)}},

  {{         -1/3.0*r*om-1/6.0*r,     1/3.0*(-r-3)*om-1/6.0*r,                             -1}, 
   {  1/3.0*(r+3)*om+1/6.0*(r+6),                1/2.0*(-r-2),         2/3.0*r*om+1/3.0*(r+3)}, 
   {1/2.0*(-r-2)*om+1/2.0*(-r-2), 1/2.0*(-r-2)*om+1/2.0*(r+2), 1/3.0*(-r-3)*om+1/3.0*(-2*r-6)}},

  {{                   1/3.0*r*om+1/6.0*r,  1/3.0*(-r-3)*om-1/6.0*r,                         -om}, 
   {           1/3.0*(r+3)*om+1/6.0*(r+6),             1/2.0*(-r-2),                       -om-2}, 
   {                          1/2.0*(r+2), 1/3.0*(2*r+6)*om-1/6.0*r, 1/3.0*(r+3)*om+1/3.0*(-r-3)}},

  {{     1/6.0*(-r-6)*om+1/6.0*r,     1/2.0*(-r-2)*om+1/2.0*(-r-2),                         om-1}, 
   {1/2.0*(-r-2)*om+1/2.0*(-r-2),         1/6.0*(5*r+6)*om+1/6.0*r, 1/3.0*(r+3)*om+1/3.0*(2*r+3)}, 
   { -1/6.0*r*om+1/6.0*(-5*r-12), 1/6.0*(7*r+18)*om+1/6.0*(5*r+12),                          r+2}},

  {{              1/2.0*(r+2)*om,   1/6.0*(-5*r-6)*om+1/3.0*(-2*r-3), 1/3.0*r*om+1/3.0*(2*r+3)}, 
   {1/6.0*(-r-6)*om+1/3.0*(-r-3),                    1/2.0*(-r-2)*om,                   2*om+1}, 
   {   1/6.0*(5*r+12)*om+1/6.0*r, 1/6.0*(-5*r-12)*om+1/6.0*(-7*r-18),                      r+2}},

  {{                     1/2.0*(-r-2), 1/3.0*(2*r+3)*om+1/6.0*(5*r+12), 1/3.0*(-r+3)*om+1/3.0*(-2*r-3)}, 
   {     1/3.0*(-r-3)*om+1/6.0*(-r-6),              1/3.0*r*om+1/6.0*r,                          -om-1}, 
   {1/6.0*(-5*r-12)*om+1/3.0*(-2*r-6),   1/6.0*(7*r+18)*om+1/3.0*(r+3),                  (-r-2)*om-r-2}},

  {{                 1/2.0*(-r-2), 1/3.0*(r+3)*om+1/6.0*r,  1/3.0*(-r-3)*om+1/3.0*(r-3)}, 
   {1/3.0*(-2*r-3)*om+1/6.0*(r+6),    -1/3.0*r*om-1/6.0*r, 1/3.0*(2*r+3)*om+1/3.0*(r+3)},
   {     1/6.0*r*om+1/3.0*(4*r+9),        1/2.0*(-r-2)*om,                     (r+2)*om}}
  };

#define PtsDesired 5000 // (Approximately) the number of good points we want
#define PtsSiz (2*(GpKSiz+MaxNoGens)*PtsDesired) // Size of full array of points

typedef struct {double s2; double complex z[2];} Pt_t;
Pt_t Pts[PtsSiz]; // Points as generated
Pt_t SavPts[PtsDesired]; // Points from previous list

FILE *PtsFile;

void PrintPts(Pt_t Pts[],int N);
void PrintPt(Pt_t Pt);
void OutputPts(Pt_t Pts[],int N);
void OutputPt(Pt_t Pt);
void PrintMtx(Mtx_t m);
void MultMtx(Mtx_t result,Mtx_t m1,Mtx_t m2);
void AdjMtx(Mtx_t result,Mtx_t mtx);
int CheckUnitary(Mtx_t mtx);
void ApplyMtx(Pt_t *newpt,Mtx_t mtx,Pt_t pt);
void MakeGpK();
void PrintMatrices(Mtx_t Matrices[],int N,char *Name);
double sVal(Pt_t Theta);
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);
int ComparePts(const void *Ptr1,const void *Ptr2);

int main(int argc,char *(argv[])) {
  int Ix,PtIx,NOldPts,NPts2,NPts,Iter,EqualFg,SavedFg;
  Pt_t NewPt,*PtPtr;

  PtsFile=fopen(argv[1],"w");
  assert(PtsFile != NULL);

  printf("FF=\n");  PrintMtx(FF);
  printf("\nKU=\n");  PrintMtx(KU);
  printf("KV=\n");  PrintMtx(KV);
  printf("AA=\n");  PrintMtx(Gens[0]);

  MakeGpK(); // Set up the finite group
  PrintMatrices(GpK,GpKSiz,"GpK");
  for(Ix=0;Ix<GpKSiz;Ix++)assert(CheckUnitary(GpK[Ix]));
  PrintMatrices(Gens,NoGens,"Gens");
  for(Ix=0;Ix<NoGens;Ix++)assert(CheckUnitary(Gens[Ix]));

  exit;

  Pts[0].z[0]=0.0;
  Pts[0].z[1]=0.0;
  Pts[0].s2=PtAbsSq(Pts[0]);
  NOldPts=1;
  SavedFg=FALSE;

  Iter=0;
  while(TRUE) {
    // Apply non-$K$ generators of $\Gamma$ to the old points
    NOldPts=MIN(NOldPts,PtsSiz/(NoGens+1));
    NPts=NOldPts;
    for(Ix=0;Ix<NoGens;Ix++) 
      for(PtIx=0;PtIx<NOldPts;PtIx++) {
	ApplyMtx(&NewPt,Gens[Ix],Pts[PtIx]);
	PtPtr=bsearch(&NewPt,&Pts,(size_t) NOldPts,sizeof(Pt_t),ComparePts);
	if(PtPtr == NULL) {
	  assert(NPts<PtsSiz);
	  Pts[NPts]=NewPt;
	  NPts++;
	}
      }
    // Sort all the points
    qsort(&Pts,(size_t) NPts,sizeof(Pt_t),ComparePts);
    // Eliminate duplicates
    NPts2=1;
    for(PtIx=1;PtIx<NPts;PtIx++)
      if(ComparePts(&Pts[PtIx],&Pts[PtIx-1])>0) {
	Pts[NPts2]=Pts[PtIx];
	NPts2++;
      }
    printf("After Gens: Iter=%i, NOldPts=%i, NPts=%i, NPts2=%i\n",
	   Iter,NOldPts,NPts,NPts2);
    NOldPts=NPts2;

    // Apply elements of $K$ to the old points
    NOldPts=MIN(NOldPts,PtsSiz/(GpKSiz+1));
    NPts=NOldPts;
    for(Ix=0;Ix<GpKSiz;Ix++)
      for(PtIx=0;PtIx<NOldPts;PtIx++) {
	ApplyMtx(&NewPt,GpK[Ix],Pts[PtIx]);
	PtPtr=bsearch(&NewPt,&Pts,(size_t) NOldPts,sizeof(Pt_t),ComparePts);
	if(PtPtr == NULL) {
	  assert(NPts<PtsSiz);
	  Pts[NPts]=NewPt;
	  NPts++;
	}
    }

    // Sort all the points
    qsort(&Pts,(size_t) NPts,sizeof(Pt_t),ComparePts);
    // Eliminate duplicates
    NPts2=1;
    for(PtIx=1;PtIx<NPts;PtIx++)
      if(ComparePts(&Pts[PtIx],&Pts[PtIx-1])>0) {
	Pts[NPts2]=Pts[PtIx];
	NPts2++;
      }
    printf("After  KGp: Iter=%i, NOldPts=%i, NPts=%i, NPts2=%i\n",
	   Iter,NOldPts,NPts,NPts2);
    NOldPts=NPts2;

    // If we have a set of saved points available, and if our new
    // points match those, we're done.
    if(SavedFg && memcmp(&SavPts,&Pts,sizeof(SavPts)) == 0) break;

    // If we have at least PtsDesired points available, make a copy of
    // them.
    if(NOldPts >= PtsDesired) {
      memcpy(&SavPts,&Pts,sizeof(SavPts));
      SavedFg=TRUE;
    }
    Iter++;
  }
  // printf("\n\n===== * =====\n\n");
  PrintPts(SavPts,PtsDesired);
  OutputPts(SavPts,PtsDesired);
}

/* Display a list of points in the complex ball */
void PrintPts(Pt_t Pts[],int N) {
  int Ix;

  Ix=0;
  for(Ix=0;Ix<N;Ix++) {
    PrintPt(Pts[Ix]);
  }
}

/* Display a point in the complex ball */
void PrintPt(Pt_t Pt) {
  printf("% .15f (% .15f + % .15f*i,% .15f + % .15f*i)\n",
	 Pt.s2,
	 creal(Pt.z[0]),cimag(Pt.z[0]),
	 creal(Pt.z[1]),cimag(Pt.z[1]));
  return;
}

/* Output a list of points in the complex ball to PtsFile */
void OutputPts(Pt_t Pts[],int N) {
  int Ix;

  Ix=0;
  for(Ix=0;Ix<N;Ix++) {
    OutputPt(Pts[Ix]);
  }
}

/* Output a point in the complex ball to PtsFile.  Use slightly less
   precision, so that text editing can check whether two sets of
   output are equivalent. */
void OutputPt(Pt_t Pt) {
  double x1,y1,x2,y2;

  x1=creal(Pt.z[0]); y1=cimag(Pt.z[0]);
  x2=creal(Pt.z[1]); y2=cimag(Pt.z[1]);

  if(fabs(x1)<eps) x1=0.0; if(fabs(y1)<eps) y1=0.0;
  if(fabs(x2)<eps) x2=0.0; if(fabs(y2)<eps) y2=0.0;

  fprintf(PtsFile,"% .12f (% .12f + % .12f*i,% .12f + % .12f*i)\n",
	  Pt.s2,x1,y1,x2,y2);
  return;
}

/* Dispay a Matrix */
void PrintMtx(Mtx_t m) {
  int ir,ic;
  for(ir=0;ir<3;ir++) {
    for(ic=0;ic<3;ic++) 
      printf("% .14f%+.14f*i ",m[ir][ic]);
    printf("\n");
  }
  return;
}

/* Multiply 2 matrices */
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;
} 

/* Calculate the (standard) adjoint of a matrix */
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 FF */
int CheckUnitary(Mtx_t mtx) {
  int ir,ic,OK;
  Mtx_t adj,prod1,prod2;
  AdjMtx(adj,mtx);
  MultMtx(prod1,adj,FF);
  MultMtx(prod2,prod1,mtx);
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) 
      if(cabs(prod2[ir][ic]-FF[ir][ic])>eps2)return FALSE;
  return TRUE;
}  


/* Act on a point in complex hyperbolic 2-space via a matrix */
void ApplyMtx(Pt_t *newpt,Mtx_t mtx,Pt_t pt) {
  double complex f;

  f=1.0/(mtx[2][0]*(pt.z)[0]+mtx[2][1]*(pt.z)[1]+mtx[2][2]);
  (newpt->z)[0]=f*(mtx[0][0]*(pt.z)[0]+mtx[0][1]*(pt.z)[1]+mtx[0][2]);
  (newpt->z)[1]=f*(mtx[1][0]*(pt.z)[0]+mtx[1][1]*(pt.z)[1]+mtx[1][2]);
  (newpt->s2)=PtAbsSq(*newpt);
  return;
}

/* Set up the matrices for the finite subgroup $K$ */
void MakeGpK() {
  int Ix,ir,ic;
  for(ir=0;ir<3;ir++) {
    for(ic=0;ic<3;ic++)
      GpK[0][ir][ic]=0.0;
    GpK[0][ir][ir]=1.0;
  }

  for(Ix=1;Ix<24;Ix++)
    MultMtx(GpK[Ix],GpK[Ix-1],KU);
  for(Ix=24;Ix<48;Ix++)
    MultMtx(GpK[Ix],GpK[Ix-24],KV);
  return;
}
  
/* Display an array of matrices */
void PrintMatrices(Mtx_t Matrices[],int N,char *Name) { 
  int Ix;
  for(Ix=0;Ix<N;Ix++) {
    printf("%s[%i]=\n",Name,Ix);
    PrintMtx(Matrices[Ix]);
  }
  return;
}
  
/* Calculate the scalar product of two points in the complex ball */
inline double complex Scal(Pt_t Pt1,Pt_t Pt2) {
  return Pt1.z[0]*conj(Pt2.z[0])+Pt1.z[1]*conj(Pt2.z[1]);
}

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

inline double PtAbsSq(Pt_t Pt) {
  double a1,a2;
  a1=cabs(Pt.z[0]);
  a2=cabs(Pt.z[1]);
  return a1*a1+a2*a2;
}

inline double QQ(Pt_t Pt1,Pt_t Pt2) {
  double a1;

  a1=cabs(R0*R0-Scal(Pt1,Pt2));
  return (a1*a1)/((R0*R0-PtAbsSq(Pt1))*(R0*R0-PtAbsSq(Pt2)));
}

inline double QQ1(Pt_t Pt1,Pt_t Pt2) {
  double a1;

  a1=cabs(R0*R0-Scal(Pt1,Pt2));
  return (a1*a1)/(R0*R0-PtAbsSq(Pt1));
}
  
/* Compare two points: The order to use is somewhat arbitrary.  Here
   it is lexocographic, with the first number tested being distance
   from the origin. */

int ComparePts(const void *Ptr1,const void *Ptr2) {
  const Pt_t *Pt1=(const Pt_t *)Ptr1, *Pt2=(const Pt_t *)Ptr2;

  if((Pt1->s2)>(Pt2->s2)+eps)return +1;
  if((Pt2->s2)>(Pt1->s2)+eps)return -1;

  if(creal((Pt1->z)[0])>creal((Pt2->z)[0])+eps)return +1;
  if(creal((Pt2->z)[0])>creal((Pt1->z)[0])+eps)return -1;

  if(cimag((Pt1->z)[0])>cimag((Pt2->z)[0])+eps)return +1;
  if(cimag((Pt2->z)[0])>cimag((Pt1->z)[0])+eps)return -1;

  if(creal((Pt1->z)[1])>creal((Pt2->z)[1])+eps)return +1;
  if(creal((Pt2->z)[1])>creal((Pt1->z)[1])+eps)return -1;

  if(cimag((Pt1->z)[1])>cimag((Pt2->z)[1])+eps)return +1;
  if(cimag((Pt2->z)[1])>cimag((Pt1->z)[1])+eps)return -1;

  return 0;
}