#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 K1=
  {{-om, 0, 0},
   {  0, 1, 0},
   {  0, 0, 1}};

Mtx_t K2=
  {{1,   0, 0},
   {0, -om, 0},
   {0,   0, 1}};

Mtx_t K3=
  {{ 0,1, 0},
   { 1,0, 0},
   { 0,0, 1}};

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

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

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

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

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

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

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

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

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

  {{0, 1/3.0*(-r-3)*om+1/3.0*(-2*r-3), 1/3.0*r*om+2.0/3*r},
   {1/3.0*(2*r+3)*om+1/3.0*r, r+2, 1/3.0*r*om+1/3.0*(-r-6)},
   {1/3.0*(r+3)*om+1/3.0*(-r-3), 1/3.0*(4*r+9)*om+1/3.0*(5*r+12),
      1/3.0*(-2*r-6)*om+1/3.0*(-4*r-9)}},

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

#define PtsDesired 30000 // (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("\nK1=\n");  PrintMtx(K1);
  printf("K2=\n");  PrintMtx(K2);
  printf("K3=\n");  PrintMtx(K3);

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

  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<6;Ix++)
    MultMtx(GpK[Ix],GpK[Ix-1],K1);
  for(Ix=6;Ix<36;Ix++)
    MultMtx(GpK[Ix],GpK[Ix-6],K2);
  for(Ix=36;Ix<72;Ix++)
    MultMtx(GpK[Ix],GpK[Ix-36],K3);

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