#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
// Generator of $\QQ[\sqrt{6},\omega]$ 
// om = -aa^2/2; r = 3*aa^3/(4*om+2)
#define aa  ((om-1)*r/3.0)

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

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

typedef double complex Mtx_t[2][2];

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

// Generators for a big finite subgroup $K$
Mtx_t KU4=
  {{ 1,   0}, 
   { 0, -om}
  };

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

// Other generators of our group $\Gamma$
#define NoGens 1
#define MaxNoGens 1
Mtx_t Gens[MaxNoGens]={
  // This first one is _AA_
  {{      1/3.0*r*om+1/3.0*(2*r+3),            -1/3.0*r*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*(-2*r-3)}}
  };

#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[1];} Pt_t;
Pt_t Pts[PtsSiz]; // Points as generated
Pt_t SavPts[PtsDesired]; // Points from previous list

FILE *PtsFile;
FILE *GnuplotFile;

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 GnuplotOutPts(Pt_t Pts[],int N);
void GnuplotOutPt(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);
  GnuplotFile=fopen(argv[2],"w");
  assert(GnuplotFile != NULL);

  printf("FF=\n");  PrintMtx(FF);
  printf("\nKU4=\n");  PrintMtx(KU4);
  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]));

  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);
  GnuplotOutPts(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)\n",
	 Pt.s2,
	 creal(Pt.z[0]),cimag(Pt.z[0]));
  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;

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

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

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

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

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

/* Output a point in the complex disk to GnuplotFile. */
void GnuplotOutPt(Pt_t Pt) {
  fprintf(GnuplotFile,"% .15f % .15f\n",
	  creal(Pt.z[0]),cimag(Pt.z[0]));
  return;
}

/* Dispay a Matrix */
void PrintMtx(Mtx_t m) {
  int ir,ic;
  for(ir=0;ir<2;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<2;ir++)
    for(ic=0;ic<2;ic++) {
      result[ir][ic]=0;
      for(i=0;i<2;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<2;ir++)
    for(ic=0;ic<2;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<2;ir++)
    for(ic=0;ic<2;ic++) 
      if(cabs(prod2[ir][ic]-FF[ir][ic])>eps2)return FALSE;
  return TRUE;
}  


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

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

/* Set up the matrices for the finite subgroup $K$ */
void MakeGpK() {
  int Ix,ir,ic;
  for(ir=0;ir<2;ir++) {
    for(ic=0;ic<2;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],KU4);
  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 disk */
inline double complex Scal(Pt_t Pt1,Pt_t Pt2) {
  return Pt1.z[0]*conj(Pt2.z[0]);
}

/* Calculate the _real_ scalar product of two points in the complex disk */
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]);
  return a1*a1;
}

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;

  return 0;
}