/* ComplexFunctions.cp */

#include "ComplexFunctions.h"
#include <math.h>

typedef struct
{
	double		Re;
	double		Im;
	
} Comp;

static void ToPolar(double& r, double& theta, const double xRe, const double xIm);
static void ToRect(double& xRe, double& xIm, const double r, const double theta);

static inline double dabs(double x)
{
	if(x<0)
		x = -x;
	
	return x;
	
}/* dabs */

static void Mul(Comp& z, const Comp& x, const Comp& y)
{
	double		ac, bd;
	
	if(x.Re!=0 && x.Im!=0 && y.Re!=0 && y.Im!=0)
	{
		ac = x.Re*y.Re;
		bd = x.Im*y.Im;
		
		// the order here is important since z might be the same as x or y
		z.Im = (x.Re + x.Im)*(y.Re + y.Im) - ac - bd;
		z.Re = ac - bd;
		return;
	}
	
	if(x.Re==0)
	{
		ac = -x.Im*y.Im;
		z.Im = x.Im*y.Re;
		z.Re = ac;
		return;
	}
	
	if(x.Im==0)
	{
		ac = x.Re*y.Re;
		z.Im = x.Re*y.Im;
		z.Re = ac;
		return;
	}
	
	if(y.Re==0)
	{
		ac = -x.Im*y.Im;
		z.Im = x.Re*y.Im;
		z.Re = ac;
		return;
	}
	
	if(y.Im==0)
	{
		ac = x.Re*y.Re;
		z.Im = x.Im*y.Re;
		z.Re = ac;
		return;
	}
	
}/* Mul */


static void Div(Comp& z, const Comp& x, const Comp& y)
{
	double		ac, bd, den;
	
	if(y.Re==0 && y.Im==0)
	{
		z.Re = NAN;
		z.Im = NAN;
		return;
	}
	
	if(x.Im==0 && y.Im==0)
	{
		z.Re = x.Re/y.Re;
		z.Im = 0;
		return;
	}
	
	if(x.Re==0 && y.Re==0)
	{
		z.Re = x.Im/y.Im;
		z.Im = 0;
		return;
	}
	
	if(x.Re==0 && y.Im==0)
	{
		z.Re = 0;
		z.Im = x.Im/y.Re;
		return;
	}
	
	if(x.Im==0 && y.Re==0)
	{
		z.Re = 0;
		z.Im = -x.Re/y.Im;
		return;
	}
	
	den = y.Re*y.Re + y.Im*y.Im;
	
	ac = x.Re*y.Re;
	bd = x.Im*y.Im;
	
	// the order here is important since z might be the same as x or y
	z.Im = ((x.Re + x.Im)*(y.Re - y.Im) - ac + bd)/den;
	z.Re = (ac + bd)/den;
	
}/* Div */


void myExp(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp 		x, z, currentTerm;
	int         n;
	
	x.Re = xRe;
	x.Im = xIm;
	
	zRe = z.Re = currentTerm.Re = 1;
	zIm = z.Im = currentTerm.Im = 0;
	
	for(n=1; ;n++)
	{
		Mul(currentTerm, x, currentTerm);
		currentTerm.Re = currentTerm.Re/n;
		currentTerm.Im = currentTerm.Im/n;
		
		z.Re = z.Re + currentTerm.Re;
		z.Im = z.Im + currentTerm.Im;
		
		if(z.Re==zRe && z.Im==zIm)
			break;
		
		zRe = z.Re;
		zIm = z.Im;
	}
	
}/* myExp */


void myLog(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		r, theta;
	
	if(xRe==0 && xIm==0)
	{
		zRe = NAN;
		zIm = NAN;
		return;
	}
	
	ToPolar(r, theta, xRe, xIm);
	
	zRe = log(r);
	zIm = theta;
	
}/* myLog */


