/* ExtComplex.cp */

#include "myTypedef.h"
#include "fpMath.h"
#include "fpConv.h"
#include "fpFuncs.h"
#include "ExtComplex.h"
#include <stdlib.h>

bool init(ExtComplex& x)
{
	init(x.Re);
	init(x.Im);
	return true;
	
}/* init */

bool isEqualZero(ExtComplex& z)
{
	if(z.Re.i.n==0 && z.Im.i.n==0)
		return true;
	else
		return false;
	
}/* isEqualZero */


bool isEqual(ExtComplex& z, ExtComplex& x)
{
	if(!compare(z.Re, x.Re) && !compare(z.Im, x.Im))
		return true;
	else
		return false;
		
}/* isEqual */


void assignZeroC(ExtComplex& z)
{
	//z.Re.i.n = 0;
	//z.Im.i.n = 0;
	equate(z.Re, 0.);
	equate(z.Im, 0.);
	
}/* assignZeroC */


void EquateC_C(ExtComplex& z, const ExtComplex& x)
{
	equate(z.Re, x.Re);
	equate(z.Im, x.Im);
	
}/* EquateC_C */


void AbsC(fp& myAbs, ExtComplex& z)
{
	static fp		s;
	static bool		arrayGood=false;
	
	if(!arrayGood)
	{
		init(s);
		arrayGood = true;
	}
	
	mul(myAbs, z.Re, z.Re);
	mul(s, z.Im, z.Im);
	add(myAbs, myAbs, s);
	sqrt(myAbs, myAbs);
	
}/* AbsC */


void Add(ExtComplex& z, const ExtComplex& x, const ExtComplex& y)
{
	add(z.Re, x.Re, y.Re);
	add(z.Im, x.Im, y.Im);
	
}/* Add */


void Sub(ExtComplex& z, const ExtComplex& x, const ExtComplex& y)
{
	sub(z.Re, x.Re, y.Re);
	sub(z.Im, x.Im, y.Im);
	
}/* Sub */


void Mul(ExtComplex& z, const ExtComplex& x, const ExtComplex& y)
{
	static fp		ac, bd, s;
	static bool		arrayGood=false;
	
	if(!arrayGood)
	{
		init(ac);
		init(bd);
		init(s);
		arrayGood = true;
	}
	
	if(x.Im.i.n==0 && y.Im.i.n==0)
	{
		mul(z.Re, x.Re, y.Re);
		z.Im.i.n = 0;
		return;
	}
	
	if(x.Re.i.n!=0 && x.Im.i.n!=0 && y.Re.i.n!=0 && y.Im.i.n!=0)
	{
		
		mul(ac, x.Re, y.Re);
		
		mul(bd, x.Im, y.Im);
		
		// the order here is important since z might be the same as x or y
		
		add(z.Im, x.Re, x.Im);
		add(s, y.Re, y.Im);
		mul(z.Im, z.Im, s);
		add(s, ac, bd);
		sub(z.Im, z.Im, s);
		sub(z.Re, ac, bd);
		return;
	}
	
	if(x.Re.i.n==0)
	{
		mul(ac, x.Im, y.Im);
		ac.i.n = -ac.i.n;
		mul(z.Im, x.Im, y.Re);
		equate(z.Re, ac);
		return;
	}
	
	if(x.Im.i.n==0)
	{
		mul(ac, x.Re, y.Re);
		mul(z.Im, x.Re, y.Im);
		equate(z.Re, ac);
		return;
	}
	
	if(y.Re.i.n==0)
	{
		mul(ac, x.Im, y.Im);
		ac.i.n = -ac.i.n;
		mul(z.Im, x.Re, y.Im);
		equate(z.Re, ac);
		return;
	}
	
	if(y.Im.i.n==0)
	{
		mul(ac, x.Re, y.Re);
		mul(z.Im, x.Im, y.Re);
		equate(z.Re, ac);
		return;
	}
	
}/* Mul */


