#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <limits.h>
#include <math.h>
#include <complex.h>
#include <string.h>
#include <gsl/gsl_poly.h>

#include "../numeric/SU21-numeric.h"
#include "../numeric/SU21-numeric-base.c"

#define epsHuge (1E-3)

// Type for a complex hyperbolic triangle
//
// Suppose that the three vertices are given by vectors~$A$, $B$,
// and~$C$ in a hermitian space of signature $(-1,-1,+1)$, satisfying
// $(A,A)=(B,B)=(C,C)=1$.
//
// Then the invariants of the triangle can be expressed in terms of
// the matrix of inner products:
//
//   (A,A) (A,B) (A,C)      1   (A,B)  (A,C)
//   (B,A) (B,B) (B,C)  = (B,A)   1    (B,C)
//   (C,A) (C,B) (C,C)    (C,A) (C,B)    1
//
// The three side lengths, $a$, $b$, and~$c$, are related to this
// matrix via:
//
//   \cosh a = |(B,C)|     \cosh b = |(C,A)|     \cosh c = |(A,B)|
//
// A fourth invariant is Brehm's shape, $\RE (A,B)(B,C)(C,A)$.
// 
// Instead of the side lengths, one can deal with their squared
// hyperbolic cosines, $|(B,C)|^2$, $|(C,A)|^2$, $|(A,B)|^2$, or the
// the symmetric functions of the squared hyperbolic cosines.
//
// Alternatively, one can use the squared hyperbolic sines,
// $|(B,C)|^2-1$, $|(C,A)|^2-1$, $|(A,B)|^2-1$, or the symmetric
// functions in the squared hyperbolic sines.
//
// Related invariants can be calculated starting with the inverse of
// the above matrix.  Done naively, this excludes the case when the
// triangle is contained in a complex hyperbolic disc, for in that
// case the matrix of inner products is singular.
//

typedef struct {
  FLOAT_t Side[3];  // The 3 side-lengths
  FLOAT_t ShapeL1; // Brehm's shape LESS 1 --- 1 LESS than the real
	   	   // part of the product of the three inner products
	   	   // taken in cyclic order: (A,B)(B,C)(C,A) - 1
  FLOAT_t ShapeIm; // The imaginary part of (A,B)(B,C)(C,A)
  FLOAT_t Det; // Determinant of the matrix of inner products
  FLOAT_t SymSinh2[3]; // The symmetric functions of the squared
                       // sinh's of the side lengths.
  FLOAT_t SymCosh2[3]; // The symmetric functions squared hyperbolic
		       // cosh's of the side lengths; that means the
		       // symmetric functions in the squared absolute
		       // values of the off-diagonal elements of the
		       // matrix of inner products.
  FLOAT_t SymInv[3]; // The symmetric functions of the squared
		     // absolute values of the off-diagonal elements
		     // of the inverse of the matrix of inner products
  FLOAT_t TrInv; // The trace of the inverse of the matrix of inner
		 // products
  FLOAT_t InvShape; // The real part of the product of the cyclical
		    // product of the three off-diagonal elements of
		    // the inverse of the matrix of inner products
  FLOAT_t InvShapeIm; // The imaginary part of the product of the
		      // cyclical product of the three off-diagonal
		      // elements of the inverse of the matrix of
		      // inner products
  unsigned char Flags;
#define REALFG 01  // Lies in a totally real hyperbolic disc
#define CPXFG  02  // Lies in a complex hyperbolic disc
#define ISOFG  04  // Is isosceles
#define NEARCPXFG 010 // ``Nearly'' lies in a complex hyperbolic disc

} CHypTri_t;

void CHypTriFromSidesAndShape(CHypTri_t *Tri_p);
void CHypTriFromPts(CHypTri_t *Tri_p,
		    const Pt_t PtA,const Pt_t PtB,const Pt_t PtC);
void CheckCHypTri(const CHypTri_t *Tri_p);
void CHypTriResultant(const CHypTri_t *Tri_p);
void CHypTriFindU(FLOAT_t *uSol,const CHypTri_t *Tri_p);
void CHypTriSinhR2(FLOAT_t *SinhR2,int *RealCnt_p,
		   const FLOAT_t *uSol,const CHypTri_t *Tri_p);
FLOAT_t CircumRadius(Pt_t X,Pt_t A,Pt_t B,Pt_t C);
int gPolish(COMPLEX_t *g0_p,COMPLEX_t *g1_p,COMPLEX_t *g2_p,
	     const COMPLEX_t v12,const COMPLEX_t v20,const COMPLEX_t v01);
void gAdd(int *gNum_p,COMPLEX_t gList[6][3],const int Ignore,
	  COMPLEX_t g[3],COMPLEX_t v[3][3]);
void uCheck(FLOAT_t uEq[],FLOAT_t uEqComp[],COMPLEX_t uVal[]);

#define NN 20
//#define RR 3.0
#define RR 40.0

int GoodCnt=0, BadCnt=0, gNumCnt[7], NumTorusCnt[7];
FLOAT_t sigmaMaxErr[6]={0,0,0,0,0,0};

int main(int argc,char *(argv[])) {
  char *Prefix_p,*NN_p,FileNam[MAX_FILENAME_LEN]="../";
  // Input files
  FILE *FewEltsFile;
  // Output files

  int NumFewElts,NumPts,Pt1Ix,Pt2Ix;
  Pt_t *Pts,DummyPt;
  Elt_t *FewElts;

  CHypTri_t Tri,*Tri_p=&Tri;
  FLOAT_t uSol[12],*uR,*uI,SinhR2[6],SinhR2C,RRelErr,MaxRRelErr=0;
  int RealCnt,PosCnt,Ix,BadFg,SavIx;

  // Get the prefix
  Prefix_p=getenv("PREFIX");
  assert(Prefix_p != NULL);

  // Open the files
  {
    size_t len=strlen(Prefix_p);
    assert(3+len+9<MAX_FILENAME_LEN);
    strcpy(&FileNam[3],Prefix_p);
    // Input files

    strcpy(&FileNam[3]+len,"-FewElts");
    printf("%s\n",FileNam);
    FewEltsFile=fopen(FileNam,"r");
    assert(FewEltsFile != NULL);
    fflush(stdout);
  }

  // Get the few elements
  GetEltsString(FewEltsFile,"NumFewElts",&NumFewElts,&FewElts);
  qsort(FewElts,(size_t) NumFewElts,sizeof(Elt_t),CompareEltsWithWords);

  SetupPts(&Pts,&NumPts,FewElts,NumFewElts);
  printf("NumFewElts=%i NumPts=%i\n",NumFewElts,NumPts);
  fflush(stdout);

  memset(gNumCnt,0,sizeof(gNumCnt));
  memset(NumTorusCnt,0,sizeof(NumTorusCnt));

  for(Pt1Ix=1;Pt1Ix<NumPts;Pt1Ix++) { // <===
    printf("CHypTri P: Pt1Ix=%5i\n",Pt1Ix);
    for(Pt2Ix=Pt1Ix+1;Pt2Ix<NumPts;Pt2Ix++) {
      SinhR2C=CircumRadius(DummyPt,ZeroPt,Pts[Pt1Ix],Pts[Pt2Ix]);
      if(SinhR2C < HUGE_FLOAT) {
	CHypTriFromPts(Tri_p,ZeroPt,Pts[Pt1Ix],Pts[Pt2Ix]);
	//if(!(Tri_p->Flags & (NEARCPXFG | REALFG | ISOFG))) {
	if(TRUE) {
	  BadFg=FALSE;
	  CHypTriResultant(Tri_p);
	  CHypTriFindU(uSol,Tri_p);
	  CHypTriSinhR2(SinhR2,&RealCnt,uSol,Tri_p);
	  //printf("%i",RealCnt);
	  PosCnt=0;
	  for(Ix=0;Ix<RealCnt;Ix++)
	    if(SinhR2[Ix]>0) PosCnt++;
	  if(PosCnt != 1) {
	    BadFg=TRUE;
	    printf("CHypTri Q: PosCnt=%i\n",PosCnt);
	  }
	  //if(PosCnt>1 || (RealCnt!=4 && RealCnt!=6)) {
	  //  uR=&uSol[0];
	  //  uI=&uSol[1];
	  //  printf("\nPosCnt=%1i RealCnt=%1i PtIx1=%i PtIx2=%i\n",
	  //	 PosCnt,RealCnt,Pt1Ix,Pt2Ix);
	  //  for(Ix=0;Ix<6;Ix++) {
	  //    printf("% 11.5g%+11.5g*i\n",*uR,*uI);
	  //    uR+=2; uI+=2;
	  //  }
	  //  printf("hi/h=%11.5g "
	  //	 "d/hL1=%11.5g\n",
	  //	 Tri_p->ShapeIm/(1.0+Tri_p->ShapeL1),
	  //	 Tri_p->Det/Tri_p->ShapeL1);
	  //}
	  //if(PosCnt==0)
	  //  printf("N\n");
	  //else {
	  //  printf("O\n");
	  //}
	  for(Ix=0;Ix<RealCnt;Ix++)
	    if(SinhR2[Ix]>0) {
	      SavIx=Ix;
	      RRelErr=FABS(SinhR2[Ix]-SinhR2C)/SinhR2C;
	      if(RRelErr > MaxRRelErr) {
		BadFg=TRUE;
		MaxRRelErr=RRelErr;
	      }
	      if(BadFg || RRelErr > 0.1) {
		printf("CHypTri E: RRelErr=%12.6g "
		       "MaxRRelErr=%12.6g\n",
		       RRelErr,MaxRRelErr);
		BadFg=TRUE;
	      }
	    }
	  if(BadFg) {
		printf("CHypTri S: Side lengths: %14.8g %14.8g %14.8g\n",
		       Tri_p->Side[0],Tri_p->Side[1],Tri_p->Side[2]);
		printf("  SinhR2C=%23.16g  R=%14.8g\n",
		       SinhR2C,ASINH(SQRT(SinhR2C)));
		printf("   SinhR2=%23.16g  R=%14.8g\n",
		       SinhR2[SavIx],ASINH(SQRT(SinhR2[SavIx])));
	  }
	}
      }
    }
    printf("CHypTri C: GoodCnt=%10i BadCnt=%10i\n",GoodCnt,BadCnt);
    printf("      gNumCnt=");
    for(Ix=1;Ix<=6;Ix++)
      printf("%12i",gNumCnt[Ix]);
    printf("\n");
    printf("  NumTorusCnt=");
    for(Ix=1;Ix<=6;Ix++)
      printf("%12i",NumTorusCnt[Ix]);
    printf("\n");
  }
}

