import java.util.StringTokenizer;

/**
 * A class representing complex numbers, implementing calculations with them
 * and providing some numerical computations for fundamental domains.
 *
 * @author Jonathan A. Poritz
 * @version 1.7, 10 Sep 1998
 */

public class Complex extends Object {
  /**
   * The real part of this complex number.
   */
  public double x;
  /**
   * The imaginary part of this complex number.
   */
  public double y;

  /**
   * The complex constant 0.
   */
  public static final Complex ZERO = new Complex(0D, 0D);
  /**
   * The complex constant 1.
   */
  public static final Complex ONE = new Complex(1D, 0D);
  /**
   * The complex constant 2.
   */
  public static final Complex TWO = new Complex(2D,0D);
  /**
   * The complex constant 4.
   */
  public static final Complex FOUR = new Complex(4D,0D);
  /**
   * The complex constant -1.
   */
  public static final Complex MINUS_ONE = new Complex(-1D, 0D);
  /**
   * The complex constant i.
   */
  public static final Complex I = new Complex(0D, 1D);
  /**
   * The complex constant -i.
   */
  public static final Complex MINUS_I = new Complex(0D, -1D);

  /**
   * The cutoff for calculated maxpows -- above this value we arbitrarily
   * truncate down to this.  A good default is around 50; much more and
   * computations start to take way too long, much less and a reasonable
   * window on the trace plane includes lots of pixel for which we use the
   * artificial cutoff.
   */
  public static final int MAXPOWMAX = 51;
  /**
   * An cutoff for the calculation of orders of elliptic elements.  We only
   * try possible orders up to this value, then give up.
   */
  public static final int MAXELLIPORD = 1000;
  /**
   * Powers of an elliptic generator whose (1,1) entry are at this distance
   * squared or less are considered to be the identity, so enabling us to
   * compute the order of that elliptic element.
   */
  public static final double IDENTRESOL = .000001D;

  /**
   * A string containing one period for formatting <b>Double.toString</b>
   * output to a certain number of decimal digits.
   */
  public static final String period = ".";
  /**
   * The default number of digits to display in formatted coversion of
   * <b>Double</b>s to <b>String</b>s.
   */
  public static final int DOUBLEDIGITS = 7;

  /**
   * Constructor for the complex number zero.
   */
  public Complex() { 
    x = 0D; 
    y = 0D; 
  }

  /**
   * Constructor for an arbitrary complex number.
   *
   * @param xx the real part
   * @param yy the imaginary part
   */
  public Complex(double xx, double yy) { 
    x = xx; 
    y = yy; 
  }

  /**
   * Constructor for a complex number identical to one given as input.
   * @param z the complex number to duplicate
   */
  public Complex(Complex z) { 
    x = z.x; 
    y = z.y; 
  }

  /**
   * Compute the sum of two complex numbers.
   *
   * @param a one input complex number
   * @param b and another
   * @return the sum
   */
  public static Complex add(Complex a, Complex b) {
    Complex c = new Complex(a.x + b.x, a.y + b.y);
    return c;
  }

  /**
   * Compute the difference of two complex numbers.
   *
   * @param a one input complex number
   * @param b and another
   * @return the difference
   */
  public static Complex subtract(Complex a, Complex b) {
    Complex c = new Complex(a.x - b.x, a.y - b.y);
    return c;
  }

  /**
   * Multiplies one input complex number by another.
   *
   * @param a one input complex number
   * @param b and another
   * @return the product
   */
  public static Complex multiply(Complex a, Complex b) {
    Complex c = new Complex();

    c.x = a.x*b.x - a.y*b.y;
    c.y = a.x*b.y + a.y*b.x;

    return c;
  }

  /**
   * Divides one input complex number by another.
   *
   * @param a the numerator
   * @param b the denominator
   * @return the ratio
   */
  public static Complex divide(Complex a, Complex b) {
    Complex c = new Complex();

    c.x = (a.x*b.x + a.y*b.y) / (b.x*b.x + b.y*b.y);
    c.y = (-1*a.x*b.y + a.y*b.x) / (b.x*b.x + b.y*b.y);

    return c;
  }

  /**
   * Squares an input complex number.
   *
   * @param z the complex number to square
   * @return its square
   */
  public static Complex square(Complex z) {
    return multiply(z, z);
  }

  /**
   * Computes the square root of an input complex number
   *
   * @param z the number whose square root to compute
   * @returns the square root
   */
  public static Complex sqrt(Complex z) {
    double mod = modulus(z);
    // return zero if input was zero, since we cannot take the
    // argument of the complex number zero
    if (mod == 0D)
      return new Complex(0D, 0D);
    mod = Math.sqrt(mod);
    double arg= argument(z);
    //    if (arg == -Math.PI)
    //      arg = Math.PI/2;
    //    else
          arg /= 2D;
    return new Complex(mod*Math.cos(arg), mod*Math.sin(arg));
  }

  /**
   * Compute the argument of a complex number, in the range (-pi,pi].
   *
   * @param z the complex number to examine
   * @return its argument
   */
  public static double argument(Complex z) {
    return Math.atan2(z.y, z.x);
  }

  /**
   * Compute this's argument.
   *
   * @return this's argument
   */
  public double argument() { 
    return argument(this); 
  }

  /**
   * Computes the complex conjugate of the input complex number.
   *
   * @param z the complex number to conjugate
   * @return the conjugate
   */
  public static Complex conj(Complex z) {
    return new Complex(z.x, -z.y);
  }

  /**
   * Computes this's complex conjugate.
   *
   * @return this's conjugate
   */
  public Complex conj() { 
    return conj(this); 
  }