/*
	might use sin(x+iy) = sin(x)*cosh(y) + icos(x)*sinh(y)

*/
void mySin(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp 		x, z, currentTerm, xSqr;
	long		n;
	
	z.Re = x.Re = currentTerm.Re = xRe;
	z.Im = x.Im = currentTerm.Im = xIm;
	
	Mul(xSqr, x, x); // here in case z=x
	
	zRe = z.Re;
	zIm = z.Im;
	
	for(n=2; ;n+=2)
	{
		Mul(currentTerm, xSqr, currentTerm);
		currentTerm.Re = -currentTerm.Re/(n*(n+1));
		currentTerm.Im = -currentTerm.Im/(n*(n+1));
		
		z.Re = z.Re + currentTerm.Re;
		z.Im = z.Im + currentTerm.Im;
		
		if(z.Re==zRe && z.Im==zIm)
			break;
		
		zRe = z.Re;
		zIm = z.Im;
	}
	
}/* mySin */


/*
	might use cos(x+iy) = cos(x)*cosh(y) - isin(x)*sinh(y)

*/
void myCos(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp 		x, z, currentTerm, xSqr;
	long		n;
	
	x.Re = xRe;
	x.Im = xIm;
	
	Mul(xSqr, x, x); // here in case z=x
	
	zRe = z.Re = currentTerm.Re = 1;
	zIm = z.Im = currentTerm.Im = 0;
	
	for(n=1; ;n+=2)
	{
		Mul(currentTerm, xSqr, currentTerm);
		currentTerm.Re = -currentTerm.Re/(n*(n+1));
		currentTerm.Im = -currentTerm.Im/(n*(n+1));
		
		z.Re = z.Re + currentTerm.Re;
		z.Im = z.Im + currentTerm.Im;
		
		if(z.Re==zRe && z.Im==zIm)
			break;
		
		zRe = z.Re;
		zIm = z.Im;
	}
	
}/* myCos */


void myTan(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp 		z, sinx, cosx;
	
	myCos(cosx.Re, cosx.Im, xRe, xIm);
	
	if(cosx.Re==0 && cosx.Im==0)
	{
		zRe = NAN;
		zIm = NAN;
		return;
	}
		
	mySin(sinx.Re, sinx.Im, xRe, xIm);
	
	Div(z, sinx, cosx);
	
	zRe = z.Re;
	zIm = z.Im;
	
}/* myTan */


static void ToPolar(double& r, double& theta, const double xRe, const double xIm)
{
	r = sqrt(xRe*xRe + xIm*xIm);
	theta = atan2(xIm, xRe);

}/* ToPolar */


static void ToRect(double& xRe, double& xIm, const double r, const double theta)
{
	xRe = r*cos(theta);
	xIm = r*sin(theta);

}/* ToRect */


void mySqrt(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		r, costheta;
	bool		isNegative;
	
	if(xRe==0 && xIm==0)
	{
		zRe = zIm = 0;
		return;
	}
	
	// in case z=x
	if(xIm<0)
		isNegative = true;
	else
		isNegative = false;
	
	r = sqrt(xRe*xRe + xIm*xIm);
	
	costheta = xRe/r;
	
	zRe = sqrt(r*(1 + costheta)/2);
	zIm = sqrt(r*(1 - costheta)/2);
	
	if(isNegative)
		zIm = -zIm;
	
}/* mySqrt */


void mySqr(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		temp;
	
	temp = xRe*xRe - xIm*xIm;
	
	zIm = 2*xRe*xIm;
	zRe = temp;
	
}/* mySqr */


void mySinh(double& z, const double& x)
{
	if(dabs(x)>=.1)
   {
   		double temp;
   		temp = x;
   		temp = exp(x);
   		z = (temp - 1/temp)/2;
   		return;
   }

   // now abs(x)<.1 so use power series
   double t2, s, s1, t;
   long i;

   s = 0;
   t = 1;
   t2 = x*x;

   for (i = 2; ; i+=2)
   {
      s1 = s + t;
      if (s == s1)
      	break;
      s = s1;
      t = t*t2;
      t = t/(i*(i+1));
   }
 
   z = x*s;
   
}/* mySinh */


