#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) )

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

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

/* R0 = 1 
   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 R0sq (1)
#define R0 (1)

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

typedef double complex Mtx_t[3][3];

// 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 conjugated matrices of the group are
// unitary
Mtx_t FF=
  {{ 1, 0,     0},
   { 0, 1,     0},
   { 0, 0, -R0sq}};

// Generators for (trivial) finite subgroup $K$

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

// Other generators of our group $\Gamma$
#define NoGens 3
#define MaxNoGens 3
Mtx_t Gens[MaxNoGens];   // These will be the conjugated generators
Mtx_t Gens0[MaxNoGens]={ // These are the given generators
  {{1/18.0*(ZZ*ZZ+3*ZZ+3)*s+1/6.0*(ZZ*ZZ+ZZ-3),  // This is U
      1/126.0*(11*ZZ*ZZ+18.0*ZZ-36)*s+1/6.0*(-ZZ*ZZ-4*ZZ-2),
      1/126.0*(-2*ZZ*ZZ-9*ZZ-3)*s+1/6.0*(ZZ+1)
    },
    {1/63.0*(-2*ZZ*ZZ-9*ZZ-3)*s+1/3.0*(-ZZ-1),
     1/18.0*(ZZ*ZZ+3*ZZ+3)*s+1/6.0*(-ZZ*ZZ-ZZ+1),
     -1/18.0*ZZ*ZZ*s+1/6.0*(-ZZ*ZZ-2*ZZ)
    },
    {1/63.0*(-4*ZZ*ZZ-18.0*ZZ-27)*s+1/3.0*(2*ZZ*ZZ+6*ZZ-1),
     1/63.0*(4*ZZ*ZZ+18.0*ZZ+6)*s,
     1/18.0*(-2*ZZ*ZZ-6*ZZ+3)*s-1/6.0
    }
  },
  {{1/126.0*(2*ZZ*ZZ+9*ZZ+24)*s+1/6.0*(ZZ-2), // This is A
    1/126.0*(-8*ZZ*ZZ-15*ZZ+9)*s+1/6.0*(-ZZ+1),
    1/126.0*(2*ZZ*ZZ+9*ZZ+3)*s+1/6.0*(ZZ+3)
  },
   {1/63.0*(-4*ZZ*ZZ-18*ZZ-6)*s+2/3.0,
    1/126.0*(-4*ZZ*ZZ-18*ZZ+15)*s-1/2.0,
    1/18.0*(ZZ*ZZ+3*ZZ)*s+1/6.0*(ZZ*ZZ+3*ZZ+2)
   },
   {1/63.0*(ZZ*ZZ-6*ZZ-9)*s+1/3.0*(-ZZ*ZZ-2*ZZ+3),
    1/63.0*(2*ZZ*ZZ+9*ZZ+3)*s+1/3.0*(-ZZ+1),
    1/126.0*(2*ZZ*ZZ+9*ZZ+24)*s+1/6.0*(-ZZ-4)
    }
  },
  {{1/126.0*(8*ZZ*ZZ+15*ZZ-9)*s+1/6.0*(ZZ+5), // This is B
    1/63.0*(2*ZZ*ZZ+9*ZZ+3)*s,
    1/252.0*(-11*ZZ*ZZ-18*ZZ-6)*s+1/12.0*(-ZZ*ZZ-2*ZZ+2)
  },
   {1/126.0*(-2*ZZ*ZZ-9*ZZ-24)*s+1/6.0*(2*ZZ*ZZ+5*ZZ),
    1/18.0*(-ZZ*ZZ-3*ZZ)*s+1/6.0*(-ZZ*ZZ-3*ZZ+4),
    1/126.0*(-2*ZZ*ZZ-9*ZZ-3)*s+1/6.0*(ZZ+1)
   },
   {1/63.0*(-2*ZZ*ZZ-9*ZZ-3)*s+1/3.0*(-ZZ-1),
    1/126.0*(13*ZZ*ZZ+27*ZZ-33)*s+1/6.0*(-ZZ*ZZ-3*ZZ+1),
    1/126.0*(-ZZ*ZZ+6*ZZ+9)*s+1/6.0*(ZZ*ZZ+2*ZZ+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 ConjugateMtx(Mtx_t result,Mtx_t m);
void AdjMtx(Mtx_t result,Mtx_t mtx);
int CheckUnitary(Mtx_t mtx, Mtx_t form);
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("FF0=\n");  PrintMtx(FF0);
  printf("FF=\n");  PrintMtx(FF);
  printf("\nU0=\n");  PrintMtx(Gens0[0]);
  printf("A0=\n");  PrintMtx(Gens0[1]);
  printf("B0=\n");  PrintMtx(Gens0[2]);

  // PrintMatrices(Gens0,NoGens,"Gens0");
  for(Ix=0;Ix<NoGens;Ix++)assert(CheckUnitary(Gens0[Ix],FF0));

  for(Ix=0;Ix<NoGens;Ix++)
    ConjugateMtx(Gens[Ix],Gens0[Ix]);
  PrintMatrices(Gens,NoGens,"Gens");

  for(Ix=0;Ix<NoGens;Ix++)
    assert(CheckUnitary(Gens[Ix],FF));

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

  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 conjugate a matrix by ConjMtx */
void ConjugateMtx(Mtx_t result,Mtx_t m) {
  Mtx_t tmp;
  MultMtx(tmp,ConjMtx,m);
  MultMtx(result,tmp,ConjMtxInv);
  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 a sesquilinear form */
int CheckUnitary(Mtx_t mtx, Mtx_t form) {
  int ir,ic,OK;
  Mtx_t adj,prod1,prod2;
  AdjMtx(adj,mtx);
  MultMtx(prod1,adj,form);
  MultMtx(prod2,prod1,mtx);
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) 
      if(cabs(prod2[ir][ic]-form[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;
  }
  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;
}