  /**
   * Computes the modulus of an input complex number.
   *
   * @param z the number whose modulus to compute
   * @return the modulus
   */
  public static double modulus(Complex z) {
    return Math.sqrt(z.x*z.x + z.y*z.y);
  }

  /**
   * Compute this's modulus.
   *
   * @return this's modulus
   */
  public double modulus() {
    return modulus(this);
  }

  /**
   * Compute the modulus squared of a complex number.
   *
   * @param z the number whose modulus squared to compute
   * @return the modulus square
   */
  public static double modsq(Complex z) {
    return z.x*z.x + z.y*z.y;
  }

  /**
   * Compute this's modulus squared.
   *
   * @return this's modulus squared
   */
  public double modsq() {
    return modsq(this);
  }

  /**
   * Convert an input string to a complex number.  Really should throw a
   * NumberFormatException if the string has an inappropriate form (which
   * would then have to be caught wherever this method is used!); instead,
   * returns the complex number i in such cases.  Valid format of a string
   * is two numbers separated by a comma, indicating the real and imaginary
   * parts, respectively, of the complex number.
   *
   * @param the string to convert
   * @return the complex number corresponding to that string
   */
  public static Complex stringToComplex(String s) {
    StringTokenizer st = new StringTokenizer(s, ",", false);
    double x = 0D, y = 0D;
    try {
      String ss = st.nextToken();
      String tt = st.nextToken();
      x = Double.valueOf(ss).doubleValue(); 
      y = Double.valueOf(tt).doubleValue(); 
    }
    catch(Exception e) {
      x = 0D;
      y = 1D;
    }

    return new Complex(x, y);
  }

  /**
   * Computes the square of the distance between two complex numbers.
   *
   * @param a one input complex number
   * @param b and another
   * @return the square of the distance
   */
  public static double distsq(Complex a, Complex b) {
    double xdiff = a.x - b.x;
    double ydiff = a.y - b.y;
    return xdiff*xdiff + ydiff*ydiff;
  }

  /**
   * The default invocation of DoubleToString has <b>DOUBLEDIGITS</b> of
   * precision.
   *
   * @param d the double to be stringified
   * @return the string version
   */
  public static String DoubleToString(double d) {
    return DoubleToString(d, DOUBLEDIGITS);
  }

  /**
   * Converts a double to a string with a certain number of digits of
   * precision.
   *
   * @param d   the double to convert
   * @param res the number of digits of precision to use
   */
  public static String DoubleToString(double d, int res) {
    String s = Double.toString(d);
    int k = s.indexOf("e");
    String epart;
    int useres = res;
    if (k >= 0) {
      // there is some scientific notation, strip it off
      epart = s.substring(k, s.length());
      s = s.substring(0, k);
      if (res > 2)
	useres = res - 2;
    }
    else
      epart = new String("");
    int i = s.indexOf(period);
    if (i >= 0) {
      // there is a decimal point, take only res digits after it
      int j = i + useres + 1;
      if (j < s.length())
	return s.substring(0, j)+epart;
      else
	return s+epart;
    }
    return s+epart;
  }

  /**
   * Convert this complex number to a string, with the real and imaginary
   * parts separated by a comma.
   *
   * @return the string version of this complex number
   */
  public String toString() {
    return new String(DoubleToString(this.x)+","+DoubleToString(this.y));
  }

  /**
   * Compute the (1,1) entry of the canonical matrix with given trace.
   * Answer will be the x such that x+x^{-1}=trace, in fact the one chosen
   * to have modulus greater than or equal to one.
   *
   * @param trace the trace of the matrix we are interested in
   * @return the (1,1) entry of the corresponding canonical matrix
   */
  public static Complex make_z11(Complex trace) {
    Complex a = subtract(square(trace), FOUR);
    a = sqrt(a);
    Complex b = add(trace, a);
    // add or subtract the discriminant to have modulus greater than one
    if ((b.x*b.x + b.y*b.y) > 4D)
      return divide(b, TWO);
    a = subtract(trace, a);
    return divide(a, TWO);
  }

  /**
   * Compute the maxpow corresponding to the current trace and previously
   * calculated (1,1) entry.  For purely hyerpbolic and parabolic elements,
   * always returns 1; for elliptic, returns half the order of the element,
   * up to value <b>MAXELLIPORD</b>; for loxodromic elements, the result is
   * valid for Ford domains, but is not known to be correct for Dirichlet
   * domains -- although we use it in that case as well.
   *
   * @param trace the current trace
   * @param z11   the (1,1) entry of the generator with that trace
   * @return the corresponding Ford domain maxpow
   */
  public static int make_maxpow(Complex trace, Complex z11) {
    int mp;
    if (trace.y == 0D) {
      // not purely loxodromic
      if (Math.abs(trace.x) >= 2D) // parabolic or purely hyperbolic
        return 1;
      // elliptic
      Complex z11powers = new Complex(z11);
      for (mp=2; mp<=MAXELLIPORD; mp++) {
	z11powers = multiply(z11powers, z11);
        if ((distsq(z11powers, ONE) < IDENTRESOL) || (distsq(z11powers, MINUS_ONE) < IDENTRESOL)) {
	  // we're within the allowed resolution of the identity,
	  // so this is the order
	  mp /= 2;
          break;
	}
      }
    }
    else {
      // purely loxodromic
      double numer = modulus(subtract(square(z11), ONE));
      double denom = modulus(z11);
      mp = ((int)Math.floor(1D+Math.log(1D+numer/(denom-1D))/Math.log(denom)));
    }
    // return no more than MAXPOWMAX
    if (mp > MAXPOWMAX)
      return MAXPOWMAX;
    else
      return mp;
  }

}