import java.awt.*;

public class Complex {

  /*
	* Paul Garrett, 18 July 98
	*
	* basic Complex number arithmetic, linear fractional transformations, etc.
	*
	* with static + corresponding non-static methods
	*
	* (NB: This does not match ~/w/a02/Complex.java ... )
	*/

  double x, y;

  final static Complex I = new Complex(0D, 1D); // convenient constant

  public Complex(double x, double y) {
	 this.x = x;
	 this.y = y;
  }

  public Complex(double x) {
	 this(x,0D);
  }

  /**********************************/

  public final static synchronized double re(Complex alpha) {
	 return alpha.x;
  }

  public final synchronized double re() {
	 return Complex.re(this);
  }

  /********************/

  public static final synchronized double im(Complex alpha) {
	 return alpha.y;
  }

  public final synchronized double im() {
	 return Complex.im(this);
  }

  /********************/
  
  public static final synchronized Complex add(Complex alpha, Complex beta) {
	 return new Complex(alpha.re() + beta.re(), alpha.im() + beta.im());
  }

  public final synchronized Complex add(Complex beta) {
	 return Complex.add(this, beta);
  }

  /*********************/

  public static final synchronized Complex sub(Complex alpha, Complex beta) {
	 return new Complex(alpha.re() - beta.re(), alpha.im() - beta.im());
  }

  public final synchronized Complex sub(Complex beta) {
	 return Complex.sub(this, beta);
  }

  /*********************/
    
  public static final synchronized Complex mul(Complex alpha, Complex beta) {
	 return new Complex(alpha.re() * beta.re() - alpha.im() * beta.im(),
							  alpha.re() * beta.im() + alpha.im() * beta.re());
  }

  public final synchronized Complex mul(Complex beta) {
	 return Complex.mul(this, beta);
  }

  /**********************/
  
  public static final synchronized Complex realMul(double r, Complex alpha) {
	 return new Complex(r * alpha.x, r * alpha.y);
  }

  public final synchronized Complex realMul(double r) {
	 return Complex.realMul(r, this);
  }
  
  /************************/

  public static final synchronized Complex realDiv(double r, Complex alpha) {
	 return new Complex(alpha.re()/r, alpha.im()/r);
  }
  
  public final synchronized Complex realDiv(double r) {
	 return Complex.realDiv(r, this);
  }

  /****************************/
  
  public static final synchronized Complex conj(Complex alpha) {
	 return new Complex(alpha.re(), -alpha.im());
  }
  
  public final synchronized Complex conj() {
	 return Complex.conj(this);
  }

  /**************************/

  public static final synchronized double normSq(Complex alpha) {
	 return alpha.x * alpha.x + alpha.y * alpha.y;
  }

  public final synchronized double normSq() {
	 return Complex.normSq(this);
  }

  /*********************/

  public static final synchronized double norm(Complex alpha) {
	 return Math.sqrt(alpha.x*alpha.x + alpha.y*alpha.y);
  }

  public final synchronized double norm() {
	 return Complex.norm(this);
  }

  /********************/

  public static final synchronized Complex inv(Complex alpha) {
	 double x = alpha.x;
	 double y = alpha.y;
	 double nsq = x*x+y*y;
	 if (nsq <= .0000001) { //__EDIT__??
		nsq = .0000001;
	 }
	 return new Complex(x/nsq, -y/nsq);
  }

  public final synchronized Complex inv() {
	 return Complex.inv(this);
  }

  /**********************/

  public static final synchronized Complex div(Complex alpha, Complex beta) {
	 return alpha.mul(beta.inv());
  }

  public final synchronized Complex div(Complex beta) {
	 return Complex.div(this, beta);
  }

  /************************/

  public static final synchronized Complex neg(Complex alpha) {
	 return new Complex(-alpha.re(), -alpha.im());
  }

  public final synchronized Complex neg() {
	 return Complex.neg(this);
  }

  /*********************/

  /*
	* This returns and "arg" between -PI and +PI
	*
	* NB arcTangent is "stabilized" against division by 0
	*/
  public static final synchronized double arg(Complex alpha) {
	 return arcTangent2(alpha.y, alpha.x);
  }

  public final synchronized double arg() {
	 return arcTangent2(this.y, this.x); // see below: custom atan2()
  }

  /*****************************
	*
	* Apply a linear fractional transformation to complex number
	*/

  public static final synchronized Complex LFT(Complex[][] g, Complex z) {
	 return g[0][0].mul(z).add(g[0][1]).div(g[1][0].mul(z).add(g[1][1]));
	 // just az+b / cz+d form...
  }

  public final synchronized Complex LFT(Complex[][] g) {
	 return Complex.LFT(g, this);
  }

  /**********************************************
	* 
	* For command-line output, such as diagnostics
	*/ 

  public static final synchronized void print(Complex z) {
	 System.out.print("(" + z.re() + ") + i(" + z.im() + ")");
  }

  public final synchronized void print() {
	 System.out.print("(" + this.re() + ") + i(" + this.im() + ")");
  }

  public static final synchronized void println(Complex z) {
	 System.out.println("(" + z.re() + ") + i(" + z.im() + ")");
  }

  public final synchronized void println() {
	 System.out.println("(" + this.re() + ") + i(" + this.im() + ")");
  }

  /****************************	 
	*
	* Some compilers, such as VisualAge1.1.2, do not have a good atan2()
	*/
  
  public static synchronized final double arcTangent2(double y, double x) {
	 if (Math.abs(x) < .0001 * Math.abs(y)) { // effectively 0!?
		if (y > .0001) {
		  return Math.PI/2;
		}
		else {
		  return -Math.PI/2;
		}
	 }
	 else { // case that x is not too tiny by comparison to y
		double angle = Math.atan(Math.abs(y)/Math.abs(x)); // not Math.atan2()
		if (y<0 && x>=0) {
		  angle = -angle;
		}
		else if (y>=0 && x<0) {
		  angle = Math.PI - angle;
		}
		else if (y<=0 && x<0) {
		  angle = angle - Math.PI;
		}
		return angle;
	 }
  }

