#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 sqm7 (2.645751311064590590501615753639*I)

// 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)
#define NN (400)

#define NumPts 161

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 (trivial) finite subgroup $K$

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

FILE *PtsFile;

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

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 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 main(int argc,char *(argv[])) {
  int Ix,Th1Ix,Th2Ix,Th3Ix;
  double Vol,VolCpt,Th1,Th2,Th3,r1,r2,hh0,hh1,hh2,hh3,
    sMax,ss,tt,VolFact,fTmp,
    Q,JacFact;
  Pt_t Pt0;

  MakeGpK();
  // PrintMatrices(GpK,GpKSiz,"GpK");
  for(Ix=0;Ix<GpKSiz;Ix++)assert(CheckUnitary(GpK[Ix]));

  PtsFile=fopen(argv[1],"r");
  assert(PtsFile != NULL);
  assert(ReadPts(Pts)==NumPts);

  // PrintPts(Pts,NumPts);

  memset(PtsCnt,0,sizeof(PtsCnt));
  sMax=0.0;
  Vol=0.0;
  VolCpt=0.0;
  hh0=((M_PI/2.0)/NN)*((2.0*M_PI)/NN)*((2.0*M_PI)/NN);

  /* Run through all possible directions on the 3-sphere. */
  Th3=0.0;
  for(Th3Ix=0;Th3Ix<=NN;Th3Ix++) {
    if(Th3Ix % 3 == 0)
      printf("A: Th3Ix=%i Th3/(pi/2)=%f sMax=%f\n",Th3Ix,
	     Th3/(M_PI/2.0),sMax);
    if(Th3Ix==0 || Th3Ix==NN)
      hh1=hh0*0.5;
    else
      hh1=hh0;
    r1=cos(Th3);
    r2=sin(Th3);
    JacFact=r1*r2;
//Th1=0.0;  // Here we use symmetry under complex conjugation
    Th1=(-M_PI);
    for(Th1Ix=0;Th1Ix<=NN;Th1Ix++) {
      // printf("B: Th3Ix=%i Th3=%f Th1Ix=%i Th1=%f sMax=%f\n",
      //     Th3Ix,Th3,Th1Ix,Th1,sMax);
      if(Th1Ix==0 || Th1Ix==NN)
	hh2=hh1*0.5;
      else
	hh2=hh1;
      Pt0.z[0]=r1*(cos(Th1)+I*sin(Th1));
      Th2=(-M_PI);
      for(Th2Ix=0;Th2Ix<=NN;Th2Ix++) {
	if(Th2Ix==0 || Th2Ix==NN)
	  hh3=hh2*0.5;
	else
	  hh3=hh2;
	Pt0.z[1]=r2*(cos(Th2)+I*sin(Th2));
	ss=sVal(Pt0);
	sMax=MAX(ss,sMax);
	tt=ss/R0;
	fTmp=tt*tt/(1.0-tt*tt);
	Vol+=hh3*JacFact*fTmp*fTmp/4.0;
	VolCpt+=hh3*JacFact/4.0;
	Th2+=(2.0*M_PI)/NN;
      }
      Th1+=(2.0*M_PI)/NN;
    }
    Th3+=(M_PI/2.0)/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=%5i Cnt=%10i Pt=",Ix,PtsCnt[Ix]);
      PrintPt(Pts[Ix]);
    }
  }
    
  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 z0R,z0I,z1R,z1I,s2;

  assert(PtsFile != NULL);
  Ix=0;
  while(fscanf(PtsFile,"%*[ \n\t]") != EOF) {
    assert(fscanf(PtsFile,"%lf (%lf + %lf*i, %lf + %lf*i) ",
		  &s2,&z0R,&z0I,&z1R,&z1I)==5);
    assert(Ix<NumPts);
    Pts[Ix].s2=s2;
    Pts[Ix].z[0]=z0R+z0I*I;
    Pts[Ix].z[1]=z1R+z1I*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",
	 Pt.s2,
	 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;
} 

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