//int main(int argc,char *(argv[])) {
//  CHypTri_t Tri,*Tri_p=&Tri;
//  int PerIx,s0Ix,s1Ix,s2Ix,ShapeIx,j;
//  FLOAT_t SumSideIx,Per,ShapeMax,ShapeMin;
//  FLOAT_t Cosh[3],Sinh2[3],SinhTmp;
//
//  FLOAT_t uSol[12],*solR,*solI,SinhR2[6];
//  int RealCnt,PosCnt,Ix;
//
//  for(PerIx=1;PerIx<=NN;PerIx++) {
//    Per=PerIx*RR/NN;
//    for(s0Ix=1;s0Ix<=NN;s0Ix++) {
//      for(s1Ix=s0Ix+1;s1Ix<NN;s1Ix++) {
//	//for(s1Ix=s0Ix;s1Ix<=NN;s1Ix++) {
//	s2Ix=NN;
//	if(s2Ix>=s0Ix+s1Ix)continue;
//	SumSideIx=s0Ix+s1Ix+s2Ix;
//	Tri_p->Side[0]=s0Ix*Per/SumSideIx;
//	Tri_p->Side[1]=s1Ix*Per/SumSideIx;
//	Tri_p->Side[2]=s2Ix*Per/SumSideIx;
//	for(j=0;j<3;j++) {
//	  Cosh[j]=COSH(Tri_p->Side[j]);
//	  SinhTmp=SINH(Tri_p->Side[j]);
//	  Sinh2[j]=SinhTmp*SinhTmp;
//	}
//	ShapeMax=Cosh[0]*Cosh[1]*Cosh[2];
//	ShapeMin=0.5*(Sinh2[0]+Sinh2[1]+Sinh2[2]+2.0);
//	for(ShapeIx=1;ShapeIx<NN;ShapeIx++) {
//	  printf("%4i;%4i%4i%4i;%4i: ",PerIx,s0Ix,s1Ix,s2Ix,ShapeIx);
//	  Tri_p->Shape=ShapeMin+ShapeIx*(ShapeMax-ShapeMin)/NN;
//	  printf("%12g %12g %12g %12g ",
//		 Tri_p->Side[0],
//		 Tri_p->Side[1],
//		 Tri_p->Side[2],
//		 Tri_p->Shape);
//	  CHypTriFromSidesAndShape(Tri_p);
//	  CHypTriResultant(Tri_p);
//	  CHypTriFindU(uSol,Tri_p);
//	  CHypTriSinhR2(SinhR2,&RealCnt,uSol,Tri_p);
//	  assert(RealCnt == 4 || RealCnt == 6);
//	  printf("%i",RealCnt);
//	  PosCnt=0;
//	  for(Ix=0;Ix<RealCnt;Ix++)
//	    if(SinhR2[Ix]>0) PosCnt++;
//	  assert(PosCnt<=1);
//	  if(PosCnt==0)
//	    printf("N\n");
//	  else
//	    printf("O\n");
//
//	  {
//	    for(Ix=0;Ix<RealCnt;Ix++)
//	      if(SinhR2[Ix]>0) {
//		printf("  R=%15.8g \\sinh^2 R=%15.8g\n",
//		       ASINH(SQRT(SinhR2[Ix])),SinhR2[Ix]);
//	    }
//	  }
//	}
//      }
//    }
//  }
//}

// Calculates the triangle invariants on the basis of the three side
// lengths and Brehm's shape.
void CHypTriFromSidesAndShape(CHypTri_t *Tri_p) {
  FLOAT_t ProdCosh;

  FLOAT_t cosh0,cosh1,cosh2; // Hyperbolic cosines of the side lengths
  FLOAT_t s0,s1,s2; // Symmetric functions in the sinh^2's
  FLOAT_t c0,c1,c2; // Symmetric functions in the cosh^2's
  // Shape, the real part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hL1=Tri_p->ShapeL1;
  const FLOAT_t h=1.0+hL1;
  FLOAT_t hi; // The imaginary part of $(A,B)(B,C)(C,A)$
  FLOAT_t d; // The determinant of the matrix of inner products
  FLOAT_t v0,v1,v2; // Symmetric functions in the squared absolute
		    // values of the off-diagonal entries of the
		    // inverted matrix of inner products
  FLOAT_t vt; // The trace of the inverted matrix of inner products
  FLOAT_t vh; // Real part of the cyclic inner product of the
	      // off-diagonal entries of the inverted matrix of inner
	      // products --- CREAL(AInv[0][1]*AInv[1][2]*AInv [2][0])
  FLOAT_t vhi; // Imaginary part of the cyclic inner product of the
	       // off-diagonal entries of the inverted matrix of inner
	       // products
 
  assert(Tri_p->Side[0]<=Tri_p->Side[1]+Tri_p->Side[2]);
  assert(Tri_p->Side[1]<=Tri_p->Side[2]+Tri_p->Side[0]);
  assert(Tri_p->Side[2]<=Tri_p->Side[0]+Tri_p->Side[1]);

  cosh0=COSH(Tri_p->Side[0]);
  cosh1=COSH(Tri_p->Side[1]);
  cosh2=COSH(Tri_p->Side[2]);

  // Calculate the imaginary part of (A,B)(B,C)(C,A)
  ProdCosh=cosh0*cosh1*cosh2;
  assert(h<=ProdCosh);
  hi=SQRT(ProdCosh*ProdCosh-h*h);


  // Calculate the symmetric polynomials in the squared cosh's
  {
    const FLOAT_t chSq0=cosh0*cosh0,chSq1=cosh1*cosh1,chSq2=cosh2*cosh2;
    c0=chSq0+chSq1+chSq2;
    c1=chSq0*chSq1+chSq1*chSq2+chSq2*chSq0;
    c2=chSq0*chSq1*chSq2;
  }

  // Calculate the symmetric polynomials in the squared sinh's
  {
    const FLOAT_t sinh0=SINH(Tri_p->Side[0]);
    const FLOAT_t sinh1=SINH(Tri_p->Side[1]);
    const FLOAT_t sinh2=SINH(Tri_p->Side[2]);
    const FLOAT_t shSq0=sinh0*sinh0,shSq1=sinh1*sinh1,shSq2=sinh2*sinh2;
    s0=shSq0+shSq1+shSq2;
    s1=shSq0*shSq1+shSq1*shSq2+shSq2*shSq0;
    s2=shSq0*shSq1*shSq2;
  }

  // Calculate the other invariants
  d=(2.0*(h-1.0)-s0);
  vt=(-s0)/d;
  vh=(h*(d-2.0) + s0 - s2 + 2.0) / (d*d*d);
  vhi= (-hi)/(d*d);
  v0=( - 6.0*(h-1.0) + 3.0*s0 + s1) / (d*d);
  v1=(4.0*(3.0*h - (3.0*s0 + s1 + 6.0))*h + (3.0*s0 + 2.0*s1 + s2 + 12.0)*s0 
      + 4.0*s1 + 12.0) / (d*d*d*d);
  v2=( - 2.0*(2.0*h*(2.0*h - (3.0*s0 + s1 + 6.0))
	      + s0*(3.0*s0 + 2.0*s1 + s2 + 12.0) + 4*(s1 + 3.0))*h 
       + (s0*(s0 + s1 + s2 + 6.0) + 2.0*(2.0*s1 + s2 + 6.0))*s0
       + 4.0*s1 + s2*s2 + 8.0) / (d*d*d*d*d*d);

  // Load *Tri_p with the calculated values
  (Tri_p->ShapeIm)=hi;
  (Tri_p->Det)=d;
  (Tri_p->SymSinh2)[0]=s0; (Tri_p->SymSinh2)[1]=s1; (Tri_p->SymSinh2)[2]=s2;
  (Tri_p->SymCosh2)[0]=c0; (Tri_p->SymCosh2)[1]=c1; (Tri_p->SymCosh2)[2]=c2;
  (Tri_p->SymInv)[0]=v0; (Tri_p->SymInv)[1]=v1; (Tri_p->SymInv)[2]=v2;
  (Tri_p->TrInv)=vt;
  (Tri_p->InvShape)=vh;
  (Tri_p->InvShapeIm)=vhi;

  // Set the flags
  (Tri_p->Flags)=0;
  if(FABS(hi)<epsBigger*FABS(h))
    (Tri_p->Flags)|=REALFG;
  if(FABS(Tri_p->Side[0]-Tri_p->Side[1])/Tri_p->Side[1] < epsBigger
     || FABS(Tri_p->Side[1]-Tri_p->Side[2])/Tri_p->Side[2] < epsBigger
     || FABS(Tri_p->Side[2]-Tri_p->Side[0])/Tri_p->Side[0] < epsBigger) {
    (Tri_p->Flags)|=ISOFG;
  }
  if(d<epsBig*hL1)
    (Tri_p->Flags)|=CPXFG;
  if(d<epsHuge*hL1)
    (Tri_p->Flags)|=NEARCPXFG;

  CheckCHypTri(Tri_p);
}

