/*
 *  fpComplexFuncs.cpp
 *  fp
 *
 *  Created by Robert Delaney on 11/8/09.
 *  Copyright 2009 Bob Delaney's Science Software. All rights reserved.
 *
 */

#include "fpComplex.h"


// calculates both sinh(x) and cosh(x) for real x using only one evaluation of exp(x)
// needed for calculation of sine and cosine for complex arguments
static void bothSandC(fp& S, fp& C, const fp& x)
{
	static fp		pointOne;
	fp				xt, St, temp, temp1, tempSave, t2, s, s1, t;
	static bool		initGood=false;
	INT32			i;
	
	if(!initGood)
	{
		init(pointOne, .1);
		initGood = true;
	}
	
	equate(xt, x);
	if(xt.i.n<0)
		xt.i.n = -xt.i.n;
	if(compare(xt, pointOne)==1)
	{
		exp(temp, x);
		tempSave = temp;
		div(temp1, 1., temp);
		sub(temp1, temp, temp1);
		div(St, temp1, 2.);
		goto coshCalc;
	}
	
	// now abs(x)<=.1 so use power series
	s.i.n = 0;
	s.e = 0;
	equate(t, 1.);
	mul(t2, x, x);
	
	for (i=2; ;i+=2)
	{
		add(s1, s, t);
		if(!compare(s, s1))
			break;
		equate(s, s1);
		mul(t, t, t2);
		div(t, t, i*(i+1));
	}
	mul(St, x, s);
	
coshCalc:
	
	if(compare(xt, pointOne)<1)
		exp(tempSave, x);
	
	equate(C, tempSave);
	div(tempSave, 1., tempSave);
	add(C, C, tempSave);
	div(C, C, 2.);
	
	S = St;
	
}/* bothSandC */


static void toPolar(fp& r, fp& theta, const fpComplex& x)
{
	r = abs(x);
	atan2(theta, x);
	
}/* toPolar */


static void ToRect(fpComplex& z, const fp& r, const fp& theta)
{
	z.re = r*cos(theta);
	z.im = r*sin(theta);
	
}/* ToRect */



fp abs(const fpComplex& x)
{
	fp		z;
	
	z = sqrt(x.re*x.re + x.im*x.im);
	
	return z;
	
}/* abs */


fpComplex conj(const fpComplex& x)
{
	fpComplex		z;
	
	z.re = x.re;
	z.im = -x.im;
	
	return z;
	
}/* conj */


void exp(fpComplex& z, const fpComplex& x)
{
	fp			temp;
	fpComplex	zt;
	
	exp(temp, x.re);
	zt.re = temp*cos(x.im);
	zt.im = temp*sin(x.im);
	
	z = zt;

}/* exp */


fpComplex exp(const fpComplex& x)
{
	fpComplex		z;
	
	exp(z, x);
	
	return z;
	
}/* exp */


void log(fpComplex& z, const fpComplex& x)
{
	fp		r, theta;
	
	if(!x)
	{
		z = 0.; // really undetermined
		return;
	}
	
	if(!x.im.i.n && x.re.i.n>0)
	{
		z.re = log(x.re);
		z.im = 0.;
		return;
	}
	 
	
	toPolar(r, theta, x);
	
	z.re = log(r);
	z.im = theta;
	
}/* log */


fpComplex log(const fpComplex& x)
{
	fpComplex		z;
	
	log(z, x);
	
	return z;
	
}/* log */


// sin(a + ib) = sin(a) cosh(b) + i cos(a) sinh(b)
void sin(fpComplex& z, const fpComplex& x)
{
	fp			S, C;
	fpComplex	xt;
	
	xt = x;
	
	bothSandC(S, C, xt.im);
	
	z.re = sin(xt.re)*C;
	z.im = cos(xt.re)*S;
	
}/* sin */


fpComplex sin(const fpComplex& x)
{
	fpComplex		z;
	
	sin(z, x);
	
	return z;
	
}/* sin */


