import java.util.*;

/**
 * A class which computes the combinatorics of a fundamental domain.
 *
 * @author Jonathan Poritz
 * @version 1.7, 17 Sep 1998
 */

public class DomainCombinatorics {
  /**
   * A Vector of BoundaryArcs which form the boundary of the
   * current fundamentaldomain
   */
  Vector boundary;

  /**
   * The precision we require to consider two circles in the plane
   * to be overlapping.
   */
  public static final double CIRCLAPRES = .000001;

  /**
   * The precision we require to consider two angles to be different.
   */
  public static final double ANGDIFPRES = .000001;

  /**
   * Constructor for DomainCombinatorics's using only the trace, axis
   * distance and domain type.
   *
   * @param trace the trace whose domain we are to examine
   * @param dist the axis distance for the domain in question
   * @param domaintype the domaintype, given by the selected index of
   *                   <b>outputpanel</b>'s Choice <b>fordordirichlet</b>
   */
  public DomainCombinatorics(Complex trace, double dist, int domaintype) {
    this(trace, dist, domaintype, Complex.make_z11(trace));
  }

  /**
   * Constructor for DomainCombinatorics's using only the trace, axis
   * distance, domain type and (1,1) element of the generator.
   *
   * @param trace the trace whose domain we are to examine
   * @param dist the axis distance for the domain in question
   * @param domaintype the domaintype, given by the selected index of
   *                   <b>outputpanel</b>'s Choice <b>fordordirichlet</b>
   * @param z11 the (1,1) element of the generator matrix at the given trace
   */
  public DomainCombinatorics(Complex trace, double dist, int domaintype, Complex z11) {
    this(trace, dist, domaintype, z11, Complex.make_maxpow(trace, z11));
  }