// Calculates the triangle invariants on the basis of three given
// vertices in the unit ball$B(\CC^2)$
void CHypTriFromPts(CHypTri_t *Tri_p,
		    const Pt_t PtA,const Pt_t PtB,const Pt_t PtC) {

  FLOAT_t shSq0,shSq1,shSq2; // sinh^2's of the side lengths

  FLOAT_t s0,s1,s2; // Symmetric functions in the sinh^2's
  FLOAT_t c0,c1,c2; // Symmetric functions in the cosh^2's
  // Shape, the real part of $(A,B)(B,C)(C,A)$
  FLOAT_t h;
  FLOAT_t hL1; // h, LESS 1
  FLOAT_t hi; // The imaginary part of $(A,B)(B,C)(C,A)$
  FLOAT_t d; // The determinant of the matrix of inner products
  FLOAT_t v0,v1,v2; // Symmetric functions in the squared absolute
		    // values of the off-diagonal entries of the
		    // inverted matrix of inner products
  FLOAT_t vt; // The trace of the inverted matrix of inner products
  FLOAT_t vh; // Real part of the cyclic inner product of the
	      // off-diagonal entries of the inverted matrix of inner
	      // products --- CREAL(AInv[0][1]*AInv[1][2]*AInv [2][0])
  FLOAT_t vhi; // Imaginary part of the cyclic inner product of the
	       // off-diagonal entries of the inverted matrix of inner
	       // products

  const FLOAT_t PtAPtA=PtAbsSq(PtA);
  const FLOAT_t PtBPtB=PtAbsSq(PtB);
  const FLOAT_t PtCPtC=PtAbsSq(PtC);

  const COMPLEX_t PtAPtB=Scal(PtA,PtB);
  const COMPLEX_t PtBPtC=Scal(PtB,PtC);
  const COMPLEX_t PtCPtA=Scal(PtC,PtA);

  const FLOAT_t ANorm2L1=PtAPtA/(1.0-PtAPtA);
  const FLOAT_t BNorm2L1=PtBPtB/(1.0-PtBPtB);
  const FLOAT_t CNorm2L1=PtCPtC/(1.0-PtCPtC);

  const FLOAT_t ANorm2=1.0+ANorm2L1;
  const FLOAT_t BNorm2=1.0+BNorm2L1;
  const FLOAT_t CNorm2=1.0+CNorm2L1;

  // The three vectors in $\CC^3$ are:
  //
  //                    ( PtA )
  //   A = SQRT(ANorm2) (     )  , B = ..., C = ...
  //                    (  1  )
  //
  // This is (A,B)(B,C)(C,A)
  const COMPLEX_t hCpx=
    ((1.0-PtAPtB)*(1.0-PtBPtC)*(1.0-PtCPtA))*ANorm2*BNorm2*CNorm2;
  // This is the same, LESS~1, calculated to avoid round off
  // error when PtA, PtB, and PtC are small
  const COMPLEX_t hCpxL1=
    -PtAPtB*(1.0-PtBPtC)*(1.0-PtCPtA)*ANorm2*BNorm2*CNorm2
    -PtBPtC*(1.0-PtCPtA)*ANorm2*BNorm2*CNorm2
    -PtCPtA*ANorm2*BNorm2*CNorm2
    +ANorm2L1*BNorm2*CNorm2
    +BNorm2L1*CNorm2
    +CNorm2L1;

  // Real and complex parts of (A,B)(B,C)(C,A)
  hL1=CREAL(hCpxL1);
  h=1.0+hL1;
  hi=CIMAG(hCpxL1);

  // Squared sinh's of the side lengths.  The first one, for instance, is
  // given as:
  //
  //   |PtB - PtC|^2 / ((1 - |PtB|^2) (1-|PtC|^2))

  {
    Pt_t Diff; 
    FLOAT_t AbsDet;
    SubPts(Diff,PtB,PtC);
    AbsDet=CABS(PtB[0]*PtC[1]-PtB[1]*PtC[0]);
    shSq0=(PtAbsSq(Diff)-AbsDet*AbsDet)*BNorm2*CNorm2;
    SubPts(Diff,PtC,PtA);
    AbsDet=CABS(PtC[0]*PtA[1]-PtC[1]*PtA[0]);
    shSq1=(PtAbsSq(Diff)-AbsDet*AbsDet)*CNorm2*ANorm2;
    SubPts(Diff,PtA,PtB);
    AbsDet=CABS(PtA[0]*PtB[1]-PtA[1]*PtB[0]);
    shSq2=(PtAbsSq(Diff)-AbsDet*AbsDet)*ANorm2*BNorm2;
  }

  // The side lengths themselves
  (Tri_p->Side)[0]=ASINH(SQRT(shSq0));
  (Tri_p->Side)[1]=ASINH(SQRT(shSq1));
  (Tri_p->Side)[2]=ASINH(SQRT(shSq2));

  // The symmetric polynomials in the squared sinh's
  {
    s0=shSq0+shSq1+shSq2;
    s1=shSq0*shSq1+shSq1*shSq2+shSq2*shSq0;
    s2=shSq0*shSq1*shSq2;
  }

  // The symmetric polynomials in the squared cosh's
  {
    const FLOAT_t chSq0=shSq0+1.0,chSq1=shSq1+1.0,chSq2=shSq2+1.0;
    c0=chSq0+chSq1+chSq2;
    c1=chSq0*chSq1+chSq1*chSq2+chSq2*chSq0;
    c2=chSq0*chSq1*chSq2;
  }

  // Calculate the other invariants
  d=(2.0*hL1-s0);
  vt=(-s0)/d;
  vh=((2*hL1 - s0)*hL1 - s2) / (d*d*d);
  vhi=(-hi)/(d*d);
  v0=(- 6.0*hL1 + 3.0*s0 + s1) / (d*d);
  v1=(- 4*(3*s0 + s1 - 3*hL1)*hL1 + (2*s1 + s2 + 3*s0)*s0) / (d*d*d*d);
  v2=(2*(2*(3*s0 + s1 - 2*hL1)*hL1 - (2*s1 + s2 + 3*s0)*s0)*hL1
      + (s1 + s2 + s0)*s0*s0 + s2*s2) / (d*d*d*d*d*d);

  // Load *Tri_p with the calculated values
  (Tri_p->ShapeL1)=hL1;
  (Tri_p->ShapeIm)=hi;
  (Tri_p->Det)=d;
  (Tri_p->SymSinh2)[0]=s0; (Tri_p->SymSinh2)[1]=s1; (Tri_p->SymSinh2)[2]=s2;
  (Tri_p->SymCosh2)[0]=c0; (Tri_p->SymCosh2)[1]=c1; (Tri_p->SymCosh2)[2]=c2;
  (Tri_p->SymInv)[0]=v0; (Tri_p->SymInv)[1]=v1; (Tri_p->SymInv)[2]=v2;
  (Tri_p->TrInv)=vt;
  (Tri_p->InvShape)=vh;
  (Tri_p->InvShapeIm)=vhi;

  // Set the flags
  (Tri_p->Flags)=0;
  if(FABS(hi)<epsBigger*FABS(h))
    (Tri_p->Flags)|=REALFG;
  if(FABS(Tri_p->Side[0]-Tri_p->Side[1])/Tri_p->Side[1] < epsBig
     || FABS(Tri_p->Side[1]-Tri_p->Side[2])/Tri_p->Side[2] < epsBig
     || FABS(Tri_p->Side[2]-Tri_p->Side[0])/Tri_p->Side[0] < epsBig) {
    (Tri_p->Flags)|=ISOFG;
  }
  if(d<epsBig*hL1)
    (Tri_p->Flags)|=CPXFG;
  if(d<epsHuge*hL1)
    (Tri_p->Flags)|=NEARCPXFG;

  CheckCHypTri(Tri_p);
}

