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

#define phi (1.61803398874989484)
#define R0  (1.27201964951406896)

#define Zeta10       (0.8090169943749474+0.5877852522924731*I)
#define Zeta10Sq     (0.3090169943749474+0.9510565162951535*I)
#define Zeta10Inv    (0.8090169943749474-0.5877852522924731*I)
#define Zeta10SqInv  (0.3090169943749474-0.9510565162951535*I)
#define Zeta20       (0.9510565162951535+0.3090169943749474*I)
#define Zeta20Cub    (0.5877852522924731+0.8090169943749474*I)
#define Zeta20Inv    (0.9510565162951535-0.3090169943749474*I)
#define Zeta20CubInv (0.5877852522924731-0.8090169943749474*I)

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

#define NumPts 6368

typedef struct{double complex z[2];} Pt_t;
Pt_t Pts[NumPts];
Pt_t NPts[NumPts];
int NumAllPts,AllPtsCnt[NumPts*200];
Pt_t AllPts[NumPts*200];

int ReadPts(Pt_t Pts[]);
void PrintPts(Pt_t Pts[],int N);
void PrintPt(Pt_t Pt);
inline void NormalizePt(Pt_t p1,Pt_t *p1Norm);
inline double complex NormalizeZ(double complex z);
void MakeAllPts();
double sVal(Pt_t Theta);
inline complex Scal(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,Th1Ix,Th2Ix,Th3Ix;
  double Th1,Th2,Th3,r1,r2,Vol,hh0,hh1,hh2,hh3,
    sMax,ss,tt,VolFact,fTmp;
  Pt_t Theta;

  // These lines are for initial runs, where we wish to identify
  // the relevant points in the orbit of zero.
  //
  // assert(ReadPts(Pts)==NumPts);
  // PrintPts(Pts,NumPts);
  // for(Ix=0;Ix<NumPts;Ix++) {
  //   NormalizePt(Pts[Ix],&NPts[Ix]);
  // }
  // MakeAllPts();

  // These lines are for when we have already identified the
  // relevant points in the orbit of zero.
  //
  NumAllPts=ReadPts(AllPts);
  PrintPts(AllPts,NumAllPts);

  memset(&AllPtsCnt,0,sizeof(AllPtsCnt));
  sMax=0.0;
  Vol=0.0;
  hh0=(M_PI/(4.0*NTheta))*(M_PI/(5.0*NTheta))*(M_PI/(5.0*NTheta));

  Th3=0.0;
  for(Th3Ix=0;Th3Ix<=NTheta;Th3Ix++) {
    if(Th3Ix==0 || Th3Ix==NTheta)
      hh1=hh0*0.5;
    else
      hh1=hh0;
    r1=cos(Th3);
    r2=sin(Th3);
    // Th1=(-M_PI/10.0);
    Th1=0.0;
    for(Th1Ix=0;Th1Ix<=NTheta;Th1Ix++) {
      if(Th1Ix==0 || Th1Ix==NTheta)
	hh2=hh1*0.5;
      else
	hh2=hh1;
      Theta.z[0]=r1*(cos(Th1)+I*sin(Th1));
      // Th2=(-M_PI/10.0);
      Th2=0.0;
      for(Th2Ix=0;Th2Ix<=NTheta;Th2Ix++) {
	if(Th2Ix==0 || Th2Ix==NTheta)
	  hh3=hh2*0.5;
	else
	  hh3=hh2;
	Theta.z[1]=r2*(cos(Th2)+I*sin(Th2));
	ss=sVal(Theta);
	if(ss>sMax)sMax=ss;
	tt=ss/R0;
	fTmp=tt*tt/(1.0-tt*tt);
	Vol+=hh3*r1*r2*fTmp*fTmp/4.0;
	// Th2+=(M_PI/5.0)/NTheta;
	Th2+=(M_PI/10.0)/NTheta;
      }
      // Th1+=(M_PI/5.0)/NTheta;
      Th1+=(M_PI/10.0)/NTheta;
    }
    Th3+=(M_PI/4.0)/NTheta;
  }

  for(Ix=0;Ix<NumAllPts;Ix++) {
    if(AllPtsCnt[Ix]>0) {
      printf("Ix=%6i Cnt=%6i Pt=",Ix,AllPtsCnt[Ix]);
      PrintPt(AllPts[Ix]);
    }
  }
  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;
  FILE *PtFile;
  PtFile=fopen("Pts6368Selected.txt","r");
  assert(PtFile != NULL);

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

void PrintPts(Pt_t Pts[],int N) {
  int Ix;

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

void PrintPt(Pt_t Pt) {
  printf("(% .15f%+.15f*i,% .15f%+.15f*i)\n",
	 creal(Pt.z[0]),cimag(Pt.z[0]),
	 creal(Pt.z[1]),cimag(Pt.z[1]));
  return;
}

inline void NormalizePt(Pt_t p1,Pt_t *p1Norm) {

  if(cabs(p1.z[0]) >= cabs(p1.z[1])) {
    p1Norm->z[0]=NormalizeZ(p1.z[0]);
    p1Norm->z[1]=NormalizeZ(p1.z[1]);
  }
  else {
    p1Norm->z[0]=NormalizeZ(p1.z[1]);
    p1Norm->z[1]=NormalizeZ(p1.z[0]);
  }
}  

inline double complex NormalizeZ(double complex z) {
  double complex zN;
  double test;

  if(fabs(creal(z))<eps) {
    if(cimag(z)>0)
      zN=z*Zeta10SqInv;
    else
      zN=z*Zeta10Sq;
  }
  else {
    if(creal(z)>0)
      zN=z;
    else
      zN=(-z);

    test=cimag(zN)/creal(zN);
    if(test>0) {
      if(test>cimag(Zeta20Cub)/creal(Zeta20Cub))
	zN*=Zeta10SqInv;
      else
	if(test>cimag(Zeta20)/creal(Zeta20))
	  zN*=Zeta10Inv;
    }
    else {
      if(test<cimag(Zeta20CubInv)/creal(Zeta20CubInv))
	zN*=Zeta10Sq;
      else
	if(test<cimag(Zeta20Inv)/creal(Zeta20Inv))
	  zN*=Zeta10;
    }
  }
  return zN;
}
    
void MakeAllPts() {
  int Ix,TenIx1,TenIx2,EqFg;
  double complex z1,z2;

  NumAllPts=0;
  for(Ix=1;Ix<NumPts;Ix++) {
    if(cabs(NPts[Ix].z[1])<eps) {
      assert(cabs(NPts[Ix].z[0])>eps);
      z1=NPts[Ix].z[0];
      for(TenIx1=0;TenIx1<10;TenIx1++) {
	AllPts[NumAllPts].z[0]=z1;
	AllPts[NumAllPts].z[1]=0;
	NumAllPts++;
	AllPts[NumAllPts].z[0]=0;
	AllPts[NumAllPts].z[1]=z1;
	NumAllPts++;
	z1*=Zeta10;
      }
    }
    else {
      EqFg=((cabs(NPts[Ix].z[0]-NPts[Ix].z[1])<eps));
      z1=NPts[Ix].z[0];
      for(TenIx1=0;TenIx1<10;TenIx1++) {
	z2=NPts[Ix].z[1];
	for(TenIx2=0;TenIx2<10;TenIx2++) {
	  AllPts[NumAllPts].z[0]=z1;
	  AllPts[NumAllPts].z[1]=z2;
	  NumAllPts++;
	  if(!EqFg) {
	    AllPts[NumAllPts].z[0]=z2;
	    AllPts[NumAllPts].z[1]=z1;
	    NumAllPts++;
	  }
	  z2*=Zeta10;
	}
	z1*=Zeta10;
      }
    }
  }
  // for(Ix=0;Ix<NumAllPts;Ix++) {
  //   printf("(% .5f%+.5f*i,% .5f%+.5f*i) args=(% .5f,% .5f) abs=(%.5f,%.5f)\n",
  // 	   creal(AllPts[Ix].z[0]),cimag(AllPts[Ix].z[0]),
  // 	   creal(AllPts[Ix].z[1]),cimag(AllPts[Ix].z[1]),
  // 	   carg(AllPts[Ix].z[0])/M_PI,carg(AllPts[Ix].z[1])/M_PI,
  // 	   cabs(AllPts[Ix].z[0]),cabs(AllPts[Ix].z[1]));
  // }
}

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

  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<NumAllPts;Ix++) {
    // printf("sVal 100: Theta="); PrintPt(Theta);
    // printf("          z="); PrintPt(AllPts[Ix]);
    if(QQ1(AllPts[Ix],LPt)<limit) {
      assert(QQ(AllPts[Ix],LPt)-QQSav < eps*QQSav);
      u1=Scal(Theta,AllPts[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(AllPts[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: Theta="); PrintPt(Theta);
	  // printf("              z="); PrintPt(AllPts[Ix]);
	  // printf("              s=% .15f\n",s);
	  if(s<R0) {
	    Ix0=Ix;
	    assert(s<s0+eps2);
	    s0=s;
	    LPt.z[0]=s*Theta.z[0];
	    LPt.z[1]=s*Theta.z[1];
	    QQ0=QQ(ZeroPt,LPt);
            QQz=QQ(AllPts[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(AllPts[Ix],LPt);
	  }  
	}
      }
    }
  }
  // assert(Ix0 != -1);
  if(Ix0 == -1) {
    printf("Bad Theta: ");
    PrintPt(Theta);
    printf("           (% .5f%+.5f*i,% .5f%+.5f*i) "
	   "args=(% .5f,% .5f) abs=(%.5f,%.5f)\n",
   	   creal(Theta.z[0]),cimag(Theta.z[0]),
   	   creal(Theta.z[1]),cimag(Theta.z[1]),
   	   carg(Theta.z[0])/M_PI,carg(Theta.z[1])/M_PI,
   	   cabs(Theta.z[0]),cabs(Theta.z[1]));
  }
  AllPtsCnt[Ix0]++;
  // printf("Ix0=%i s0=%f\n",Ix0,s0);
  return s0;
}

inline complex Scal(Pt_t Pt1,Pt_t Pt2) {
  return Pt1.z[0]*conj(Pt2.z[0])+Pt1.z[1]*conj(Pt2.z[1]);
}

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