  /**
   * Constructor for DomainCombinatorics's using the trace, axis distance,
   * domain type, (1,1) element of the generator and corresponding maxpow.
   *
   * @param trace the trace whose domain we are to examine
   * @param dist the axis distance for the domain in question
   * @param domaintype the domaintype, given by the selected index of
   *                   <b>outputpanel</b>'s Choice <b>fordordirichlet</b>
   * @param z11 the (1,1) element of the generator matrix at the given trace
   * @param maxpow the maxpow corresponding to the other input data
   */
  public DomainCombinatorics(Complex trace, double dist, int domaintype, Complex z11, int maxpow) {
    // create the vector to hold the boundary's BoundaryArcs
    boundary = new Vector();
    // the +/- 1 circles are a good starting point, we'll always need them
    BoundaryArc plusone;
    BoundaryArc minusone;
    if (trace.y == 0D) {
      // here we are either elliptic, parabolic or hyperbolic
      // make an empty BoundaryArc
      plusone = new BoundaryArc();
      if (Math.abs(trace.x) >= 2D) // hyerbol or parabol: only need indx 1
	plusone.index = 1;
      else {
	// here be elliptic elements
	// this z11 is e^theta, compute theta
	double theta = Math.atan2(-z11.y, z11.x);
	// now find the power of z11 for which the maximum height of the
	// corresponding circle is greatest, up to the power maxpow
	double maxheight = 0D;
	int maxheightindex = 0;
	for (double k=1; k <= maxpow; k++) {
	  double ktheta = k*theta;
	  double skt = Math.sin(ktheta);
	  double ckt = Math.cos(ktheta);
	  // try both the positive and negative powers, to see which
	  // gives the greatest height
	  double height = Math.abs(1D/skt)-(ckt/skt);
	  if (height > maxheight) {
	    maxheight = height;
	    maxheightindex = (int) k;
	  }
	  height = Math.abs(-1D/skt)+(ckt/skt);
	  if (height > maxheight) {
	    maxheight = height;
	    maxheightindex = -((int) -k);
	  }
	}
	// here, then, is the index of the highest circle, note that
	// this may not be a positive index, despite the name
	plusone.index = maxheightindex;
      }
      // make a negative copy of the plusone BoundaryArc, put both
      // in the Vector boundary and return
      minusone = BoundaryArc.makeNeg(plusone);
      boundary.addElement(plusone);
      boundary.addElement(minusone);
      return;
    }
    // now loxodromic elements.  make the +/- 1 BoundaryArcs
    // start and end angles are set to the nonsensical values of 0 and -1;
    // these angles should either never be important, if the domain is of
    // type {1}, or should be reset to reasonable values if there are
    // intersections
    plusone = new BoundaryArc(z11, dist, 1, domaintype, 0D, -1D);
    minusone = BoundaryArc.makeNeg(plusone);
    // put them in the boundary Vector
    boundary.addElement(plusone);
    boundary.addElement(minusone);
    // if we are outside the circle of radius 2 in the trace plane or have
    // zero axis distance in a Dirichlet domain, +/- 1 is the entire boundary
    if ((trace.modulus() >= 2D) || ((domaintype == 1) && (dist == 0D)))
      return;
    // also the entire boundary is the +/- 1 circles if they do not intersect
    double cendist = BoundaryCircle.distance(plusone, minusone);
    if (cendist > 2*plusone.rad)
      return;
    double angoffset = Math.acos(cendist/(2D*plusone.rad));
    double cenlineang = Math.atan2((plusone.cen.y - minusone.cen.y),(plusone.cen.x - minusone.cen.x));
    // +/- 1 boundary circles intersect, set the beginning and end of their
    // arcs which are outside their negatives
    minusone.start = fixAngle(cenlineang + angoffset);
    minusone.end = minusone.start + 2D*Math.PI - 2D*angoffset;
    plusone.start = fixAngle(minusone.start + Math.PI);
    plusone.end = plusone.start + 2D*Math.PI - 2D*angoffset;
    Complex z11powers = new Complex(z11);
    // we'll now loop through all powers of the generator up to maxpow and
    // see if the corresponding positive or negative boundary circle
    // contributes to the boundary
g:  for(int i=2; i <= maxpow; i++) {
      // make the next power of z11, the corresp pos boundary arc and its neg
      z11powers = Complex.multiply(z11powers, z11);
      BoundaryArc nextpos = new BoundaryArc(z11powers, dist, i, domaintype, 0, -1);
      BoundaryArc nextneg = BoundaryArc.makeNeg(nextpos);
      // the following are the number and values of the angles on the
      // boundary circle nextpos where it intersects other boundary circles
      int anglecnt = 0;
      double[] angles = new double[2*boundary.size()+3];
      // first see if and where nextpos intersects nextneg
      cendist = BoundaryCircle.distance(nextpos, nextneg);
      if (cendist <= 2*nextpos.rad) {
        angoffset = Math.acos(cendist/(2D*nextpos.rad));
        cenlineang = Math.atan2((nextneg.cen.y - nextpos.cen.y),(nextneg.cen.x - nextpos.cen.x));
        if (angoffset == 0D) {
          angles[0] = cenlineang;
          anglecnt = 1;
        }
        else {
          angles[0] = cenlineang - angoffset;
          angles[1] = cenlineang + angoffset;
          anglecnt = 2;
        }
      }
      // next loop through the current list of boundary circles to
      // see if and where they intersect nextpos
      for(int j=0; j < boundary.size(); j++) {
	// working on the jth BoundaryArc in the Vector boundary
        BoundaryArc arcj = ((BoundaryArc) boundary.elementAt(j));
        double cendistnptoj = BoundaryCircle.distance(nextpos, arcj);
	// if arcj contains nextpos completely, arcj does not contribute to
	// the boundary at all, and we can go to the next step of the loop g
        if ((cendistnptoj + nextpos.rad) <= arcj.rad)
          continue g;
	// if nextpos contains arcj or they are too far apart to intersect,
	// continue the current loop
        if (((cendistnptoj + arcj.rad) < nextpos.rad) || ((arcj.rad + nextpos.rad) < cendistnptoj))
          continue;
	// have at least one intersection; compute and put in array of angles
        angoffset = Math.acos((nextpos.rad*nextpos.rad+cendistnptoj*cendistnptoj-arcj.rad*arcj.rad)/(2D*nextpos.rad*cendistnptoj));
        cenlineang = Math.atan2((arcj.cen.y - nextpos.cen.y),(arcj.cen.x - nextpos.cen.x));
        if (angoffset == 0D)
          angles[anglecnt++] = cenlineang;
        else {
          angles[anglecnt++] = cenlineang - angoffset;
          angles[anglecnt++] = cenlineang + angoffset;
        }
      }
      // fix all the angles we just made
      for (int j=0; j<anglecnt; j++)
	angles[j] = fixAngle(angles[j]);
      // bubblesort the angles in increasing order
      for (int j=anglecnt; --j>=0; ) {
        for (int jj=0; jj<j; jj++) {
          if (angles[jj] > angles[jj+1]) {
            double switcheroo = angles[jj];
            angles[jj] = angles[jj+1];
            angles[jj+1] = switcheroo;
          }
        }
      }
      // put on a last angle, so that between successive angles we have
      // arcs of nextpos which do not intersect any other boundary circle
      angles[anglecnt] = angles[0]+2*Math.PI;
      // we will loop through and see if these arcs of nextpos lie inside
      // some boundary circle and thus can be discarded from the boundary
      boolean[] discard = new boolean[anglecnt];
      for (int j=0; j<anglecnt; j++) {
	// examine each arc of nextpos by computing the location of its midpt
        double midangle = (angles[j]+angles[j+1])/2D;
        Complex midpoint=new Complex(nextpos.cen.x+nextpos.rad*Math.cos(midangle), nextpos.cen.y+nextpos.rad*Math.sin(midangle));
        if (Complex.distsq(midpoint, nextneg.cen) > nextneg.rad*nextneg.rad) {
	  // if this arc not contained in nextneg, it may contribute to the bdy
          discard[j] = false;
          for (int jj=0; jj < boundary.size(); jj++) {
	    // now check if this arc is inside any of the circles
	    // already known to contribute to the boundary
            BoundaryCircle circlejj = ((BoundaryCircle) boundary.elementAt(jj));
            if (Complex.distsq(midpoint, circlejj.cen) < circlejj.rad*circlejj.rad) {
	      // yes, this arc is hidden, set to discard and break out of loop
              discard[j] = true;
              break;
            }
          }
        }
        else {
	  // this arc hidden inside nextneg, discard
          discard[j] = true;
        }
      }
      // start putting together the next idea of the boundary
      Vector newboundary = new Vector();
      for (int j=0; j<anglecnt; j++) {
        if (!discard[j]) {
	  // for every non-discarded arc of nextpos, make a copy with the
	  // correct angles and add to the Vector newboundary
          BoundaryArc next = nextpos.copyArc();
          next.start = angles[j];
          next.end = angles[j+1];
          newboundary.addElement(next);
	  // likewise for nextneg
          next = nextneg.copyArc();
          next.start = fixAngle(angles[j]+Math.PI);
          next.end = next.start + angles[j+1] - angles[j];
          newboundary.addElement(next);
        }
      }
      // if there was nothing saved, continue the g loop
      if (newboundary.size() == 0)
	continue;
      // now that we have some new arcs, we may have to discard part or all
      // of one of circles we previously thought were on the boundary
      for (int j=0; j < boundary.size(); j++) {
	// as before, take the next arc and compute the number and values
	// of its intersects with other arcs
	BoundaryArc arcj = ((BoundaryArc) boundary.elementAt(j));
	anglecnt = 0;
	angles = new double[6];
	// here we look for intersections with nextpos
	double cendistnptoj = BoundaryCircle.distance(nextpos, arcj);
        if (((cendistnptoj + nextpos.rad) >= arcj.rad) && ((cendistnptoj + arcj.rad) >= nextpos.rad) && ((arcj.rad + nextpos.rad) >= cendistnptoj)) {
	  angoffset = Math.acos((arcj.rad*arcj.rad+cendistnptoj*cendistnptoj-nextpos.rad*nextpos.rad)/(2D*arcj.rad*cendistnptoj));
	  cenlineang = Math.atan2((nextpos.cen.y - arcj.cen.y),(nextpos.cen.x - arcj.cen.x));
	  if (angoffset == 0D)
	    angles[anglecnt++] = cenlineang;
	  else {
	    angles[anglecnt++] = cenlineang - angoffset;
	    angles[anglecnt++] = cenlineang + angoffset;
	  }
	}
	// here we look for intersections with nextneg
	double cendistnntoj = BoundaryCircle.distance(nextneg, arcj);
        if (((cendistnntoj + nextneg.rad) >= arcj.rad) && ((cendistnntoj + arcj.rad) >= nextneg.rad) && ((arcj.rad + nextneg.rad) >= cendistnntoj)) {
	  angoffset = Math.acos((arcj.rad*arcj.rad+cendistnntoj*cendistnntoj-nextneg.rad*nextneg.rad)/(2D*arcj.rad*cendistnntoj));
	  cenlineang = Math.atan2((nextneg.cen.y - arcj.cen.y),(nextneg.cen.x - arcj.cen.x));
	  if (angoffset == 0D)
	    angles[anglecnt++] = cenlineang;
	  else {
	    angles[anglecnt++] = cenlineang - angoffset;
	    angles[anglecnt++] = cenlineang + angoffset;
	  }
	}
	// now loop through and fix all these angles.  this involves putting
	// them in the interval [0,2Pi), then seeing if they are in the
	// range (start,end) of arcj, or likewise with the angle + 2Pi
	int newj = 0;
	for (int jj=0; jj<anglecnt; jj++) {
	  double currangle = fixAngle(angles[jj]);
	  if ((currangle > arcj.start) && (currangle < arcj.end)) {
	    angles[newj++] = currangle;
	    continue;
	  }
	  currangle += 2*Math.PI;
	  if ((currangle > arcj.start) && (currangle < arcj.end))
	    angles[newj++] = currangle;
	}
	anglecnt = newj;
	// add the beginning and end of arcj and then bubblesort
	// the angles in increasing order
	angles[anglecnt++] = arcj.start;
	angles[anglecnt++] = arcj.end;
	for (int jj=anglecnt; --jj>=0; ) {
	  for (int jjj=0; jjj<jj; jjj++) {
	    if (angles[jjj] > angles[jjj+1]) {
	      double switcheroo = angles[jjj];
	      angles[jjj] = angles[jjj+1];
	      angles[jjj+1] = switcheroo;
	    }
	  }
	}
	// now remove any possible duplicate angles
	newj = 1;
	for (int jj=1; jj<anglecnt; jj++) {
	  double currangle = angles[jj];
	  if (Math.abs(currangle-angles[newj-1]) > ANGDIFPRES)
	    angles[newj++] = currangle;
	}
	anglecnt = newj - 1;
	// now we can compute midpoints of all of the subarcs of arcj
	// whose endpoints are successive angles of the array angles
	for (int jj=0; jj<anglecnt; jj++) {
	  double midangle = (angles[jj]+angles[jj+1])/2D;
	  Complex midpoint = new Complex(arcj.cen.x+arcj.rad*Math.cos(midangle), arcj.cen.y+arcj.rad*Math.sin(midangle));
	  if ((Complex.distsq(midpoint, nextneg.cen) > nextneg.rad*nextneg.rad) && (Complex.distsq(midpoint, nextpos.cen) > nextpos.rad*nextpos.rad)) {
	    // if the midpoint is not inside nextpos or nextneg this subarc
	    // does contribute to the bdy, so add it to the Vector newboundary
	    BoundaryArc next = arcj.copyArc();
	    next.start = fixAngle(angles[jj]);
	    next.end = next.start + angles[jj+1] - angles[jj];;
	    newboundary.addElement(next);
	  }
	}
      }
      // done fixing up newboundary, use it as the actual boundary
      boundary = newboundary;
    }
  }