void CheckCHypTri(const CHypTri_t *Tri_p) {
  Mtx_t A,AInv; // The matrix of inner products and its inverse
  COMPLEX_t omega;
  int RIx,CIx;

  // Hyperbolic cosines of the side lengths
  const FLOAT_t cosh0=COSH((Tri_p->Side)[0]);
  const FLOAT_t cosh1=COSH((Tri_p->Side)[1]);
  const FLOAT_t cosh2=COSH((Tri_p->Side)[2]);

  const FLOAT_t hL1=Tri_p->ShapeL1; // Shape LESS 1, the real part of
				    // $(A,B)(B,C)(C,A) - 1$
  const FLOAT_t h=hL1+1.0;
  const FLOAT_t hi=Tri_p->ShapeIm; // The imaginary part of
				   // $(A,B)(B,C)(C,A)$
  FLOAT_t d=Tri_p->Det; // The determinant of the matrix of inner products
  // Symmetric functions in the sinh^2's
  const FLOAT_t s0=(Tri_p->SymSinh2)[0];
  const FLOAT_t s1=(Tri_p->SymSinh2)[1];
  const FLOAT_t s2=(Tri_p->SymSinh2)[2];
  // Symmetric functions in the cosh^2's
  const FLOAT_t c0=(Tri_p->SymCosh2)[0];
  const FLOAT_t c1=(Tri_p->SymCosh2)[1];
  const FLOAT_t c2=(Tri_p->SymCosh2)[2];
  // Symmetric functions in the squared absolute values of the
  // off-diagonal entries of the inverted matrix of inner products.
  const FLOAT_t v0=(Tri_p->SymInv)[0];
  const FLOAT_t v1=(Tri_p->SymInv)[1];
  const FLOAT_t v2=(Tri_p->SymInv)[2];
  // The trace of the inverted matrix of inner products
  const FLOAT_t vt=(Tri_p->TrInv);
  // Real part of the cyclic inner product of the off-diagonal entries
  // of the inverted matrix of inner products ---
  // CREAL(AInv[0][1]*AInv[1][2]*AInv[2][0])
  const FLOAT_t vh=(Tri_p->InvShape);
  // Imaginary part of the cyclic inner product of the off-diagonal
  // entries of the inverted matrix of inner products
  const FLOAT_t vhi=(Tri_p->InvShapeIm);

  // Checks --- the triangle inequality
  assert(Tri_p->Side[0]<=Tri_p->Side[1]+Tri_p->Side[2]+epsSmall);
  assert(Tri_p->Side[1]<=Tri_p->Side[2]+Tri_p->Side[0]+epsSmall);
  assert(Tri_p->Side[2]<=Tri_p->Side[0]+Tri_p->Side[1]+epsSmall);

  // Checks --- Isosceles triangle
  if(Tri_p->Flags & ISOFG) {
    FLOAT_t Err;
    Err=MIN(FABS((Tri_p->Side[0]-Tri_p->Side[1])/Tri_p->Side[1]),
	    FABS((Tri_p->Side[1]-Tri_p->Side[2])/Tri_p->Side[2]));
    Err=MIN(Err,
	    FABS((Tri_p->Side[2]-Tri_p->Side[0])/Tri_p->Side[0]));	    ;
    printf("Isosceles, side lengths: %12.6g %12.6g %12.6g"
	   "  Relative error= %14.8g\n",
	   Tri_p->Side[0],
	   Tri_p->Side[1],
	   Tri_p->Side[2],
	   Err);
  }

  // Checks --- the symmetric functions in the sinh's of the sides
  {
    const FLOAT_t sinh0=SINH((Tri_p->Side)[0]);
    const FLOAT_t sinh1=SINH((Tri_p->Side)[1]);
    const FLOAT_t sinh2=SINH((Tri_p->Side)[2]);
    const FLOAT_t sinh0Sq=sinh0*sinh0;
    const FLOAT_t sinh1Sq=sinh1*sinh1;
    const FLOAT_t sinh2Sq=sinh2*sinh2;
    const FLOAT_t s0Tmp=sinh0Sq+sinh1Sq+sinh2Sq;
    const FLOAT_t s1Tmp=sinh0Sq*sinh1Sq+sinh1Sq*sinh2Sq+sinh2Sq*sinh0Sq;
    const FLOAT_t s2Tmp=sinh0Sq*sinh1Sq*sinh2Sq;
    assert(FABS(s0-s0Tmp)<epsBig*s0);
    assert(FABS(s1-s1Tmp)<epsBig*s1);
    assert(FABS(s2-s2Tmp)<epsBig*s2);
  }

  // Checks --- the symmetric functions in the cosh's of the sides
  {
    const FLOAT_t cosh0Sq=cosh0*cosh0;
    const FLOAT_t cosh1Sq=cosh1*cosh1;
    const FLOAT_t cosh2Sq=cosh2*cosh2;
    const FLOAT_t c0Tmp=cosh0Sq+cosh1Sq+cosh2Sq;
    const FLOAT_t c1Tmp=cosh0Sq*cosh1Sq+cosh1Sq*cosh2Sq+cosh2Sq*cosh0Sq;
    const FLOAT_t c2Tmp=cosh0Sq*cosh1Sq*cosh2Sq;
    assert(FABS(c0-c0Tmp)<epsBig*c0);
    assert(FABS(c1-c1Tmp)<epsBig*c1);
    assert(FABS(c2-c2Tmp)<epsBig*c2);
  }

  // Checks --- h and hi
  {
    assert(FABS(h*h+hi*hi-c2)<epsBig*c2);
    if(Tri_p->Flags & REALFG) {
      printf("Real triangle: ");
      printf("hi=%16.8g h=%16.8g hi/h=%16.8g\n",hi,h,hi/h);
    }
  }

  if(Tri_p->Flags & CPXFG) {
    printf("Complex triangle: ");
    printf("d=%16.8g hL1=%16.8g d/hL1=%16.8g\n",d,hL1,d/hL1);
    return;
  }

  // Calculate (one version) of the matrix of inner products
  omega=(h+I*hi) / (cosh0*cosh1*cosh2);
  assert(FABS(CABS(omega)-1.0) <= epsSmall);

  A[0][0]=A[1][1]=A[2][2]=1.0;
  A[0][1]=A[1][0]=cosh2;
  A[2][0]=A[0][2]=cosh1;
  A[1][2]=cosh0*omega;  A[2][1]=CONJ(A[1][2]);

  // Calculate its inverse
  InvMtx(AInv,A);

  // Checks --- the matrix AInv
  {
    for(RIx=0;RIx<3;RIx++) 
      for(CIx=0;CIx<3;CIx++)
	assert(CABS(AInv[RIx][CIx]-CONJ(AInv[CIx][RIx]))
	       <epsBigger*CABS(AInv[RIx][CIx]));
  }

  // Checks ---- the determinant of the matrix
  {
    const COMPLEX_t DetTmp=Det3By3(A);
    //printf("d=%.8g, DetTmp=%.8g%+.8g*i\n",d,DetTmp);
    assert(d>epsSmall); // This excludes the case of a triangle lying in a
                        // complex hyperbolic disk
    assert(CABS(d-DetTmp)<epsBigger*d);
  }
  {
    const FLOAT_t DetTmp=(2.0*(h-1.0)-s0);
    assert(CABS(d-DetTmp)<epsBig*d);
  }
  {
    const FLOAT_t DetTmp=2*h - c0 + 1.0;
    assert(CABS(d-DetTmp)<epsHuge*d);
  }

  // Checks --- the product A[0][1]*A[1][2]+A[2][0]
  {
    const COMPLEX_t hCpx=A[0][1]*A[1][2]*A[2][0];
    const FLOAT_t hTmp=CREAL(hCpx);
    const FLOAT_t hiTmp=CIMAG(hCpx);
    assert(FABS(h-hTmp) < epsBig*FABS(h));
    if(FABS(hi) > epsSmall)
      assert(FABS(hi-hiTmp) < epsBigger*FABS(hi));
  }

  // Checks --- the trace of the inverse matrix
  {
    const COMPLEX_t TrTmp=AInv[0][0]+AInv[1][1]+AInv[2][2];
    //printf("vt=%.8g TrTmp=%.8g\n",vt,TrTmp);
    assert(CABS(vt-TrTmp) < epsBigger*FABS(TrTmp));
  }
  {
    const FLOAT_t TrTmp=(-s0)/d;
    assert(FABS(vt-TrTmp) < epsBig*FABS(TrTmp));
  }
  {
    const FLOAT_t TrTmp=(3.0-c0)/d;
    assert(FABS(vt-TrTmp) < epsBigger*FABS(TrTmp));
  }

  // Checks --- the product AInv[0][1]*AInv[1][2]+AInv[2][0]
  {
    const COMPLEX_t vhCpx=AInv[0][1]*AInv[1][2]*AInv[2][0];
    const FLOAT_t vhTmp=CREAL(vhCpx);
    //const FLOAT_t vhiTmp=CIMAG(vhCpx);
    const FLOAT_t vhiTmp=
      +CREAL(AInv[0][1])*(+CREAL(AInv[1][2])*CIMAG(AInv[2][0])
			  +CIMAG(AInv[1][2])*CREAL(AInv[2][0]))
      +CIMAG(AInv[0][1])*(+CREAL(AInv[1][2])*CREAL(AInv[2][0])
			  -CIMAG(AInv[1][2])*CIMAG(AInv[2][0]));
    assert(FABS(vh-vhTmp) < epsHuge*FABS(vh));
    if(!(Tri_p->Flags & REALFG))
      if(FABS(vhi-vhiTmp)>epsBigger*FABS(vhi)) {
	printf("h=%.8g hi=%.8g d=%.8g hL1=%.8g d/hL1=%.8g\n",
	       h,hi,d,hL1,d/hL1);
	printf("vhi=%15.8g vhiTmp=%16.8g Relative Error=%15.8g\n",
	       vhi,vhiTmp,(vhi-vhiTmp)/vhi);
	if(!(Tri_p->Flags & NEARCPXFG))assert(FALSE);
      }
  }
  {
    const FLOAT_t vhTmp=(h*(d-2.0) + s0 - s2 + 2.0) / (d*d*d);
    assert(FABS(vh-vhTmp) < epsBigger*FABS(vh));
  }
  {
    const FLOAT_t vhTmp=(- c2 + h*(d - 2.0) + c1) / (d*d*d);
    assert(FABS(vh-vhTmp) < epsHuge*FABS(vh));
  }

  // Checks --- the symmetric polynomials in the squared absolute values
  // of the off-diagonal elements of AInv
  {
    const FLOAT_t abs01=CABS(AInv[0][1]);
    const FLOAT_t abs12=CABS(AInv[1][2]);
    const FLOAT_t abs20=CABS(AInv[2][0]);
    const FLOAT_t absSq01=abs01*abs01;
    const FLOAT_t absSq12=abs12*abs12;
    const FLOAT_t absSq20=abs20*abs20;
    const FLOAT_t v0tmp=absSq01+absSq12+absSq20;
    const FLOAT_t v1tmp=absSq01*absSq12+absSq12*absSq20+absSq20*absSq01;
    const FLOAT_t v2tmp=absSq01*absSq12*absSq20;
    //printf("v0=%.8g v0Tmp=%.8g\n",v0,v0tmp);
    assert(FABS(v0-v0tmp)<epsHuge*FABS(v0));
    //printf("v1=%.8g v1Tmp=%.8g\n",v1,v1tmp);
    assert(FABS(v1-v1tmp)<epsHuge*FABS(v1));
    //printf("v2=%.8g v2Tmp=%.8g\n",v2,v2tmp);
    assert(FABS(v2-v2tmp)<epsHuge*FABS(v2));
  }
  {
    const FLOAT_t v0tmp=( - 6.0*(h-1.0) + 3.0*s0 + s1)/(d*d);
    const FLOAT_t v1tmp=
      (4.0*(3.0*h - (3.0*s0 + s1 + 6.0))*h
       + (3.0*s0 + 2.0*s1 + s2 + 12.0)*s0 + 4.0*s1 + 12.0) / (d*d*d*d);
    const FLOAT_t v2tmp=
      ( - 2.0*(2.0*h*(2.0*h - (3.0*s0 + s1 + 6.0))
	       + s0*(3.0*s0 + 2.0*s1 + s2 + 12.0) + 4*(s1 + 3.0))*h 
	+ (s0*(s0 + s1 + s2 + 6.0) + 2.0*(2.0*s1 + s2 + 6.0))*s0
	+ 4.0*s1 + s2*s2 + 8.0) / (d*d*d*d*d*d);
    assert(FABS(v0-v0tmp)<epsBig*FABS(v0));
    assert(FABS(v1-v1tmp)<epsBig*FABS(v1));
    assert(FABS(v2-v2tmp)<epsBigger*FABS(v2));
  }
  {
    const FLOAT_t v0tmp=(- 6.0*h + c0 + c1)/(d*d);
    const FLOAT_t v1tmp=
      (4.0*h*(3.0*h - (c0 + c1)) + c0*(c1 + c2) + c1 - 3.0*c2) / (d*d*d*d);
    const FLOAT_t v2tmp=
      (- 2.0*h*(2.0*h*(2.0*h - (c0 + c1)) + c0*(c1 + c2) + c1 - 3.0*c2)
       + c0*c2*(c0 - 2.0) + c1*(c1 - 2*c2) + c2*(c2 + 1.0)) / (d*d*d*d*d*d);
    assert(FABS(v0-v0tmp)<epsBig*FABS(v0));
    assert(FABS(v1-v1tmp)<epsBigger*FABS(v1));
    //printf("v2=%20.15g, v2tmp=%20.15g\n",v2,v2tmp);
    assert(FABS(v2-v2tmp)<epsHuge*FABS(v2));
  }
}

