/* 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("%.14f (% .14f%+.14f*i,% .14f%+.14f*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(FILE *file,Pt_t Pts[],int N) {
  int Ix;

  Ix=0;
  for(Ix=0;Ix<N;Ix++) {
    OutputPt(file,Pts[Ix]);
  }
}

/* Output a point in the complex ball to a file. */
void OutputPt(FILE *file,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(file,"%.10f (% .15f%+.15f*i,% .15f%+.15f*i)\n",
	  Pt.s2,x1,y1,x2,y2);
  return;
}

/* Display a list of points in the complex ball */
void PrintElts(Elt_t Elts[],int N) {
  int Ix;

  Ix=0;
  for(Ix=0;Ix<N;Ix++) {
    PrintElt(Elts[Ix]);
  }
}

/* Display a point in the complex ball */
void PrintElt(Elt_t Elt) {
  PrintPt(Elt.p);
  PrintMtx(Elt.m);
  return;
}

/* Output a list of elements to a file */
void OutputElts(FILE *file,Elt_t Elts[],int N) {
  int Ix;

  Ix=0;
  for(Ix=0;Ix<N;Ix++) {
    OutputElt(file,Elts[Ix]);
  }
}

/* Output an element to a file */
void OutputElt(FILE *file,Elt_t Elt) {
  double x1,y1,x2,y2;

  OutputPt(file,Elt.p);
  OutputMtx(file,Elt.m);
  if(Elt.wLength==0)
    fprintf(file,"1\n");
  else
    fprintf(file,"%s\n",Elt.w);
  return;
}

/* Output the squared HS-norm of an element.  This would be for an
   FF-unitary element, _not_ for the normalized version of the
   matrix.  */
inline void OutputHS2(FILE *file,Elt_t Elt) {
  double sum,sum2,HS2,sumHS2;
  int ir,ic;
  
  // Calculate the squared HS-norm of the normalized matrix
  sum=sum2=0.0;
  for(ir=0;ir<3;ir++) 
    for(ic=0;ic<3;ic++) {
      sum+=cabs(Elt.m[ir][ic])*cabs(Elt.m[ir][ic]);
      sum2+=conj(Elt.m[ir][2])*FF[ir][ic]*Elt.m[ic][2];
    }
  // ... and adjust it so it would correspond to the FF-unitary
  // version of the matrix

  sumHS2=sum*FF[2][2]/sum2;

  HS2=3.0+(R0+1.0/R0)*(R0+1.0/R0)*(Elt.p.s2/R0)/(1.0-Elt.p.s2/R0);
  assert(fabs(sumHS2-HS2)<eps);

  fprintf(file,"%.14f\n",
	  3.0+(R0+1.0/R0)*(R0+1.0/R0)*(Elt.p.s2/R0)/(1.0-Elt.p.s2/R0));

  return;
  }

// Read a list of points in the complex ball of radius R0 from a file
int ReadPts(FILE *file,Pt_t Pts[],int MaxSiz) {
  int Ix;
  double z0R,z0I,z1R,z1I,s2;

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

/* Read a list of elements from a file */
int ReadElts(FILE *file,Elt_t Elts[],int MaxSiz) {
  int Ix,ir,ic;
  double z0R,z0I,z1R,z1I,s2,x,y;

  assert(file != NULL);
  Ix=0;
  while(fscanf(file,"%*[ \n\t]") != EOF) {

    assert(fscanf(file,"%lf (%lf%lf*i, %lf%lf*i) ",
		  &s2,&z0R,&z0I,&z1R,&z1I)==5);
    assert(Ix<MaxSiz);
    Elts[Ix].p.s2=s2;
    Elts[Ix].p.z[0]=z0R+z0I*I;
    Elts[Ix].p.z[1]=z1R+z1I*I;
    for(ir=0;ir<3;ir++)
      for(ic=0;ic<3;ic++) {
	assert(fscanf(file,"%lf %lf*i ",&x,&y)==2);
	Elts[Ix].m[ir][ic]=x+y*I;
      }
    assert(fscanf(file,ReadWordDirective,&Elts[Ix].w)==1);
    Elts[Ix].wLength=strlen(Elts[Ix].w);
    Ix++;
  }

  return Ix;
}

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

/* Output a matrix to a file. */
void OutputMtx(FILE *file, Mtx_t m) {
  double x,y;
  int ir,ic;
  for(ir=0;ir<3;ir++) {
    for(ic=0;ic<3;ic++) {
      x=creal(m[ir][ic]); y=cimag(m[ir][ic]);
      if(fabs(x)<eps) x=0.0; if(fabs(y)<eps) y=0.0;
      fprintf(file,"  % .15f%+.15f*i ",x,y);
    }
    fprintf(file,"\n");
  }
  return;
}

// Checks 2 points for equality
int EqPt(Pt_t p1,Pt_t p2) {
  int ix;

  for(ix=0;ix<2;ix++)
      if(cabs(p1.z[ix]-p2.z[ix]) > eps)return FALSE;
  return TRUE;
}

// Checks 2 points for equality
int EqPtEps2(Pt_t p1,Pt_t p2) {
  int ix;

  for(ix=0;ix<2;ix++)
      if(cabs(p1.z[ix]-p2.z[ix]) > eps2)return FALSE;
  return TRUE;
}

// Checks 2 matrices for equality
int EqMtx(Mtx_t m1,Mtx_t m2) {
  int ir,ic;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      if(cabs(m1[ir][ic]-m2[ir][ic]) > eps)return FALSE;
  return TRUE;
}

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

// Invert a matrix
void InvMtx(Mtx_t result,Mtx_t m) {
  double complex det,invdet;
  int ir,ic,ix;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      result[ir][ic]=
         m[(ic+1)%3][(ir+1)%3]*m[(ic+2)%3][(ir+2)%3]
	-m[(ic+1)%3][(ir+2)%3]*m[(ic+2)%3][(ir+1)%3];
  det=0;
  for(ix=0;ix<3;ix++)det+=m[1][ix]*result[ix][1];
  invdet=1.0/det;
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++)
      result[ir][ic]*=invdet;
  return;
}