  /**
   * Make a string suitable for use in <b>inputpanel</b>'s
   * <b>currdomaincomblabel</b> describing this domain combinatorics.
   * The string is either the domain type in the sense of Drumm and Poritz
   * or a comma-separated list of indices (beginning and ending with two
   * question marks and enclosed in square brackets) of sides that have been
   * computed to appear on the boundary of the fundamental domain at
   * infinity.
   *
   * @return the domain type or error side listing string
   */
  public String makeLabel() {
    // how many sides in the boundary
    int siz = boundary.size();
    // keep a list of just the indices of the sides
    int[] indices = new int[siz];
    // and the location in the Vector boundary of these indices
    int[] location = new int[siz];
    for(int i=0; i < siz; i++) {
      indices[i] = ((BoundaryArc) boundary.elementAt(i)).index;
      location[i] = i;
    }
    if (siz == 2) {
      // if two faces, one should be the negative of the other;
      // if not, convert boundary to String and return with questions
      if (indices[1] != -indices[0])
	return "?? "+boundary.toString()+" ??";
      else
	return "{"+Math.abs(indices[0])+"}";
    }
    // should always have at least an even number of faces;
    // if not, convert boundary to String and return with questions
    if (siz != 2*((int) Math.floor(((double) siz)/2D)))
      return "?? "+boundary.toString()+" ??";
    for (int i=siz; --i >= 0; ) {
      //bubblesort the indices, keep the location of these indices
      // in the Vector boundary up to date
      for (int j=0; j < i; j++) {
	if (indices[j] > indices[j+1]) {
	  int switcheroo = indices[j];
	  indices[j] = indices[j+1];
	  indices[j+1] = switcheroo;
	  switcheroo = location[j];
	  location[j] = location[j+1];
	  location[j+1] = switcheroo;
	}
      }
    }
    // if indices not appearing in positive/negative pairs,
    // convert boundary to String and return with questions
    for (int jj=0; jj < (siz/2); jj++)
      if (indices[jj] != - indices[siz - jj - 1])
	return "?? "+boundary.toString()+" ??";
    int j = siz-1;
    // find now the location, as a Complex number, of the counterclockwise
    // end of the arc corresponding to the biggest index
    BoundaryArc arcbig = (BoundaryArc) boundary.elementAt(location[j]);
    Complex bigend = new Complex(arcbig.cen.x+arcbig.rad*Math.cos(arcbig.end), arcbig.cen.y+arcbig.rad*Math.sin(arcbig.end));
    for ( ; --j >= 0; ) {
      // break if bigend is on the circle of the jth arc; this means
      // that the counterclockwise successor to arcbig is arcj
      BoundaryArc arcj = (BoundaryArc) boundary.elementAt(location[j]);
      if (Math.abs(Complex.distsq(arcj.cen, bigend)-arcj.rad*arcj.rad)<CIRCLAPRES)
	break;
    }
    // j should be positive; if not, the endpoint of the BoundaryArc of
    // biggest index is on no other BoundaryCircle, so convert boundary
    // to String and return with questions
    if (j < 0)
      return "?? "+boundary.toString()+" ??";
    switch (siz) {
    case 4:
      // if four sides in the Vector boundary, make the appropriate type
      // string, with the positve side indices appearing the order that
      // they arise as one traverses the boundary of the fundamental
      // domain in a clockwise fashion
      if (indices[j] < 0)
	return "{"+(-indices[j])+","+indices[3]+"}";
      else
	return "{"+indices[3]+","+indices[j]+"}";
    case 6:
      // we're expecting {n,n+m,m}, so convert boundary to string and
      // return with questionsmarks if the counterclockwise successor
      // to arcbig is negative, leave no room for the n or m; otherwise
      // return the type string
      if ((indices[j] < 0) || (indices[5] != (indices[3]+indices[4])))
	return "?? "+boundary.toString()+" ??";
      else
	return "{"+indices[7-j]+","+indices[5]+","+indices[j]+"}";
    default:
      // only get here if error, so convert boundary to String and return
      // with question marks
      return "?? "+boundary.toString()+" ??";
    }
  }

  /**
   * Converts angle variables that are in the interval [-2Pi,4Pi) to lie
   * in the interval [0,2Pi).
   *
   * @param angle the input angle to convert
   * @return the angle adjusted to lie in [0,2Pi)
   */
  public static double fixAngle(double angle) {
    if (angle >= 2D*Math.PI)
      return angle - 2D*Math.PI;
    if (angle < 0D)
      return angle + 2D*Math.PI;
    return angle;
  }

}