void CHypTriFindU(FLOAT_t *uSol,const CHypTri_t *Tri_p) {
  // Shape, the real part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hL1=Tri_p->ShapeL1;
  const FLOAT_t h=1.0+hL1;
  // The imaginary part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hi=Tri_p->ShapeIm; 
  // The determinant of the matrix of inner products
  FLOAT_t d=Tri_p->Det;
  // Symmetric functions in the sinh^2's
  const FLOAT_t s0=(Tri_p->SymSinh2)[0];
  const FLOAT_t s1=(Tri_p->SymSinh2)[1];
  const FLOAT_t s2=(Tri_p->SymSinh2)[2];
  // Symmetric functions in the cosh^2's
  const FLOAT_t c0=(Tri_p->SymCosh2)[0];
  const FLOAT_t c1=(Tri_p->SymCosh2)[1];
  const FLOAT_t c2=(Tri_p->SymCosh2)[2];
  // Symmetric functions in the squared absolute values of the
  // off-diagonal entries of the inverted matrix of inner products.
  const FLOAT_t v0=(Tri_p->SymInv)[0];
  const FLOAT_t v1=(Tri_p->SymInv)[1];
  const FLOAT_t v2=(Tri_p->SymInv)[2];
  // The trace of the inverted matrix of inner products
  const FLOAT_t vt=(Tri_p->TrInv);
  // Real part of the cyclic inner product of the off-diagonal entries
  // of the inverted matrix of inner products ---
  // CREAL(AInv[0][1]*AInv[1][2]*AInv[2][0])
  const FLOAT_t vh=(Tri_p->InvShape);
  // Imaginary part of the cyclic inner product of the off-diagonal
  // entries of the inverted matrix of inner products
  const FLOAT_t vhi=(Tri_p->InvShapeIm);

  const FLOAT_t vhSq=vh*vh,v1Sq=v1*v1,v2Sq=v2*v2;

  FLOAT_t uEq[7]=
    {16*(vhSq*(vhSq - 2*v2) + v2Sq),
     0,
     8*v1*(vhSq - v2),
     8*vhi*(vhSq + v2),
     - 4*vhSq*v0 + v1Sq,
     - 2*vhi*v1,
     v2};

  gsl_poly_complex_workspace *ws
    = gsl_poly_complex_workspace_alloc(7);

  gsl_poly_complex_solve(uEq,7,ws,uSol);

  gsl_poly_complex_workspace_free(ws);
}

// Calculates and displays the signs of the three factors of the
// resultant of the equation for~$u$.  The first two are always
// positive and negative, respectively, with the proviso that
// everything goes kablooey when the triangle is isoceles (and
// probably also if it is contained in a complex disc.)  The last sign
// comes out:
//   
//   - when there are 4~real solutions for~$u$
//   + when there are 6~real solutions for~$u$

void CHypTriResultant(const CHypTri_t *Tri_p) {
  // Shape, the real part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hL1=Tri_p->ShapeL1;
  const FLOAT_t h=hL1;
  // The imaginary part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hi=Tri_p->ShapeIm; 
  // The determinant of the matrix of inner products
  FLOAT_t d=Tri_p->Det;
  // Symmetric functions in the sinh^2's
  const FLOAT_t s0=(Tri_p->SymSinh2)[0];
  const FLOAT_t s1=(Tri_p->SymSinh2)[1];
  const FLOAT_t s2=(Tri_p->SymSinh2)[2];
  // Symmetric functions in the cosh^2's
  const FLOAT_t c0=(Tri_p->SymCosh2)[0];
  const FLOAT_t c1=(Tri_p->SymCosh2)[1];
  const FLOAT_t c2=(Tri_p->SymCosh2)[2];
  // Symmetric functions in the squared absolute values of the
  // off-diagonal entries of the inverted matrix of inner products.
  const FLOAT_t v0=(Tri_p->SymInv)[0];
  const FLOAT_t v1=(Tri_p->SymInv)[1];
  const FLOAT_t v2=(Tri_p->SymInv)[2];
  // The trace of the inverted matrix of inner products
  const FLOAT_t vt=(Tri_p->TrInv);
  // Real part of the cyclic inner product of the off-diagonal entries
  // of the inverted matrix of inner products ---
  // CREAL(AInv[0][1]*AInv[1][2]*AInv[2][0])
  const FLOAT_t vh=(Tri_p->InvShape);
  // Imaginary part of the cyclic inner product of the off-diagonal
  // entries of the inverted matrix of inner products
  const FLOAT_t vhi=(Tri_p->InvShapeIm);

  const FLOAT_t vhSq=vh*vh,v1Sq=v1*v1,v2Sq=v2*v2,s1Sq=s1*s1;
  const FLOAT_t ResFactors[3]=
    {
      v2,

      s0*(s0*(4*s0*s2 - s1Sq) - 18*s1*s2) + 4*s1*s1Sq + 27*s2*s2,

      4*v0*(4*v0*(4*v0*vhSq - 3*v1Sq)*vhSq
	    + 3*v1*(v1Sq*v1 + 18*(v2*(3*v2 - 2*vhSq)  - vhSq*vhSq)))*vhSq
      - v1Sq*v1*(v1Sq*v1 - 54*(v2*(v2 - 6*vhSq) + 5*vhSq*vhSq))
      - 27*v2*(9*v2*(3*v2Sq - 2*vhSq*vhSq) - 8*vhSq*vhSq*vhSq)
      + 27*vhSq*vhSq*vhSq*vhSq
    };
  int Ix;

  for(Ix=0;Ix<3;Ix++) {
    //if(ResFactors[Ix]<=0)
    //  printf("-");
    //else 
    //  printf("+");
  }
  //printf("\n");
}


// Let~$R$ be the circumradius of a complex hyperbolic triangle.
// Given the solutions for~$u$, this routine calculates the possible
// real values of $\sinh^2 R$.

void CHypTriSinhR2(FLOAT_t *SinhR2,int *RealCnt_p,
		   const FLOAT_t *uSol,const CHypTri_t *Tri_p) {
  // Shape, the real part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hL1=Tri_p->ShapeL1;
  const FLOAT_t h=1.0+hL1;
  // The imaginary part of $(A,B)(B,C)(C,A)$
  const FLOAT_t hi=Tri_p->ShapeIm; 
  // The determinant of the matrix of inner products
  FLOAT_t d=Tri_p->Det;
  // Symmetric functions in the sinh^2's
  const FLOAT_t s0=(Tri_p->SymSinh2)[0];
  const FLOAT_t s1=(Tri_p->SymSinh2)[1];
  const FLOAT_t s2=(Tri_p->SymSinh2)[2];
  // Symmetric functions in the cosh^2's
  const FLOAT_t c0=(Tri_p->SymCosh2)[0];
  const FLOAT_t c1=(Tri_p->SymCosh2)[1];
  const FLOAT_t c2=(Tri_p->SymCosh2)[2];
  // Symmetric functions in the squared absolute values of the
  // off-diagonal entries of the inverted matrix of inner products.
  const FLOAT_t v0=(Tri_p->SymInv)[0];
  const FLOAT_t v1=(Tri_p->SymInv)[1];
  const FLOAT_t v2=(Tri_p->SymInv)[2];
  // The trace of the inverted matrix of inner products
  const FLOAT_t vt=(Tri_p->TrInv);
  // Real part of the cyclic inner product of the off-diagonal entries
  // of the inverted matrix of inner products ---
  // CREAL(AInv[0][1]*AInv[1][2]*AInv[2][0])
  const FLOAT_t vh=(Tri_p->InvShape);
  // Imaginary part of the cyclic inner product of the off-diagonal
  // entries of the inverted matrix of inner products
  const FLOAT_t vhi=(Tri_p->InvShapeIm);

  const FLOAT_t vhSq=vh*vh,v1Sq=v1*v1;

  int Ix;
  const FLOAT_t *uR,*uI;

  uR=&(uSol[0]);
  uI=&(uSol[1]);
  (*RealCnt_p)=0;
  for(Ix=0;Ix<6;Ix++) {
    // Rather arbitrary cutoff
    if(FABS(*uI) < 1E-7*FABS(*uR)) {
      const FLOAT_t u=(*uR);
      SinhR2[*RealCnt_p]=
	( - u*(u*(u*(u*v2 - 2*v1*vhi) - (4*vhSq*v0 - v1Sq))
	       + 4*(3*vhSq + v2)*vhi)
	  - 4*vh*(vh*(vh*(vt - 1) + v1) - v2*(vt - 1)) + 4*v1*v2) / 
	(u*(u*(u*(u*v2 - 2*v1*vhi) - (4*vhSq*v0 - v1Sq))
	    + 4*(3*vhSq + v2)*vhi) + 4*vh*(vh*(vh*vt + v1) - v2*vt)
	 - 4*v1*v2);
      if(SinhR2[*RealCnt_p]>0) {
	//printf("\n1/(cvalvh+vt)-1-sinhR2 "
	//       " where vt=>%.10g,\n  "
	//       "v0=>%.10g, "
	//       "v1=>%.10g, "
	//       "v2=>%.10g,\n  "
	//       "vh=>%.10g, "
	//       "vhi=>%.10g, "
	//       "uuvh=>%.10g,\n  "
	//       "sinhR2=>%.10g;\n",
	//       vt,v0,v1,v2,vh,vhi,u,SinhR2[*RealCnt_p]);
      }
      (*RealCnt_p)++;
    }
    uR+=2; uI+=2;
  }
}