// cos(a + ib) = cos(a) cosh(b) - i sin(a) sinh(b)
void cos(fpComplex& z, const fpComplex& x)
{
	fp			S, C;
	fpComplex	xt;
	
	xt = x;
	
	bothSandC(S, C, xt.im);
	
	z.re = cos(xt.re)*C;
	z.im = -sin(xt.re)*S;
	
}/* cos */


fpComplex cos(const fpComplex& x)
{
	fpComplex		z;
	
	cos(z, x);
	
	return z;
	
}/* cos */


void tan(fpComplex& z, const fpComplex& x)
{
	fpComplex		sinx, cosx;
	
	cos(cosx, x);
	
	if(!cosx)
	{
		z = 0.;
		return;
	}
	
	sin(sinx, x);
	
	z = sinx/cosx;
	
}/* tan */


fpComplex tan(const fpComplex& x)
{
	fpComplex		z;
	
	tan(z, x);
	
	return z;
	
}/* tan */


void atan2(fp& theta, const fpComplex& x)
{
	fp			temp, myTan, myPi;
	
	if(!x)
	{
		theta = 0.; // actually undetermined
		return;
	}
	
	myPi = Pi();
	
	if(!x.re.i.n)
	{
		theta = myPi/2;
		if(x.im.i.n<0)
			theta = -theta;
		return;
	}
	
	if(!x.im.i.n)
	{
		if(x.re.i.n>0)
			theta = 0.;
		else
			theta = myPi;
		return;
	}
	
	myTan = abs(x.im/x.re);
	
	atan(theta, myTan);
	
	if(x.re.i.n<0 && x.im.i.n>=0)
		theta = myPi - theta;
	
	if(x.re.i.n<0 && x.im.i.n<0)
		theta = theta - myPi;
	
	if(x.re.i.n>=0 && x.im.i.n<0)
		theta = -theta;
	
}/* atan2 */


fp atan2(const fpComplex& x)
{
	fp			z;
	
	atan2(z, x);
	
	return z;
	
}/* atan2 */


void sqrt(fpComplex& z, const fpComplex& x)
{
	fpComplex		xt;
	fp				r, costheta;
	
	if(!x)
	{
		z = 0.;
		return;
	}
	
	xt = x; // in case z is x
	
	r = abs(x);
	costheta = x.re/r;
	
	z.re = sqrt(r*(1 + costheta)/2);
	z.im = sqrt(r*(1 - costheta)/2);
	
	if((xt.re.i.n<0 && xt.im.i.n<0) || (xt.re.i.n>=0 && xt.im.i.n<0))
		z.im = -z.im;
	
	if(xt.re>=0 && xt.im==0)
		z.im = 0.;
	
	if(xt.re<=0 && xt.im==0)
		z.re = 0.;
	
}/* sqrt */


fpComplex sqrt(const fpComplex& x)
{
	fpComplex		z;
	
	sqrt(z, x);
	
	return z;
	
}/* sqrt */


void sqr(fpComplex& z, const fpComplex& x)
{
	z = x*x;
	
}/* sqr */


fpComplex sqr(const fpComplex& x)
{
	fpComplex		z;
	
	sqr(z, x);
	
	return z;
	
}/* sqr */


// sinh(x) = -i sin(ix)
void sinh(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt;
	
	zt.re = -x.im;
	zt.im = x.re;
	
	sin(zt, zt);
	
	z.re = zt.im;
	z.im = -zt.re;
	
}/* sinh */


fpComplex sinh(const fpComplex& x)
{
	fpComplex		z;
	
	sinh(z, x);
	
	return z;
	
}/* sinh */


// cosh(x) = cos(ix)
void cosh(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt;
	
	zt.re = -x.im;
	zt.im = x.re;
	
	cos(zt, zt);
	
	z = zt;
	
}/* cosh */


fpComplex cosh(const fpComplex& x)
{
	fpComplex		z;
	
	cosh(z, x);
	
	return z;
	
}/* cosh */


// tanh(x) = -i tan(ix)
void tanh(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt;
	
	zt.re = -x.im;
	zt.im = x.re;
	
	tan(zt, zt);
	
	z.re = zt.im;
	z.im = -zt.re;
	
}/* tanh */