void myCosh(double& z, const double& x)
{
	double		temp;
	
	temp = exp(x);
	z = (temp + 1/temp)/2;
	
}/* myCosh */


void myASinh(double& z, const double& x)
{
	
	if(dabs(x)>=.1)
  	{
  		double temp = log(dabs(x) + sqrt(x*x+1));
  		if(x<0)
			z = -temp;
		else
			z = temp;
		return;
  	}
	
	// abs(x)<.1 so use power series
	
   double t2, s, s1, t;
   long i;

   s = 0;
   t = 1;
   t2 = -x*x;

   for (i = 2; ; i+=2)
   {
      s1 = s + t;
      if (s == s1)
      	break;
      s = s1;
      t = t*t2;
      t = t*(i-1)*(i-1);
      t = t/(i*(i+1));
   }

   z = x*s;
   
}/* myASinh */


void myACosh(double& z, const double& x)
{
	if(x<1)
	{
		z = NAN;
		return;
	}
		
	z = log(x + sqrt(x*x-1));
	
}/* myACosh */


void myATanh(double& z, const double& x)
{
	if(dabs(x)>=1)
	{
		z = NAN;
		return;
	}
		
	if(x>=.1)
	{
		z = log((1+x)/(1-x))/2;
		return;
	}
	
	// now x<.1 so use power series
   double t2, t3, s, s1, t;
   long i;

   s = 0;
   t = 1;
   t2 = x*x;
   t3 = 1;

   for (i = 3; ; i+=2)
   {
      s1 = s + t;
      if (s == s1)
      	break;
      s = s1;
      t3 = t3*t2;
      t = t3/i;
   }
   
   z = x*s;
}/* myATanh */

// sinh(x) = -i sin(ix)
void mySinhC(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		yRe, yIm;
	
	mySin(yRe, yIm, -xIm, xRe);
	
	zRe = yIm;
	zIm = -yRe;
	
}/* mySinhC */


// cosh(x) = cos(ix)
void myCoshC(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	myCos(zRe, zIm, -xIm, xRe);
	
}/* myCoshC */


// tanh(x) = -i tan(ix)
void myTanhC(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		yRe, yIm;
	
	myTan(yRe, yIm, -xIm, xRe);
	
	zRe = yIm;
	zIm = -yRe;
	
}/* myTanhC */

// x = sin(z)
// z = arctan(x/sqrt(1-x^2))
void myASin(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp	x, z, temp;
	
	if(xRe==1 && xIm==0)
	{
		zRe = 2*atan(1);
		zIm = 0;
		return;
	}
	
	if(xRe==-1 && xIm==0)
	{
		zRe = -2*atan(1);
		zIm = 0;
		return;
	}
	
	if(xIm==0 && dabs(xRe)<=1)
	{
		zRe = asin(xRe);
		zIm = 0;
		return;
	}
	
	x.Re = xRe;
	x.Im = xIm;
	
	mySqr(temp.Re, temp.Im, xRe, xIm);
	
	temp.Re = 1 - temp.Re;
	temp.Im = -temp.Im;
	
	mySqrt(temp.Re, temp.Im, temp.Re, temp.Im);
	Div(z, x, temp);
	
	myATan(z.Re, z.Im, z.Re, z.Im);
	
	zRe = z.Re;
	zIm = z.Im;
	
}/* myASin */


// x = cos(z)
// z = arctan(sqrt(1-x^2)/x)
void myACos(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp	x, z, temp;
	
	if(xRe==0 && xIm==0)
	{
		zRe = 2*atan(1);
		zIm = 0;
		return;
	}
	
	if(xIm==0 && dabs(xRe)<=1)
	{
		zRe = acos(xRe);
		zIm = 0;
		return;
	}
	
	x.Re = xRe;
	x.Im = xIm;
	
	mySqr(temp.Re, temp.Im, xRe, xIm);
	
	temp.Re = 1 - temp.Re;
	temp.Im = -temp.Im;
	
	mySqrt(temp.Re, temp.Im, temp.Re, temp.Im);
	Div(z, temp, x);
	
	myATan(z.Re, z.Im, z.Re, z.Im);
	
	zRe = z.Re;
	zIm = z.Im;
	
}/* myACos */