// Calculates the circumradius of three points in~$B(\CC^2)$
FLOAT_t CircumRadius(Pt_t X,Pt_t PtA,Pt_t PtB,Pt_t PtC) {

  const FLOAT_t PtAPtA=PtAbsSq(PtA);
  const FLOAT_t PtBPtB=PtAbsSq(PtB);
  const FLOAT_t PtCPtC=PtAbsSq(PtC);

  const COMPLEX_t PtAPtB=Scal(PtA,PtB);
  const COMPLEX_t PtBPtC=Scal(PtB,PtC);
  const COMPLEX_t PtCPtA=Scal(PtC,PtA);

  const FLOAT_t ANorm2L1=PtAPtA/(1.0-PtAPtA);
  const FLOAT_t BNorm2L1=PtBPtB/(1.0-PtBPtB);
  const FLOAT_t CNorm2L1=PtCPtC/(1.0-PtCPtC);

  const FLOAT_t ANorm2=1.0+ANorm2L1;
  const FLOAT_t BNorm2=1.0+BNorm2L1;
  const FLOAT_t CNorm2=1.0+CNorm2L1;

  const FLOAT_t ANorm=SQRT(ANorm2);
  const FLOAT_t BNorm=SQRT(BNorm2);
  const FLOAT_t CNorm=SQRT(CNorm2);

  // The three vectors in $\CC^3$ are:
  //
  //             ( PtA )
  //   A = ANorm (     )  , B = ..., C = ...
  //             (  1  )
  //
  // This is (A,B)(B,C)(C,A), the ``complex shape''
  const COMPLEX_t hCpx=
    ((1.0-PtAPtB)*(1.0-PtBPtC)*(1.0-PtCPtA))*ANorm2*BNorm2*CNorm2;
  // This is the same, LESS~1, calculated to avoid round off
  // error when PtA, PtB, and PtC are small
  const COMPLEX_t hCpxL1=
    -PtAPtB*(1.0-PtBPtC)*(1.0-PtCPtA)*ANorm2*BNorm2*CNorm2
    -PtBPtC*(1.0-PtCPtA)*ANorm2*BNorm2*CNorm2
    -PtCPtA*ANorm2*BNorm2*CNorm2
    +ANorm2L1*BNorm2*CNorm2
    +BNorm2L1*CNorm2
    +CNorm2L1;
  // Real and complex parts of (A,B)(B,C)(C,A)
  const FLOAT_t hL1=CREAL(hCpxL1),hi=CIMAG(hCpxL1);
  // This is Brehm's ``shape''
  const FLOAT_t h=1.0+hL1;

  // Squared sinh's of the side lengths
  FLOAT_t shSq0,shSq1,shSq2;
  // Symmetric polynomials in the squared sinh's
  FLOAT_t s0,s1,s2;

  // The entries of the companion matrix of the matrix of inner
  // products
  COMPLEX_t v[3][3];
  // The three squared absolute values of the off-diagonal entries of
  // the companion matrix of the matrix of inner products
  FLOAT_t pv01,pv12,pv20;
  // Symmetric polynomials in the squared absolute values of the
  // off-diagonal entries of the companion matrix of the matrix of
  // inner products
  FLOAT_t v0,v1,v2,v2Sq,v1Sq;
  // The cyclic product of the three off-diagonal entries of the companion
  // matrix
  COMPLEX_t vhCpx,vhCpxSq;
  // Real and imaginary parts of vhCpx
  FLOAT_t vh,vhi,vhSq;

  // The determinant of the matrix of inner products
  FLOAT_t d;

  // Equation for u
  FLOAT_t uEq[7],uEqComp[7];

  // Solutions for u
  FLOAT_t uSol[12],*uR,*uI;
  // Likewise, after polishing
  COMPLEX_t uVal[6];
  // Single solution and its square
  COMPLEX_t u,uSq;
  int uIx;
  // g-triple, list of g-triples
  COMPLEX_t g[3],gList[6][3];
  int gNum,gIx,NumTorus;
  // sin^2 of the circumradius (the answer)
  FLOAT_t sinh2R;

  int FoundCnt,BadFg;

  // Shows diagnostics
  inline void ShowDiagnostics() {
    FLOAT_t ASide=ASINH(SQRT(shSq0));
    FLOAT_t BSide=ASINH(SQRT(shSq1));
    FLOAT_t CSide=ASINH(SQRT(shSq2));
    printf("\nCircumradius D:\n");
    printf("  hi=%16.8g   h=%16.8g  hi/h=%16.8g\n",
	   hi,h,hi/h);
    printf("  d=%16.8g  hL1=%16.8g d/hL1=%16.8g\n",
	   d,hL1,d/hL1);
    printf("    a=%16.8g   b=%16.8g   c=%16.8g\n",
	   ASide,BSide,CSide);
    printf("  a-b=%16.8g b-c=%16.8g c-a=%16.8g\n",
	   ASide-BSide,BSide-CSide,CSide-ASide);
    printf("  g0*g1*g2="FOCnv FOPlCnv"*i\n",g[0]*g[1]*g[2]);
    printf("  g0*v21-v12/g0=%20.12g%+20.12g*i\n"
	   "  g1*v02-v20/g1=%20.12g%+20.12g*i\n"
	   "  g2*v10-v01/g2=%20.12g%+20.12g*i\n"
	   "            i*u=%20.12g%+20.12g*i\n",
	   g[0]*v[2][1]-v[1][2]/g[0],
	   g[1]*v[0][2]-v[2][0]/g[1],
	   g[2]*v[1][0]-v[0][1]/g[2],
	   I*u);
    printf("   g0=%20.12g%+20.12g*i\n"
	   "   g1=%20.12g%+20.12g*i\n"
	   "   g2=%20.12g%+20.12g*i\n",
	   g[0],g[1],g[2]);
    printf("  |g0|="FOCnv" |g1|="FOCnv" |g2|="FOCnv"\n",
	   CABS(g[0]),CABS(g[1]),CABS(g[2]));
  }

  BadFg=FALSE;

  // Squared sinh's of the side lengths.  The first one, for instance, is
  // given by:
  //
  //   (|PtB - PtC|^2 - |\det(PtB,PtC)|^2) / ((1 - |PtB|^2) (1-|PtC|^2))
  {
    Pt_t Diff; 
    FLOAT_t AbsDet;
    SubPts(Diff,PtB,PtC);
    AbsDet=CABS(PtB[0]*PtC[1]-PtB[1]*PtC[0]);
    shSq0=(PtAbsSq(Diff)-AbsDet*AbsDet)*BNorm2*CNorm2;
    SubPts(Diff,PtC,PtA);
    AbsDet=CABS(PtC[0]*PtA[1]-PtC[1]*PtA[0]);
    shSq1=(PtAbsSq(Diff)-AbsDet*AbsDet)*CNorm2*ANorm2;
    SubPts(Diff,PtA,PtB);
    AbsDet=CABS(PtA[0]*PtB[1]-PtA[1]*PtB[0]);
    shSq2=(PtAbsSq(Diff)-AbsDet*AbsDet)*ANorm2*BNorm2;

    v[0][0]=(-shSq0); v[1][1]=(-shSq1); v[2][2]=(-shSq2);
  }


  // The symmetric polynomials in the squared sinh's
  {
    s0=shSq0+shSq1+shSq2;
    s1=shSq0*shSq1+shSq1*shSq2+shSq2*shSq0;
    s2=shSq0*shSq1*shSq2;
  }

  // Certain other invariants
  d=(2.0*hL1-s0);
  vh=((2*hL1 - s0)*hL1 - s2);
  vhi=(-hi*d);
  vhCpx=vh+I*vhi;
  v0=(- 6.0*hL1 + 3.0*s0 + s1);
  v1=(- 4*(3*s0 + s1 - 3*hL1)*hL1 + (2*s1 + s2 + 3*s0)*s0);
  v2=(2*(2*(3*s0 + s1 - 2*hL1)*hL1 - (2*s1 + s2 + 3*s0)*s0)*hL1
      + (s1 + s2 + s0)*s0*s0 + s2*s2);
  vhSq=vh*vh; vhCpxSq=vhCpx*vhCpx; v1Sq=v1*v1; v2Sq=v2*v2;

  // Asymmetric values:
  // Off-diagonal entries of the companion matrix
  v[0][1]=ANorm*BNorm*(-CONJ(PtCPtA*(1-PtBPtC)*CNorm2)
		       -CONJ(PtBPtC)*CNorm2
		       +CNorm2L1
		       +PtAPtB);
  v[1][2]=BNorm*CNorm*(-CONJ(PtAPtB*(1-PtCPtA)*ANorm2)
		       -CONJ(PtCPtA)*ANorm2
		       +ANorm2L1
		       +PtBPtC);
  v[2][0]=CNorm*ANorm*(-CONJ(PtBPtC*(1-PtAPtB)*BNorm2)
		       -CONJ(PtAPtB)*BNorm2
		       +BNorm2L1
		       +PtCPtA);

  v[1][0]=CONJ(v[0][1]); v[2][1]=CONJ(v[1][2]); v[0][2]=CONJ(v[2][0]);

  // Their absolute values squared
  {
    FLOAT_t AbsTmp;
    AbsTmp=CABS(v[0][1]); pv01=AbsTmp*AbsTmp;
    AbsTmp=CABS(v[1][2]); pv12=AbsTmp*AbsTmp;
    AbsTmp=CABS(v[2][0]); pv20=AbsTmp*AbsTmp;
  }
			 

  //printf("Circumradius: ANorm=%20.12g BNorm=%20.12g CNorm=%20.12g\n",
  //       ANorm,BNorm,CNorm);
  //printf("Circumradius: PtAPtB=%20.12g PtBPtC=%20.12g PtCPtA=%20.12g\n",
  //	 PtAPtB,PtBPtC,PtCPtA);
  //
  //printf("Circumradius: a01=%25.15g%+25.15g*i\n"
  //       "            : a12=%25.15g%+25.15g*i\n"
  //       "            : a20=%25.15g%+25.15g*i\n",
  //       ANorm*BNorm*(1-PtAPtB),
  //       BNorm*CNorm*(1-PtBPtC),
  //       CNorm*ANorm*(1-PtCPtA));
  //printf("Circumradius: a01=%25.15g%+25.15g*i\n"
  //       "            : a12=%25.15g%+25.15g*i\n"
  //       "            : a20=%25.15g%+25.15g*i\n",
  //       ANorm*BNorm*(1-PtAPtB),
  //       BNorm*CNorm*(1-PtBPtC),
  //       CNorm*ANorm*(1-PtCPtA));
  //printf("Circumradius: v01=%25.15g%+25.15g*i\n"
  //       "            : v12=%25.15g%+25.15g*i\n"
  //       "            : v20=%25.15g%+25.15g*i\n",
  //       v[0][1],v[1][2],v[2][0]);

  // Check pv01, pv12, pv20 against v0, v1, v2
  {
    const FLOAT_t v0Tmp=pv01+pv12+pv20;
    const FLOAT_t v1Tmp=pv01*pv12+pv12*pv20+pv20*pv01;
    const FLOAT_t v2Tmp=pv01*pv12*pv20;
    //printf("Circumradius: v0=%20.12g v0Tmp=%20.12g\n"
    //	   "              v1=%20.12g v1Tmp=%20.12g\n"
    //	   "              v2=%20.12g v2Tmp=%20.12g\n",
    //	   v0,v0Tmp,v1,v1Tmp,v2,v2Tmp);
    assert(FABS(v0Tmp - v0) < epsBigger*v0);
    assert(FABS(v1Tmp - v1) < epsBigger*v1);
    assert(FABS(v2Tmp - v2) < epsBigger*v2);
  }

  // Check v01, v12, v20 against vhCpx
  {
    const COMPLEX_t vhCpxTmp=v[0][1]*v[1][2]*v[2][0];
    //printf("Circumradius:    vhCpx=%20.12g%+20.12g*i\n"
    //	   "            : vhCpxTmp=%20.12g%+20.12g*i\n",
    //	   vhCpx,vhCpxTmp);
    assert(CABS(vhCpxTmp - vhCpx) < epsBigger*CABS(vhCpx));
  }

  // See if triangle is colinear
  if(d<epsBigger*hL1 && FABS(hi)<epsBigger*FABS(h))
    return HUGE_FLOAT;

  // See if triangle is near colinear
  if(d<epsHuge*hL1 && FABS(hi)<epsHuge*FABS(h)) {
    printf("Circumradius: Almost colinear triangle:\n");
    printf("              hi=%16.8g   h=%16.8g  hi/h=%16.8g\n",
	   hi,h,hi/h);
    printf("               d=%16.8g hL1=%16.8g d/hL1=%16.8g\n",
	   d,hL1,d/hL1);
  }

  // Find possible values of $u$
  uEq[0]= 16*(vhSq*(vhSq - 2*v2) + v2Sq);
  uEq[1]= 0;
  uEq[2]= 8*v1*(vhSq - v2);
  uEq[3]= 8*vhi*(vhSq + v2);
  uEq[4]= - 4*vhSq*v0 + v1Sq;
  uEq[5]= - 2*vhi*v1;
  uEq[6]= v2;

  uEqComp[0]= 16*(vhSq*(vhSq + 2*v2) + v2Sq);
  uEqComp[1]= 0;
  uEqComp[2]= 8*v1*(vhSq + v2);
  uEqComp[3]= 8*FABS(vhi)*(vhSq + v2);
  uEqComp[4]= 4*vhSq*v0 + v1Sq;
  uEqComp[5]= 2*FABS(vhi)*v1;
  uEqComp[6]= v2;

  gsl_poly_complex_workspace *ws
    = gsl_poly_complex_workspace_alloc(7);

  gsl_poly_complex_solve(uEq,7,ws,uSol);

  gsl_poly_complex_workspace_free(ws);

  // Run through possible circumcenters
  uR=&(uSol[0]); uI=&(uSol[1]);
  gNum=0; // Number of g-triples found
  for(uIx=0;uIx<6;uIx++) {
    //if(FABS(*uI)<epsHuge*FABS(*uR)) {
    if(TRUE) {
      u=(*uR)+I*(*uI); uSq=u*u;

      // Calculates a g-value in terms of:
      //
      //   _v_, the corresponding off-diagonal entry in the inverse of
      //   the matrix of inner products,
      //
      //   _pv_, the squared absolute value of _v_,
      //
      //   _pvS_ the sum of the other two pv's,
      //
      //    and _pvP_, the product of the other two pv's.
      //
      // Also used are _vhCpx_, _v1_, _v2_, and _u_.  Different
      // values of _u_ correspond to different possibilities for the
      // circumcenter.
      inline COMPLEX_t gCalc(const COMPLEX_t v,const FLOAT_t pv,
			     const FLOAT_t pvS,const FLOAT_t pvP) {
	const FLOAT_t pvPSq=pvP*pvP;
	return
	  (v*u*(u*(u*(u*vhCpx*pvP + (vhCpxSq*pvS - pvPSq)*I)
		   - vhCpx*(vhCpxSq + v2)) - 2*pvP*(vhCpxSq - v2)*I)) /
	  ( - uSq*(u*vhCpx*v2*I - (vhCpxSq*pvS*pv - v2*pvP))
	    + vhCpxSq*(vhCpxSq - 2*v2) + v2Sq);
      }

      g[0]=gCalc(v[1][2],pv12,pv20+pv01,pv20*pv01);
      g[1]=gCalc(v[2][0],pv20,pv01+pv12,pv01*pv12);
      g[2]=gCalc(v[0][1],pv01,pv12+pv20,pv12*pv20);

      // Polish up the g-triple, and add it to the list of
      // possibilities (if necessary).
      gAdd(&gNum,gList,0,g,v);

      if(CABS(g[0]*g[1]*g[2]-1) > epsHuge
	 || CABS(g[0]*v[2][1] - v[1][2]/g[0] -I*u) > epsHuge*CABS(u)
	 || CABS(g[1]*v[0][2] - v[2][0]/g[1] -I*u) > epsHuge*CABS(u)
	 || CABS(g[2]*v[1][0] - v[0][1]/g[2] -I*u) > epsHuge*CABS(u)) {
	// ShowDiagnostics();
	BadFg=TRUE;
      }
    }
    uR+=2; uI+=2;
  }

  if(gNum<6) {
    COMPLEX_t gSimp[4][3]={{1,1,1},{1,-1,-1},{-1,1,-1},{-1,-1,1}};
    int gSimpIx;
    for(gSimpIx=0;gSimpIx<4;gSimpIx++)
      gAdd(&gNum,gList,0,gSimp[gSimpIx],v);
  }

  if(gNum<6) {
    COMPLEX_t gSol[2][3];
    int Ix0,Ix1,Ix2,SolIx1,SolIx2;
    uR=&(uSol[0]); uI=&(uSol[1]);
    for(uIx=0;uIx<6;uIx++) {
      COMPLEX_t a,c,SqrtDisc;
      const COMPLEX_t b=(-I*((*uR)+I*(*uI)));
      // Solve quadratic equations to get possible values of g[...]
      for(Ix0=0;Ix0<3;Ix0++) {
	Ix1=(Ix0+1) % 3; Ix2=(Ix0+2) % 3;
	a=v[Ix2][Ix1]; c=(-v[Ix1][Ix2]);
	SqrtDisc=CSQRT(b*b-4*a*c);
	gSol[0][Ix0]=(-b+SqrtDisc)/(2*a);
	gSol[1][Ix0]=(-b-SqrtDisc)/(2*a);
      }
      // Try these out in various combination.  Calculations are based
      // on solutions for two of the g-values, e.g. g1 and g2,
      // ignoring the third one, e.g. g0
      for(Ix0=0;Ix0<3;Ix0++) {
	Ix1=(Ix0+1) % 3; Ix2=(Ix0+2) % 3;
	for(SolIx1=0;SolIx1<2;SolIx1++)
	  for(SolIx2=0;SolIx2<2;SolIx2++) {
	    g[Ix1]=gSol[SolIx1][Ix1]; g[Ix2]=gSol[SolIx2][Ix2];
	    gAdd(&gNum,gList,Ix0,g,v);
	  }
      }
      uR+=2; uI+=2;
    }
  }
    
  FoundCnt=0; // Number of circumcenters found
  NumTorus=0;
  for(gIx=0;gIx<gNum;gIx++) {
    memcpy(g,gList[gIx],sizeof(g));
    if(FABS(CABS(g[0])-1) < epsBig &&
       FABS(CABS(g[1])-1) < epsBig &&
       FABS(CABS(g[2])-1) < epsBig) {
      const COMPLEX_t ACoeff=ANorm*(v[0][0]+g[2]*v[1][0]+CONJ(g[1])*v[2][0]);
      const COMPLEX_t BCoeff=BNorm*(v[0][1]+g[2]*v[1][1]+CONJ(g[1])*v[2][1]);
      const COMPLEX_t CCoeff=CNorm*(v[0][2]+g[2]*v[1][2]+CONJ(g[1])*v[2][2]);
      const FLOAT_t X3=ACoeff+BCoeff+CCoeff;

      NumTorus++;
      // The trial circumcenter
      COMPLEX_t XTrial[2]=
	{(ACoeff*PtA[0]+BCoeff*PtB[0]+CCoeff*PtC[0])/X3,
	 (ACoeff*PtA[1]+BCoeff*PtB[1]+CCoeff*PtC[1])/X3};
      // Is it in $B(\CC^2)?
      if(PtAbsSq(XTrial)<1.0) {
	const FLOAT_t sinh2B=QQ(PtB,XTrial)-1;
	const FLOAT_t sinh2C=QQ(PtC,XTrial)-1;
	sinh2R=QQ(PtA,XTrial)-1;

	FoundCnt++;
	memcpy(X,XTrial,sizeof(X));
	 printf("Circumradius S: sinh^2(d(A,X))=%23.16g\n"
	 	"                sinh^2(d(B,X))=%23.16g\n"
	 	"                sinh^2(d(C,X))=%23.16g\n",
	 	sinh2R,sinh2B,sinh2C);

	if(FABS(sinh2R-sinh2B)>epsHuge*sinh2R ||
	   FABS(sinh2B-sinh2C)>epsHuge*sinh2B) {
	  //printf("Circumradius: sinh^2(d(A,X))=%20.12g\n"
	  //	 "              sinh^2(d(B,X))=%20.12g\n"
	  //	 "              sinh^2(d(C,X))=%20.12g\n",
	  //	 sinh2R,sinh2B,sinh2C);
	  assert(FALSE); // <=== !!!
	  BadFg=TRUE;
	}
      }
    }
  }

  assert(gNum == 6);
  { 
  for(uIx=0;uIx<6;uIx++)
    uVal[uIx]=(gList[uIx][0]*v[2][1] - v[1][2]/gList[uIx][0])/I;
  uCheck(uEq,uEqComp,uVal);
  }
  assert(NumTorus == 4 || NumTorus == 6);
  assert(FoundCnt<=1); 

  NumTorusCnt[NumTorus]++;
  gNumCnt[gNum]++;
  if(BadFg) BadCnt++;

  if(FoundCnt>0) {
    GoodCnt++;
    return sinh2R;
  }
  else
    return HUGE_FLOAT;
}

