#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

/* 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 y1Lo (0.0)     // Removes symmetry under conjugation
#define y1Hi (+0.597)
//#define y1Hi (+.57735026918962576451) // sqrt(1/3)
//#define y1Lo (0.31783724)
//#define y1Hi (0.31783725)

#define x2Lo (0.0)     // Takes care of an extra relection in GpK
#define x2Hi (+0.655)
//#define x2Hi (+0.63299316185545206546) // sqrt(8/3)-1
//#define x2Lo (0.44948974)
//define x2Hi (0.44948975)

#define y2Lo (-1.2)
#define y2Hi (+1.2)
//#define y2Lo (-1.09637631717731280408) // sqrt(3)-sqrt(8)
//#define y2Hi (+1.09637631717731280408) // sqrt(8)-sqrt(3)
//#define y2Lo (-0.63567450)
//#define y2Hi (-0.63567449)

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

// #define NumPts 4537
#define NumPts 12

typedef double complex Mtx_t[3][3];

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

#define GpKSiz 48
Mtx_t GpK[GpKSiz];

typedef struct{double complex z[2];} Pt_t;
Pt_t Pts[NumPts];
int PtsCnt[NumPts];

#define KOrbSiz 24
Pt_t KOrb[KOrbSiz];

int ReadPts(Pt_t Pts[]);
void PrintPts(Pt_t Pts[],int N);
void PrintPt(Pt_t Pt);
void PrintMtx(Mtx_t m);
void MultMtx(Mtx_t result,Mtx_t m1,Mtx_t m2);
void ApplyMtx(Pt_t *newpt,Mtx_t mtx,Pt_t pt);
void MakeGpK();
void PrintGpK();
void MakeKOrb();
void PrintKOrb();
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 main() {
  int Ix,y1Ix,x2Ix,y2Ix,OK,Cnt;
  double Vol,VolCpt,hhy1,hhx2,hhy2,hh0,
    sMax,ss,tt,VolFact,fTmp,
    y1,x2,y2,TestScal,y1Max,y1Min,x2Max,x2Min,y2Max,y2Min,
    Q,NormFact,JacFact,
    y1AtMax,x2AtMax,y2AtMax;
  Pt_t Pt0;

  printf("KU=\n");  PrintMtx(KU);
  printf("KV=\n");  PrintMtx(KV);

  MakeGpK();
  // PrintGpK();
  MakeKOrb();
  PrintKOrb();

  assert(ReadPts(Pts)==NumPts);
  // PrintPts(Pts,NumPts);

  memset((void *) PtsCnt,(int) 0,(size_t) sizeof(PtsCnt));
  sMax=0.0;
  Vol=0.0;
  VolCpt=0.0;

  hhy1=((y1Hi-y1Lo)/NN);
  hhx2=((x2Hi-x2Lo)/NN);
  hhy2=((y2Hi-y2Lo)/NN);
  hh0=hhy1*hhx2*hhy2;

  y1Min=y1Hi; x2Min=x2Hi; y2Min=y2Hi; 
  y1Max=y1Lo; x2Max=x2Lo; y2Max=y2Lo; 
  Cnt=0;

  /* Run through all possible directions on the 3-sphere. */

  for(y1Ix=0;y1Ix<NN;y1Ix++) {
    y1=y1Lo+(y1Ix+0.5)*hhy1;
    for(x2Ix=0;x2Ix<NN;x2Ix++) {
      x2=x2Lo+(x2Ix+0.5)*hhx2;
      for(y2Ix=0;y2Ix<NN;y2Ix++) {
	y2=y2Lo+(y2Ix+0.5)*hhy2;
	Pt0.z[0]=1.0+y1*I;
	Pt0.z[1]=x2+y2*I;
	/* Check whether the direction chosen lies in the
	   fundamental domain for the finite group $K$ acting
           on the 3-sphere. */
	TestScal=RealScal(Pt0,KOrb[0]);
	OK=TRUE;
	for(Ix=1;Ix<KOrbSiz;Ix++) {
	  if(RealScal(Pt0,KOrb[Ix])>TestScal) {
	    OK=FALSE;
	    break;
	  }
	}
	if(OK) {
	  y1Min=(y1<y1Min ? y1 : y1Min);
	  x2Min=(x2<x2Min ? x2 : x2Min);
	  y2Min=(y2<y2Min ? y2 : y2Min);
	  y1Max=(y1>y1Max ? y1 : y1Max);
	  x2Max=(x2>x2Max ? x2 : x2Max);
	  y2Max=(y2>y2Max ? y2 : y2Max);
	  Cnt++;

	  Q=1.0/(1.0+y1*y1+x2*x2+y2*y2);
	  NormFact=sqrt(Q);
	  JacFact=Q*Q;
	  Pt0.z[0] *= NormFact; Pt0.z[1] *= NormFact;
	   
	  ss=sVal(Pt0);
	  if(ss>sMax) {
	    sMax=ss;
	    y1AtMax=y1; x2AtMax=x2; y2AtMax=y2;
	  }
	  tt=ss/R0;
	  fTmp=tt*tt/(1.0-tt*tt);
	  Vol+=hh0*JacFact*fTmp*fTmp/4.0;
	  VolCpt+=hh0*JacFact/4.0;
	}
      }
    }      
  }

  printf("y1Min,y1Max= %20.15f,%20.15f\n",y1Min,y1Max);
  printf("x2Min,x2Max= %20.15f,%20.15f\n",x2Min,x2Max);
  printf("y2Min,y2Max= %20.15f,%20.15f\n",y2Min,y2Max);
  printf("Fraction of directions in the GpK-fundamental domain=%f\n\n",
	 (double)(Cnt)/(double)(NN*NN*NN));

  Vol*=2.0; // Because we've divided out by complex conjugation
  VolCpt*=2.0;

  for(Ix=0;Ix<NumPts;Ix++) {
    if(PtsCnt[Ix]>0) {
      printf("Ix=%6i Cnt=%10i Pt=",Ix,PtsCnt[Ix]);
      PrintPt(Pts[Ix]);
    }
  }
    
  printf("y1AtMax, x2AtMax, y2AtMax=%.15f, %.15f, %.15f\n",y1AtMax, x2AtMax, y2AtMax);
  printf("VolCpt=%.15f Normalized VolCpt=%.15f\n",
  	 VolCpt,2.0*VolCpt/(M_PI*M_PI),sMax);
  printf("Vol=%.15f Normalized Vol=%.15f sMax=%.15f\n",
  	 Vol,2.0*Vol/(M_PI*M_PI),sMax);
}