/* Calculate the conjugate of 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;
}  

// Apply matrix, on the left, to a point of the complex ball.
inline void ApplyMtxToPt(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;
}


/* Apply a matrix, on the left, to a group element. */
void ApplyMtxToElt(Elt_t *newElt,Mtx_t mtx,Elt_t Elt) {
  double complex f;

  ApplyMtxToPt(&(newElt->p),mtx,Elt.p);
  MultMtx(newElt->m,mtx,Elt.m);
  NormalizeMtx(newElt->m);
  
  return;
}

/* Apply a generator, on the left to a group element. */
void ApplyGenToElt(Elt_t *newElt,int GenIx,Elt_t Elt,
		   Mtx_t Gens[],char letters[]) { 

  ApplyMtxToElt(newElt,Gens[GenIx],Elt);
  assert(Elt.wLength<MaxWordLength);
  (newElt->wLength)=Elt.wLength+1;
  (newElt->w)[0]=letters[GenIx];
  memcpy(&(newElt->w[1]),&(Elt.w[0]),Elt.wLength*sizeof(unsigned char));
  (newElt->w)[Elt.wLength+1]=0;
  
  return;
}

/* Normalize a matrix, mutiplying it by a scalar. */
inline void NormalizeMtx(Mtx_t mtx) {
  double complex f;
  int ir,ic;

  f=1.0/mtx[2][2];
  for(ir=0;ir<3;ir++)
    for(ic=0;ic<3;ic++) 
      mtx[ir][ic]*=f;
  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;
}


// See whether a given point is in the Dirichlet fundamental domain
// about zero
int InDirichlet(Pt_t P,int NumPts) {
  Pt_t P0;
  double Norm;
  int dummy;
  
  Norm=sqrt(PtAbsSq(P));
  if(Norm<eps)return TRUE;

  // Divide by the norm and obtain a unit vector
  P0.z[0]=P.z[0]/Norm;
  P0.z[1]=P.z[1]/Norm;
  P0.s2=PtAbsSq(P0);
  assert(fabs(P0.s2-1.0)<eps);
  
  return Norm < sVal(P0,NumPts,&dummy)+1000.0*eps;
}