// Applies Newton's method to the system:
//
//   g0*g1*g2 = 1
//   v12/g0-v21*g0 = v20/g1-v02*g1 = v01/g2-v10*g2
//
// for the variables g0, g1, and g2 in terms of the given parameters
// v12, v20, and v01 and their conjugates, v21, v02, and v10.
//
// The iteration starts with the given values of g1 and g2, ignoring
// the given value of g0.
//
int gPolish(COMPLEX_t *g0_p,COMPLEX_t *g1_p,COMPLEX_t *g2_p,
	     const COMPLEX_t v12,const COMPLEX_t v20,const COMPLEX_t v01) {
  const COMPLEX_t v21=CONJ(v12),v02=CONJ(v20),v10=CONJ(v01);
  COMPLEX_t g0,g1,g2;
  COMPLEX_t p,q,p1,p2,q1,q2,u,u1,u2;
  COMPLEX_t Det,g1Delta,g2Delta;
  FLOAT_t uAbs;
  int Cnt,DoneFg,Limit=30;

  DoneFg=FALSE;
  g1=*g1_p;
  g2=*g2_p;

  for(Cnt=0;!DoneFg && Cnt<Limit;Cnt++) {
    g0=1/(g1*g2);

    // Three quantities
    u=v12/g0-v21*g0;
    uAbs=CABS(u);
    p=(v01/g2-v10*g2)-u;
    q=(v20/g1-v02*g1)-u;
    if(CABS(p) < epsSmall*uAbs && CABS(q) < epsSmall*uAbs) {
      DoneFg=TRUE;
    }

    // Their derivatives w.r.t. g1
    u1=v12*g2+v21*g0/g1;
    p1=(-u1);
    q1=(-v20/(g1*g1)-v02)-u1;
    // Their derivatives w.r.t. g1
    u2=v12*g1+v21*g0/g2;
    p2=(-v01/(g2*g2)-v10)-u2;
    q2=(-u2);

    // Solve:
    //
    //   ( p1  p2 ) (g1Delta) = -p
    //   ( q1  q2 ) (g2Delta) = -q
    //
    Det=p1*q2-p2*q1;
    g1Delta=(-q2*p+p2*q)/Det;
    g2Delta=( q1*p-p1*q)/Det;

    g1+=g1Delta; g2+=g2Delta;

    if(CABS(g1Delta/g1) < epsSmall && CABS(g2Delta/g2) <epsSmall) {
      DoneFg=TRUE;
    }
    Cnt++;
  }
  //if(Cnt == 30)
  //  printf("gPolish: Cnt=%i\n",Cnt);
  *g0_p=1/(g1*g2); *g1_p=g1; *g2_p=g2;
  return DoneFg;
}