int ReadPts(Pt_t Pts[]) {
  int Ix;
  double z1R,z1I,z2R,z2I,s2;
  FILE *PtFile;
  // PtFile=fopen("C18a-5000SelectedPts","r");
  PtFile=fopen("C18a-12SelectedPts","r");
  assert(PtFile != NULL);

  Ix=0;
  while(fscanf(PtFile,"%*[ \n\t]") != EOF) {
    assert(fscanf(PtFile,"%lf (%lf + %lf*i, %lf + %lf*i) ",
		  &s2,&z1R,&z1I,&z2R,&z2I)==5);
    assert(Ix<NumPts);
    Pts[Ix].z[0]=z1R+z1I*I;
    Pts[Ix].z[1]=z2R+z2I*I;
    Ix++;
  }
  return Ix;
}

/* Output 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",
	 PtAbsSq(Pt),
	 creal(Pt.z[0]),cimag(Pt.z[0]),
	 creal(Pt.z[1]),cimag(Pt.z[1]));
  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;
} 

/* 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]);
  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=0;Ix<24;Ix++)
    MultMtx(GpK[Ix+24],GpK[Ix],KV);
  return;
}
  
/* Display the matrices for the finite subgroup $K$ */
void PrintGpK() { 
  int Ix;
  for(Ix=0;Ix<GpKSiz;Ix++) {
    printf("GpK[%i]=\n",Ix);
    PrintMtx(GpK[Ix]);
  }
  return;
}
  