fpComplex tanh(const fpComplex& x)
{
	fpComplex		z;
	
	tanh(z, x);
	
	return z;
	
}/* tanh */


// x = sin(z)
// z = arctan(x/sqrt(1-x^2))
// we have to worry about the sign of the sqrt (cos(z))
void asin(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt, temp;
	fp				myPi, one=1.;
	
	if(!x)
	{
		z = 0.;
		return;
	}
	
	myPi = Pi();
	
	if(x.re==one && !x.im.i.n)
	{
		z.re = myPi/2;
		z.im = 0.;
		return;
	}
	
	if(x.re==-one && !x.im.i.n)
	{
		z.re = -myPi/2;
		z.im = 0.;
		return;
	}
	
	sqr(temp, x);
	temp.re = one - temp.re;
	temp.im = -temp.im;
	
	sqrt(temp, temp);
	div(zt, x, temp);
	
	atan(z, zt);
	
}/* asin */


fpComplex asin(const fpComplex& x)
{
	fpComplex		z;
	
	asin(z, x);
	
	return z;
	
}/* asin */


// x = cos(z)
// z = arctan(sqrt(1-x^2)/x)
// we have to worry about the sign of the sqrt (sin(z))
void acos(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt, temp;
	fp				myPi, one=1.;
	
	if(x.re==one && !x.im.i.n)
	{
		z = 0.;
		return;
	}
	
	myPi = Pi();
	if(x.re==-one && !x.im.i.n)
	{
		z.re = myPi;
		z.im = 0.;
		return;
	}
	
	if(!x)
	{
		z.re = myPi/2;
		z.im = 0.;
		return;
	}
	
	sqr(temp, x);
	temp.re = one - temp.re;
	temp.im = -temp.im;
	
	sqrt(temp, temp);
	div(zt, temp, x);
	
	atan(z, zt);
	
}/* acos */


fpComplex acos(const fpComplex& x)
{
	fpComplex		z;
	
	acos(z, x);
	
	return z;
	
}/* acos */


/*
 x = tan(z), want to calculate z
 Use:
 z = (-i/2)*log((1+ix)/(1-ix))
 singular at x = -i
 
 */
void atan(fpComplex& z, const fpComplex& x)
{
	fpComplex		xt, zt, temp1, temp2;
	fp				one=1.;
	
	if(!x)
	{
		z = 0.;
		return;
	}
	
	if(!x.re.i.n && x.im==-one)
	{
		z = 0.; // should be negative infinity
		return;
	}
	
	temp1.re = one - x.im;
	temp1.im = x.re;
	
	temp2.re = one + x.im;
	temp2.im = -x.re;
	
	div(temp1, temp1, temp2);
	
	log(zt, temp1);
	
	xt = x;
	
	z.re = zt.im/2;
	z.im = -zt.re/2;
	
	if(!xt.im.i.n)
		z.im = 0.;

}/* atan */


fpComplex atan(const fpComplex& x)
{
	fpComplex		z;
	
	atan(z, x);
	
	return z;
	
}/* atan */


// asinh(x) = -i*asin(ix)
void asinh(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt;
	
	zt.re = -x.im;
	zt.im = x.re;
	
	asin(zt, zt);
	
	z.re = zt.im;
	z.im = -zt.re;
	
}/* asinh */


fpComplex asinh(const fpComplex& x)
{
	fpComplex		z;
	
	asinh(z, x);
	
	return z;
	
}/* asinh */


// acosh(x) = i*acos(x)
void acosh(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt;
	
	zt = x;
	
	acos(zt, zt);
	
	z.re = -zt.im;
	z.im = zt.re;
	
}/* acosh */


fpComplex acosh(const fpComplex& x)
{
	fpComplex		z;
	
	acosh(z, x);
	
	return z;
	
}/* acosh */