// Polishes a given triple of g-values, and then adds it to a list of
// g-value triples, if it isn't already there.
//
// Polishing is done on the basis of two of the given g-values, ignoring
// g[Ignore].
//
void gAdd(int *gNum_p,COMPLEX_t gList[6][3],const int Ignore,
	  COMPLEX_t g[3],COMPLEX_t v[3][3]) {
  int gTripIx,Ix,FoundFg,EqFg;
  int Ix0,Ix1,Ix2;

  Ix0=(Ignore) % 3; Ix1=(Ignore+1) % 3; Ix2=(Ignore+2)% 3;

  if(gPolish(&(g[Ix0]),&(g[Ix1]),&(g[Ix2]),
	     v[Ix1][Ix2],v[Ix2][Ix0],v[Ix0][Ix1])) {
    FoundFg=FALSE;
    for(gTripIx=0;gTripIx<*gNum_p;gTripIx++) {
      EqFg=TRUE;
      for(Ix=0;Ix<3;Ix++)
	if(CABS(g[Ix]/gList[gTripIx][Ix] - 1) > epsBig) {
	  EqFg=FALSE;
	  break;
	}
      if(EqFg) {
	FoundFg=TRUE;
	break;
      }
    }
    if(!FoundFg) {
      assert(*gNum_p<6);
      gList[*gNum_p][0]=g[0];
      gList[*gNum_p][1]=g[1];
      gList[*gNum_p][2]=g[2];
      (*gNum_p)++;
    }
  }
}

// Checks whether 6 values of u seem to be the 6 roots of uEq
void uCheck(FLOAT_t uEq[],FLOAT_t uEqComp[],COMPLEX_t uVal[]) {
  int Ix,Ix0,Ix1,Ix2,Ix3,Ix4,Ix5;
  COMPLEX_t sigma[6],Prod,Prod1,Prod2,Prod3;
  FLOAT_t AbsProd,Err;

  // Make the equation monic; likewise for the comparison values
  for(Ix=0;Ix<6;Ix++) {
    uEq[Ix]/=uEq[6];
    uEqComp[Ix]/=uEq[6];
  }
  
  // Zero out the symmetric polynomials
  memset(sigma,0,sizeof(sigma));

  // Calculate sigma[5]
  for(Ix0=0;Ix0<6;Ix0++) {
    Prod=uVal[Ix0];
    sigma[5]-=Prod;
    AbsProd=CABS(Prod);
    if(AbsProd > uEqComp[5]) uEqComp[5]=AbsProd;
  }

  // Calculate sigma[4]
  for(Ix0=0;Ix0<5;Ix0++)
    for(Ix1=Ix0+1;Ix1<6;Ix1++) {
      Prod=uVal[Ix0]*uVal[Ix1];
      sigma[4]+=Prod;
      AbsProd=CABS(Prod);
      if(AbsProd > uEqComp[4]) uEqComp[4]=AbsProd;
    }

  // Calculate sigma[3]
  for(Ix0=0;Ix0<4;Ix0++)
    for(Ix1=Ix0+1;Ix1<5;Ix1++) {
      Prod1=uVal[Ix0]*uVal[Ix1];
      for(Ix2=Ix1+1;Ix2<6;Ix2++) {
	Prod=Prod1*uVal[Ix2];
	sigma[3]-=Prod;
	AbsProd=CABS(Prod);
	if(AbsProd > uEqComp[3]) uEqComp[3]=AbsProd;
      }
    }

  // Calculate sigma[2]
  for(Ix0=0;Ix0<3;Ix0++)
    for(Ix1=Ix0+1;Ix1<4;Ix1++) {
      Prod1=uVal[Ix0]*uVal[Ix1];
      for(Ix2=Ix1+1;Ix2<5;Ix2++) {
	Prod2=Prod1*uVal[Ix2];
	for(Ix3=Ix2+1;Ix3<6;Ix3++) {
	  Prod=Prod2*uVal[Ix3];
	  sigma[2]+=Prod;
	  AbsProd=CABS(Prod);
	  if(AbsProd > uEqComp[2]) uEqComp[2]=AbsProd;
	}
      }
    }

  // Calculate sigma[1]
  for(Ix0=0;Ix0<2;Ix0++)
    for(Ix1=Ix0+1;Ix1<3;Ix1++) {
      Prod1=uVal[Ix0]*uVal[Ix1];
      for(Ix2=Ix1+1;Ix2<4;Ix2++) {
	Prod2=Prod1*uVal[Ix2];
	for(Ix3=Ix2+1;Ix3<5;Ix3++) {
	  Prod3=Prod2*uVal[Ix3];
	  for(Ix4=Ix3+1;Ix4<6;Ix4++) {
	    Prod=Prod3*uVal[Ix4];
	    sigma[1]-=Prod;
	    AbsProd=CABS(Prod);
	    if(AbsProd > uEqComp[1]) uEqComp[1]=AbsProd;
	  }
	}
      }
    }
  
  // Calculate sigma[0]
  sigma[0]=1;
  for(Ix=0;Ix<6;Ix++) sigma[0]*=uVal[0];
  AbsProd=CABS(sigma[0]);
  if(AbsProd > uEqComp[0]) uEqComp[0]=AbsProd;

  printf("uCheck: sigma=");
  for(Ix=0;Ix<6;Ix++)
    printf("%12.5g",CREAL(sigma[Ix]));
  printf("\n");

  printf("          uEq=");
  for(Ix=0;Ix<6;Ix++)
    printf("%12.5g",uEq[Ix]);
  printf("\n");

  printf("      uEqComp=");
  for(Ix=0;Ix<6;Ix++)
    printf("%12.5g",uEqComp[Ix]);
  printf("\n");

  printf("  Err/uEqComp=");
  for(Ix=0;Ix<6;Ix++) {
    Err=CABS(sigma[Ix]-uEq[Ix])/uEqComp[Ix];
    printf("%12.5g",Err);
    if(Err > sigmaMaxErr[Ix]) sigmaMaxErr[Ix]=Err;
    }
  printf("\n");

  printf("  sigmaMaxErr=");
  for(Ix=0;Ix<6;Ix++)
    printf("%12.5g",sigmaMaxErr[Ix]);
  printf("\n");


}