void Div(ExtComplex& z, const ExtComplex& x, const ExtComplex& y)
{
	static fp		ac, bd, den, s;
	static bool		arrayGood=false;
	
	if(!arrayGood)
	{
		init(ac);
		init(bd);
		init(den);
		init(s);
		arrayGood = true;
	}
	
	if(y.Re.i.n==0 && y.Im.i.n==0)
	{
		z.Re.i.n = 0;
		z.Im.i.n = 0;
		return;
	}
	
	if(x.Im.i.n==0 && y.Im.i.n==0)
	{
		div(z.Re, x.Re, y.Re);
		z.Im.i.n = 0;
		return;
	}
	
	if(x.Re.i.n==0 && y.Re.i.n==0)
	{
		div(z.Re, x.Im, y.Im);
		z.Im.i.n = 0;
		return;
	}
	
	if(x.Re.i.n==0 && y.Im.i.n==0)
	{
		z.Re.i.n = 0;
		div(z.Im, x.Im, y.Re);
		return;
	}
	
	if(x.Im.i.n==0 && y.Re.i.n==0)
	{
		z.Re.i.n = 0;
		div(z.Im, x.Re, y.Im);
		z.Im.i.n = -z.Im.i.n;
		return;
	}
	
	mul(den, y.Re, y.Re);
	mul(s, y.Im, y.Im);
	add(den, den, s);
	mul(ac, x.Re, y.Re);
	mul(bd, x.Im, y.Im);
	
	// the order here is important since z might be the same as x or y
	
	add(z.Im, x.Re, x.Im);
	sub(s, y.Re, y.Im);
	mul(z.Im, z.Im, s);
	sub(s, ac, bd);
	sub(z.Im, z.Im, s);
	div(z.Im, z.Im, den);
	add(z.Re, ac, bd);
	div(z.Re, z.Re, den);
	
}/* Div */


// z = sqrt(x)
// uses Newton's Method z = (z + x/z)/2
void mySqrtCC(ExtComplex& z, ExtComplex& x)
{
	static fp		r, costheta, temp1, temp2, one, two;
	static bool		arrayGood=false;
	bool			isNegative;
	
	if(!arrayGood)
	{
		init(r);
		init(costheta);
		init(temp1);
		init(temp2);
		one.i.n = one.i.nn = 1;
		one.i.b = (UINT32*)malloc(one.i.nn*sizeof(UINT32));
		one.i.b[0] = 1;
		one.e = 0;
		two.i.n = two.i.nn = 1;
		two.i.b = (UINT32*)malloc(two.i.nn*sizeof(UINT32));
		two.i.b[0] = 2;
		two.e = 0;
		
		arrayGood = true;
	}
	
	if(x.Re.i.n==0 && x.Im.i.n==0)
	{
		assignZeroC(z);
		return;
	}
	
	// in case z=x
	if(x.Im.i.n < 0)
	{
		isNegative = true;
		x.Im.i.n = -x.Im.i.n;
	}
	else
		isNegative = false;
	
	
	AbsC(r, x);
	
	div(costheta, x.Re, r);
	
	add(temp1, one, costheta);
	sub(temp2, one, costheta);
	
	div(r, r, two);
	
	mul(temp1, temp1, r);
	mul(temp2, temp2, r);
	
	sqrt(z.Re, temp1);
	sqrt(z.Im, temp2);
	
	if(isNegative)
		z.Im.i.n = -z.Im.i.n;

}/* mySqrtCC */

void InverseC(ExtComplex& z, ExtComplex& x)
{
	static fp		den, s;
	static bool		arrayGood=false;
	
	if(!arrayGood)
	{
		init(den);
		init(s);
		arrayGood = true;
	}
	
	if(x.Re.i.n==0 && x.Im.i.n==0)
		return;
	
	mul(den, x.Re, x.Re);
	mul(s, x.Im, x.Im);
	add(den, den, s);
	sqrt(den, den);
	div(z.Re, x.Re, den);
	div(z.Im, x.Im, den);
	z.Im.i.n = -z.Im.i.n;

}/* InverseC */