/* 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 NumPts,int *PtIx_p) {
  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)/R0sq;
	// 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]));
  }
  *PtIx_p=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;
}


// This calculates, for two points in the complex ball of radius R0,
// $R0^2 \cosh^2(\hypdist(Pt1,Pt2))$.
inline double QQ(Pt_t Pt1,Pt_t Pt2) {
  double a1;

  a1=cabs(R0sq-Scal(Pt1,Pt2));
  return (a1*a1)/((R0sq-PtAbsSq(Pt1))*(R0sq-PtAbsSq(Pt2)));
}

inline double QQ1(Pt_t Pt1,Pt_t Pt2) {
  double a1;

  a1=cabs(R0sq-Scal(Pt1,Pt2));
  return (a1*a1)/(R0sq-PtAbsSq(Pt1));
}

// This calculates, for two points in the complex ball of radius R0,
// $R0 \tanh(\hypdist(Pt1,Pt2))$.  If one of the points is the origin,
// that value is the Euclidean norm of the other.
inline double sDist(Pt_t Pt1,Pt_t Pt2) {
  return R0*sqrt(1.0 - R0sq/QQ(Pt1,Pt2));
}

// Applies an isometry which sends the point $U$ to the point zero,
// and the point zero to the point $-U$ .  The point to be mapped is
// $V$. The calculated image is *VNew.  The map calculated is the
// inverse of the map calculated by _ZeroToU_.

inline void UToZero(Pt_t *VNew,Pt_t V,Pt_t U) {
  double scal,root,den,factV,factU;

  scal=Scal(V,U);
  root=sqrt(R0sq-U.s2);
  den=(R0sq-scal);
  factV=R0*root/den;
  factU=R0*(-R0+scal/(R0+root))/den;

  (VNew->z[0])= factV*V.z[0]+factU*U.z[0];
  (VNew->z[1])= factV*V.z[1]+factU*U.z[1];
  (VNew->s2)= PtAbsSq(*VNew);
  return;
}

// Applies an isometry which sends the point zero to the point $U$,
// and the point $-U$ to the point zero.  The point to be mapped is
// $V$. The calculated image is *VNew.  The map calculated is the
// inverse of the map calculated by _UToZero_.

inline void ZeroToU(Pt_t *VNew,Pt_t V,Pt_t U) {
  double scal,root,den,factV,factU;
  
  scal=Scal(V,U);
  root=sqrt(R0sq-U.s2);
  den=(R0sq+scal);
  factV=R0*root/den;
  factU=R0*(R0+scal/(R0+root))/den;

  (VNew->z[0])= factV*V.z[0]+factU*U.z[0];
  (VNew->z[1])= factV*V.z[1]+factU*U.z[1];
  (VNew->s2)= PtAbsSq(*VNew);
  return;
}

// Find the (hyperbolic) midpoint between two points of the
// (hyperbolic) complex ball of radius R0.  The formula can be deduced
// from the one on p.~70 of [Goldman, Complex Hyperbolic Geometry]

//   midpoint(U,V) =
//     \frac{  |R0^2-\langle U,V \rangle| \sqrt{R0^2-|V|^2} U
//           + (R0^2-\langle U,V \rangle) \sqrt{R0^2-|U|^2} V }
//          {  |R0^2-\langle U,V \rangle| \sqrt{R0^2-|V|^2}
//           + (R0^2-\langle U,V \rangle) \sqrt{R0^2-|U|^2} }

inline void Midpoint(Pt_t *result, Pt_t U, Pt_t V) {
  double complex VFact,UFact,den,scal;
  double halfQQ,UVQQ;  

  scal=R0sq-Scal(U,V);
  UFact=cabs(scal)*sqrt(R0sq-V.s2);
  VFact=scal*sqrt(R0sq-U.s2);
  den=(UFact+VFact);
  //printf("Midpoint A: scal=%f+i%f Ufact=%f\n",
  //	 creal(scal),cimag(scal),creal(UFact));
  //printf("Midpoint A: VFact=%f+i%f den=%f+i%f\n",
  //	 creal(VFact),cimag(VFact),creal(den),cimag(den));
  UFact/=den; VFact/=den;
  //printf("Midpoint B: UFact=%f+i%f Vfact=%f+i%f\n",
  //	 creal(UFact),cimag(UFact),creal(VFact),cimag(VFact));

  (result->z[0]) = UFact*U.z[0]+VFact*V.z[0];
  (result->z[1]) = UFact*U.z[1]+VFact*V.z[1];
  (result->s2) = PtAbsSq(*result);

  halfQQ=QQ(*result,U);
  UVQQ=QQ(U,V);
  //printf("Midpoint C: U, V, *result\n");
  //PrintPt(U);
  //PrintPt(V);
  //PrintPt(*result);
  assert(fabs(QQ(*result,V)-halfQQ) < eps);
  assert(fabs(R0sq*(2.0/R0sq*halfQQ-1.0)*(2.0/R0sq*halfQQ-1.0)
	      -UVQQ) < eps);

  return;
}

// Check whether a matrix has a fixed point for its action on the
// complex ball.
//
// Assume as usual that the matrix is unitary relative to
// $FF=\diag(1,1,-R_0^2)$, so the complex ball is of radius $R_0$
//
// Return 0 if there is no fixed point.
// Return -1 if the matrix has a fixed point, but is not of finite order.
// Otherwise return the order of the matrix.
//
// Find a fixed point when it exists (or they exist).
//
int Order(Mtx_t Mtx,Pt_t *Fixed_p) {
  Mtx_t PowMtx,tmpMtx;
  double trace,sum;
  int n,order,ix;
  Pt_t V,M;
  
  order=(-1);
  memcpy(PowMtx,Mtx,sizeof(Mtx_t));
  NormalizeMtx(PowMtx);
  for(n=1;n<=600;n++) {
    //printf("Order A: PowMtx ---\n");
    //PrintMtx(PowMtx);
    if(EqMtx(PowMtx,Id3)) {order=n; break;}

    // For an F-unitary matrix M, a necessary condition for having a
    // fixed point is that $|\tr M|<3$.  For a projectively unitary
    // matrix M such that $M^* F M = c^2 F$, the correct condition is
    // $|\tr M|/c < 3$, that is $|\tr M|^2 / c^2 < 9$.
    trace=sum=0;
    for(ix=0;ix<3;ix++) {
      trace+=PowMtx[ix][ix];
      sum+=conj(PowMtx[ix][2])*FF[ix][ix]*PowMtx[ix][2];
    }
    //printf("Order B: trace=%f, sum=%f, FF[2][2]=%f, PowMtx ---\n",
    // trace,sum,FF[2][2]);
    //PrintMtx(PowMtx);
    if( cabs(sum)<eps ||
	cabs(trace*trace / (sum / (FF[2][2]))) > 9.0+eps) {
      order=0;
      break;
    }
    MultMtx(tmpMtx,Mtx,PowMtx);
    memcpy(PowMtx,tmpMtx,sizeof(Mtx_t));
    NormalizeMtx(PowMtx);
  }

  if(order == 0) return order;

  if(order > 0) {
    // Find a fixed point
    memset(Fixed_p,0,sizeof(Pt_t)); // Trial point
    while(TRUE) {
      ApplyMtxToPt(&V,Mtx,*Fixed_p); // Shift our trial point by the
				     // matrix
      //printf("Order F: EqPt(*Fixed_p,V)=%i *Fixed_p, V, M\n",
      //	     EqPt(*Fixed_p,V));
      //PrintPt(*Fixed_p);
      //PrintPt(V);
      //PrintPt(M);
      if(EqPtEps2(*Fixed_p,V))return order; // If it's fixed, it's
					    // fixed
      Midpoint(&M,*Fixed_p,V); // Calculate the point midway between
			       // the trial point and its shift
      
      memcpy(Fixed_p,&M,sizeof(Pt_t)); // That's the new trial point
    }
  }

  // At this point we assume that our matrix belongs to a discrete
  // group.  Therefore, it ought not to fix a point without being of
  // finite order.  We're also assuming the order is $\leq 600$.
  // Given all that, it's unreasonable that we don't have an answer by
  // now.

  assert(FALSE);
  return -1;
}

/* 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;
  
  // Compare the squared Euclidean distances of the two points
  if((Pt1->s2)>(Pt2->s2)+eps)return +1;
  if((Pt2->s2)>(Pt1->s2)+eps)return -1;

  // Then compare the coordinates of the two points
  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;
}

/* Compare two group elements: the order to use is somewhat arbitrary.
   Here it is lexocographic, with the first number tested being
   distance from the origin. */