// atanh(x) = -i*atan(ix)
void atanh(fpComplex& z, const fpComplex& x)
{
	fpComplex		zt;
	
	zt.re = -x.im;
	zt.im = x.re;
	
	atan(zt, zt);
	
	z.re = zt.im;
	z.im = -zt.re;
	
}/* atanh */


fpComplex atanh(const fpComplex& x)
{
	fpComplex		z;
	
	atanh(z, x);
	
	return z;
	
}/* atanh */



// z = x^y
// we have a problem when xRe<0 and xIm going through zero
// this causes a discontinuous change in theta when we take the log
// so we should test xRe and adjust as necessary
// but this just moves the discontinuity elsewhere,
// so I do it this way to be consistent with my Extended Plugin
// 4/14/2010: I changed my mind.
void pow(fpComplex& z, const fpComplex& x, const fpComplex& y)
{
	fpComplex		temp;
	fp				myPi;
	
	if(!y)
	{
		z = 1.;  // so 0^0 = 1
		return;
	}
	
	if(!x)
	{
		z = 0.; // actually undefined if y.re<0
		return;
	}
	
	if(x==1)
	{
		z = 1.;
		return;
	}
	
	myPi = Pi();
	log(temp, x);
	
	/*
	if(x.re.i.n<0)
	{
		if(x.im.i.n<0)
			temp.im = temp.im + 2*myPi;
	}
	 */
	
	mul(temp, y, temp);
	
	exp(z, temp);
	
}/* pow */


fpComplex pow(const fpComplex& x, const fpComplex& y)
{
	fpComplex		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


void pow(fpComplex& z, const fpComplex& x, const fp& y)
{
	fpComplex		yt;
	
	yt = y;
	
	pow(z, x, yt);
	
}/* pow */


fpComplex pow(const fpComplex& x, const fp& y)
{
	fpComplex		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


void pow(fpComplex& z, const fpComplex& x, double y)
{
	fpComplex		yt;
	INT32			yInt;
	double			yD;
	
	yInt = y;
	yD = yInt;
	if(y==yD)
	{
		pow(z, x, yInt);
		return;
	}
	
	yt = y;
	pow(z, x, yt);
	
}/* pow */


fpComplex pow(const fpComplex& x, double y)
{
	fpComplex		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


void pow(fpComplex& z, const fpComplex& x, INT32 y)
{
	INT32				i, n, sy;
	static UINT32		mask; // 0x80000000
	static fpComplex	one;
	fpComplex			zt;
	static bool			initGood=false;
	bool				isNegative;
	
	if(!initGood)
	{
		one = 1.;
		mask = 1;
		mask = (mask<<31); // 0x80000000
		initGood = true;
	}
	
	if(x==one)
	{
		z = one;
		return;
	}
	
	if(!y)
	{
		z = one;
		return;
	}
	
	if(!x)
	{
		z = 0.;
		return;
	}
	
	if(y>0)
	{
		isNegative = false;
		sy = y;
	}
	else
	{
		isNegative = true;
		sy = -y;
	}
	
	zt = one;
	n = NumBits(sy);
	
	sy = (sy<<(32-n));  // go to leading 1 bit
	
	for(i=0;i<n;i++)
	{
		mul(zt, zt, zt); // zt^2
		if(sy & mask)
			mul(zt, zt, x);
		sy = (sy<<1);
	}
	
	z = zt;
	if(isNegative)
		z = 1/z;
	
}/* pow */


fpComplex pow(const fpComplex& x, INT32 y)
{
	fpComplex		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


void pow(fpComplex& z, const fp& x, const fpComplex& y)
{
	fpComplex		xt;
	
	xt = x;
	pow(z, xt, y);
	
}/* pow */


fpComplex pow(const fp& x, const fpComplex& y)
{
	fpComplex		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


void pow(fpComplex& z, double x, const fpComplex& y)
{
	fpComplex		xt;
	
	xt = x;
	pow(z, xt, y);
	
}/* pow */


fpComplex pow(double x, const fpComplex& y)
{
	fpComplex		z;
	
	pow(z, x, y);
	
	return z;
	
}/* pow */