/* Set up an orbit of the finite subgroup $K$ */
void MakeKOrb() {
  int GpKIx,OrbIx,OrbIx2,OK;
  (KOrb[0].z)[0]=1.0;
  (KOrb[0].z)[1]=0.0;

  OrbIx=1;
  for(GpKIx=1;GpKIx<GpKSiz && OrbIx<KOrbSiz;GpKIx++) {
    ApplyMtx(&KOrb[OrbIx],GpK[GpKIx],KOrb[0]);
    OK=TRUE;
    for(OrbIx2=0;OrbIx2<OrbIx;OrbIx2++)
      if(cabs((KOrb[OrbIx].z)[0]-(KOrb[OrbIx2].z)[0])<eps &&
	 cabs((KOrb[OrbIx].z)[1]-(KOrb[OrbIx2].z)[1])<eps) {
	break;
      }
    if(OK) {
      OrbIx++;
    }
  }
  assert(OrbIx == KOrbSiz);
}

/* Display an orbit of the finite subgroup $K$ */
void PrintKOrb() {
  int Ix;
  for(Ix=0;Ix<KOrbSiz;Ix++) {
    printf("KOrb[%02i]=",Ix);
    PrintPt(KOrb[Ix]);
  }
  return;
}

/* Calculate how far the Dirichlet fundamental domain extends along
   the ray through a given point of norm 1 in $\CC^2$. */

double sVal(Pt_t Pt0) {
  int Ix,Ix0;
  double complex u1;
  double R,limit,real1,a1,a2,zAbsSq,test,s,s0,QQSav,QQ0,QQz;
  Pt_t LPt,ZeroPt;

  assert(cabs(Scal(Pt0,Pt0)-1.0)<eps);

  ZeroPt.z[0]=ZeroPt.z[1]=0.0;
  LPt.z[0]=LPt.z[1]=0.0;
  s0=limit=HUGE_VAL;
  QQSav=HUGE_VAL;

  Ix0=(-1);
  for(Ix=0;Ix<NumPts;Ix++) {
    // printf("sVal 100: Pt0="); PrintPt(Pt0);
    // printf("          z="); PrintPt(Pts[Ix]);
    if(QQ1(Pts[Ix],LPt)<limit) {
      assert(QQ(Pts[Ix],LPt)-QQSav < eps*QQSav);
      u1=Scal(Pt0,Pts[Ix]);
      // printf("sVal 200: u1=% .15f%+.15f\n",u1);
      real1=creal(u1);
      // printf("sVal 300: real1=% .15f\n",real1);
      if(real1>0) {
	a1=cabs(u1);
	// printf("sVal 400:  a1=% .15f\n",a1);
	zAbsSq=PtAbsSq(Pts[Ix]);
	// printf("sVal 500: zAbsSq=% .15f\n",zAbsSq);
	test=real1*real1-(a1*a1*zAbsSq)/(R0*R0);
	// printf("sVal 600: test=% .15f\n",test);
	if(test>0) {
	  s=zAbsSq/(real1+sqrt(test));
	  // printf("sVal 700: Pt0="); PrintPt(Pt0);
	  // printf("              z="); PrintPt(Pts[Ix]);
	  // printf("              s=% .15f\n",s);
	  if(s<R0) {
	    Ix0=Ix;
	    assert(s<s0+eps2);
	    s0=s;
	    LPt.z[0]=s*Pt0.z[0];
	    LPt.z[1]=s*Pt0.z[1];
	    QQ0=QQ(ZeroPt,LPt);
            QQz=QQ(Pts[Ix],LPt);
	    assert(fabs((QQz-QQ0)/QQ0)<eps);
	    // printf("QQSav=%15.10f QQ=%15.10f Ix0=%4i\n",
	    // QQSav,QQ(ZeroPt,LPt),Ix0);
	    QQSav=QQ0;
	    limit=QQ1(Pts[Ix],LPt);
	  }  
	}
      }
    }
  }
  // assert(Ix0 != -1);
  if(Ix0 == -1) {
    printf("Bad Pt0: ");
    PrintPt(Pt0);
    printf("           (% .5f%+.5f*i,% .5f%+.5f*i) "
	   "args=(% .5f,% .5f) abs=(%.5f,%.5f)\n",
   	   creal(Pt0.z[0]),cimag(Pt0.z[0]),
   	   creal(Pt0.z[1]),cimag(Pt0.z[1]),
   	   carg(Pt0.z[0])/M_PI,carg(Pt0.z[1])/M_PI,
   	   cabs(Pt0.z[0]),cabs(Pt0.z[1]));
  }
  PtsCnt[Ix0]++;
  // printf("Ix0=%i s0=%f\n",Ix0,s0);
  return s0;
}

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