int CompareElts(const void *Ptr1,const void *Ptr2) {
  const Elt_t *Elt1=(const Elt_t *)Ptr1, *Elt2=(const Elt_t *)Ptr2;
  int Ix,Row,Col,tmp;
   
  // Compare the squared Euclidean distances of $g_1(O)$ to~$O$, and
  // $g_2(O)$ to~$O$, where $g_1$ and $g_2$ are the two elements to be
  // compared;
  // then compare the coordinates of $g_1(O)$ to those of $g_2(O)$;
  // then compare the word lengths of the two elements.

  tmp=ComparePts(&(Elt1->p),&(Elt2->p));
  if(tmp != 0)return tmp;

  // Compare the matrices of $g_1$ and $g_2$
  for(Row=0;Row<3;Row++)
    for(Col=0;Col<3;Col++) {
    if(creal((Elt1->m)[Row][Col])>creal((Elt2->m)[Row][Col])+eps)return +1;
    if(creal((Elt2->m)[Row][Col])>creal((Elt1->m)[Row][Col])+eps)return -1;

    if(cimag((Elt1->m)[Row][Col])>cimag((Elt2->m)[Row][Col])+eps)return +1;
    if(cimag((Elt2->m)[Row][Col])>cimag((Elt1->m)[Row][Col])+eps)return -1;
    }

  return 0;
}

/* Compare two group elements: the order to use is somewhat arbitrary.
   Here it is lexocographic, with the first number tested being
   distance from the origin.  If the two elements are equal as
   matrices, one then compares the word lengths of the given words.*/
int CompareEltsWithLength(const void *Ptr1,const void *Ptr2) {
  const Elt_t *Elt1=(const Elt_t *)Ptr1, *Elt2=(const Elt_t *)Ptr2;
  int tmp;
  
  tmp=CompareElts(Elt1,Elt2);
  if(tmp != 0)return tmp;

  // Compare the word lengths
  if(Elt1->wLength > Elt2->wLength)return +1;
  if(Elt2->wLength > Elt1->wLength)return -1;

  return 0;
}