  public final static synchronized int radToDeg(double rad) {
	 return ((((int) Math.round(180 * rad / Math.PI)) % 360) + 360) % 360;
  }

  public final static synchronized Complex pointToComplex(Point p, double scale) {
	 double x,y;
	 x = (((double) p.x) / scale) - 1;
	 y = (((double) p.y) / scale) - 1;
	 return new Complex(x,y); // in unit disk
  }

  public final static synchronized Point complexToPoint(Complex z, double scale) {
	 int x,y;
	 x = (int) Math.round(scale + scale * z.re());
	 y = (int) Math.round(scale + scale * z.im());
	 return new Point(x,y);
  }

  // non-static version of previous:
  public final Point complexToPoint(double scale) {
	 return complexToPoint(this, scale);
  }

  /*
	* The following maps a screen-coord point to a complex number in the 
	* unit disk, acts upon it by the linear fractional transformation mx,
	* and then draws it in screen coordinates, inside 
	* Rectangle(0,0,2*scale,2*scale) with radius "scale"
	*/

  public final static synchronized void drawGeodesic(Graphics g, Point[] pp,
																	  Complex[][] mx, double scale) {
	 // "scale" is in screen dimensions
	 // everything maps through "mx" as LFT
	 Complex z = Complex.LFT(mx, Complex.pointToComplex(pp[0], scale));
	 Complex w = Complex.LFT(mx, Complex.pointToComplex(pp[1], scale));

	 Complex cen = getCenter(z,w);
	 double rad = Complex.norm(z.sub(cen));

	 double[] angles = getAngles(z.sub(cen), w.sub(cen));

	 drawArc(g, cen, rad, angles[0], angles[1], scale); // last is "global"...
  }

  /*
	* The following embodies the main computational trick,
	* computing the center point for the circle passing through
	* 2 given points and orthogonal to the unit circle in the complex
	* plane
	*/

  public static final synchronized Complex getCenter(Complex z, Complex w) {
	 Complex num = z.realMul(1+w.normSq()).sub(w.realMul(1+z.normSq()));
	 Complex denom = z.mul(w.conj()).sub(z.conj().mul(w));
	 return num.div(denom);
  }

  /*
	* Given mapping matrix "mx" and "newPt" and "oldPt" (and "scale")
	* compute new mapping matrix
	*
	* This should replace the entries of "old" with their images...
	*/ 

  public final static synchronized void updateMapMx(Complex[][] mx, Point basePt, 
																	  Point newPt, double scale) {
	 Complex z = new Complex(((double) newPt.x - basePt.x)/scale, 
										 ((double) newPt.y - basePt.y)/scale);
	 Complex a = mx[0][0];
	 Complex b = mx[0][1];
	 Complex c = mx[1][0];
	 Complex d = mx[1][1];

	 mx[0][0] = a.add(z.mul(c));
	 mx[0][1] = b.add(z.mul(d));
	 mx[1][0] = c.add(z.conj().mul(a));
	 mx[1][1] = d.add(z.conj().mul(b));
  }

  /*
	* Following reduces angle to range 0 to 2*Math.PI
	*
	*/

  public static final synchronized double reduceAngle(double a) {
	 while (a > 2 * Math.PI) {
		a -= 2 * Math.PI;
	 }
	 while (a < 0) {
		a += 2*Math.PI;
	 }
	 return a;
  }


  /************************************************************8
	*
	* Should determine startAngle and arcAngle (in degrees)
	*
	* To draw piece of circle passing through 2 given points, _not_
	* passing through 1/z.conj()
	*
	* So these are two points in the unit circle, and we plan to draw an
	* arc from one to the other...
	*
	*/

  public static final synchronized double[] getAngles(Complex z, Complex w) {
	 double za = z.arg();
	 double wa = w.arg();
	 double startAngle, arcAngle;
	 if (za >= wa) {
		if (za - wa <= Math.PI) {
		  startAngle = reduceAngle(wa);
		  arcAngle = reduceAngle(za - wa);
		}
		else { // za - wa > Math.PI
		  startAngle = reduceAngle(za);
		  arcAngle = reduceAngle(wa - za);
		}
	 }
	 else { // wa > za
		if (wa - za > Math.PI) {
		  startAngle = reduceAngle(wa);
		  arcAngle = reduceAngle(za - wa);
		}
		else {
		  startAngle = reduceAngle(za);
		  arcAngle = reduceAngle(wa - za);
		}
	 }
	 double[] angles = {startAngle, arcAngle};
	 return angles;
  }

  /*
	* Standard complex exponentiation, in polar coordinates
	*/

  public static final synchronized Complex cxExp(double r, double theta) {
	 return new Complex(r*Math.cos(theta), r*Math.sin(theta));
  }

  /*
	* Avoid the whole-degrees inaccuracy in Graphics.drawArc()
	*/ 

  public static final synchronized void drawArc(Graphics g, Complex center, double r, 
																double startA, double arcAngle,
																double scale) { // latter "global"
	 double inc = arcAngle/10;
	 for (int i=0; i<10; i++) { // break arc into 10 pieces!? more!?
		Point pt1 = complexToPoint(center.add(cxExp(r, startA + i * inc)), scale);
		Point pt2 = complexToPoint(center.add(cxExp(r, startA + (i+1) * inc)), scale);
		g.drawLine(pt1.x, pt1.y, pt2.x, pt2.y);
	 }
  }


  /************************************/


} // End of Complex class

/**********************************************************************/