/*
	x = tan(z), want to calculate z
	Use:
	z = (-i/2)*log((1+ix)/(1-ix))
	singular at x = -i

*/
void myATan(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	Comp	temp1, temp2;
	double	yRe, yIm;
	
	if(xRe==0 && xIm==0)
	{
		zRe = 0;
		zIm = 0;
		return;
	}
	
	if(xRe==0 && xIm==-1)
	{
		zRe = NAN;
		zIm = NAN; //should be negativeinfinity
		return;
	}
	
	temp1.Re = 1 - xIm;
	temp1.Im = xRe;
	temp2.Re = 1 + xIm;
	temp2.Im = -xRe;
	
	Div(temp1, temp1, temp2);
	
	myLog(yRe, yIm, temp1.Re, temp1.Im);
	
	zRe = yIm/2;
	zIm = -yRe/2;
	if(xIm==0)
		zIm = 0;
	
}/* myATan */



// asinh(x) = -i*asin(ix)
void myASinhC(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		yRe, yIm;
	
	myASin(yRe, yIm, -xIm, xRe);
	
	zRe = yIm;
	zIm = -yRe;
	
}/* myASinhC */


// acosh(x) = -i*acos(x)
void myACoshC(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		yRe, yIm;
	
	myACos(yRe, yIm, xRe, xIm);
	
	zRe = yIm;
	zIm = -yRe;
	
}/* myACoshC */


// atanh(x) = -i*atan(ix)
void myATanhC(double& zRe, double& zIm, const double& xRe, const double& xIm)
{
	double		yRe, yIm;
	
	myATan(yRe, yIm, -xIm, xRe);
	
	zRe = yIm;
	zIm = -yRe;
	
}/* myATanhC */


void myPowCC(double& zRe, double& zIm,  const double& xRe, const double& xIm,
					const double& yRe, const double& yIm)
{
	Comp	y, temp;
	
	if(xRe==0 && xIm==0)
	{
		if(yRe<=0)
		{
			zRe = NAN;
			zIm = NAN;
		}
		else
			zRe = zIm = 0;
			
		return;
	}
		
	y.Re = yRe;
	y.Im = yIm;
		
	myLog(temp.Re, temp.Im, xRe, xIm);
	Mul(temp, y, temp);
	
	myExp(zRe, zIm, temp.Re, temp.Im);
	
}/* myPowCC */


void myPowDC(double& zRe, double& zIm,  const double& x,
					const double& yRe, const double& yIm)
{
	Comp	y, temp;
	double	tempD;
	
	if(x==0)
	{	
		if(yRe<=0)
		{
			zRe = NAN;
			zIm = NAN;
		}
		else
			zRe = zIm = 0;
			
		return;
	}
		
	if(x>0)
	{
		tempD = log(x);
		temp.Re = yRe*tempD;
		temp.Im = yIm*tempD;
	}
	else
	{
		temp.Re = log(-x);
		temp.Im = 4*atan(1);
		y.Re = yRe;
		y.Im = yIm;
		Mul(temp, y, temp);
	}
	
	myExp(zRe, zIm, temp.Re, temp.Im);
	
}/* myPowDC */


void myPowCD(double& zRe, double& zIm,  const double& xRe, const double& xIm,
					const double& y)
{
	Comp	temp;
	
	if(xRe==0 && xIm==0)
	{
		if(y<=0)
		{
			zRe = NAN;
			zIm = NAN;
		}
		else
			zRe = zIm = 0;
			
		return;
	}
		
	myLog(temp.Re, temp.Im, xRe, xIm);
	
	temp.Re = y*temp.Re;
	temp.Im = y*temp.Im;
	
	myExp(zRe, zIm, temp.Re, temp.Im);
	
}/* myPowCD */