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

#include "fpComplex.h"

// z = x + y
void add(fpComplex& z, const fpComplex& x, const fpComplex& y)
{
	z.re = x.re + y.re;
	z.im = x.im + y.im;

}/* add */


// z = x + y
void add(fpComplex& z, const fpComplex& x, const fp& y)
{
	z.re = x.re + y;
	z.im = x.im;
	
}/* add */


// z = x + y
void add(fpComplex& z, const fpComplex& x, const mb& y)
{
	z.re = x.re + y;
	z.im = x.im;
	
}/* add */


// z = x + y
void add(fpComplex& z, const fpComplex& x, double y)
{
	z.re = x.re + y;
	z.im = x.im;
	
}/* add */


// z = x + y
void add(fpComplex& z, const fp& x, const fpComplex& y)
{
	z.re = x + y.re;
	z.im = y.im;
	
}/* add */


// z = x + y
void add(fpComplex& z, const mb& x, const fpComplex& y)
{
	z.re = x + y.re;
	z.im = y.im;
	
}/* add */


// z = x + y
void add(fpComplex& z, double x, const fpComplex& y)
{
	z.re = x + y.re;
	z.im = y.im;
	
}/* add */


// z = x - y
void sub(fpComplex& z, const fpComplex& x, const fpComplex& y)
{
	z.re = x.re - y.re;
	z.im = x.im - y.im;
	
}/* sub */


// z = x - y
void sub(fpComplex& z, const fpComplex& x, const fp& y)
{
	z.re = x.re - y;
	z.im = x.im;
	
}/* sub */


// z = x - y
void sub(fpComplex& z, const fpComplex& x, const mb& y)
{
	z.re = x.re - y;
	z.im = x.im;
	
}/* sub */


// z = x - y
void sub(fpComplex& z, const fpComplex& x, double y)
{
	z.re = x.re - y;
	z.im = x.im;
	
}/* sub */


// z = x - y
void sub(fpComplex& z, const fp& x, const fpComplex& y)
{
	z.re = x - y.re;
	z.im = -y.im;
	
}/* sub */


// z = x - y
void sub(fpComplex& z, const mb& x, const fpComplex& y)
{
	z.re = x - y.re;
	z.im = -y.im;
	
}/* sub */


// z = x - y
void sub(fpComplex& z, double x, const fpComplex& y)
{
	z.re = x - y.re;
	z.im = -y.im;
	
}/* sub */


// z = x * y
void mul(fpComplex& z, const fpComplex& x, const fpComplex& y)
{
	fp		ac, bd; // x=a+ib  y = c+id
	
	if(x.re.i.n && x.im.i.n && y.re.i.n && y.im.i.n)
	{
		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.i.n)
	{
		ac = -x.im*y.im;
		z.im = x.im*y.re;
		z.re = ac;
		return;
	}
	
	if(!x.im.i.n)
	{
		ac = x.re*y.re;
		z.im = x.re*y.im;
		z.re = ac;
		return;
	}
	
	if(!y.re.i.n)
	{
		ac = -x.im*y.im;
		z.im = x.re*y.im;
		z.re = ac;
		return;
	}
	
	if(!y.im.i.n)
	{
		ac = x.re*y.re;
		z.im = x.im*y.re;
		z.re = ac;
	}
	
}/* mul */


// z = x * y
void mul(fpComplex& z, const fpComplex& x, const fp& y)
{
	fpComplex		yC;
	
	yC.re = y;
	yC.im = 0;
	
	mul(z, x, yC);
	
}/* mul */


// z = x * y
void mul(fpComplex& z, const fpComplex& x, const mb& y)
{
	fpComplex		yC;
	
	yC.re = y;
	yC.im = 0;
	
	mul(z, x, yC);
	
}/* mul */


// z = x * y
void mul(fpComplex& z, const fpComplex& x, double y)
{
	fpComplex		yC;
	
	yC.re = y;
	yC.im = 0;
	
	mul(z, x, yC);
	
}/* mul */


// z = x * y
void mul(fpComplex& z, const fp& x, const fpComplex& y)
{
	fpComplex		xC;
	
	xC.re = x;
	xC.im = 0;
	
	mul(z, xC, y);
	
}/* mul */


// z = x * y
void mul(fpComplex& z, const mb& x, const fpComplex& y)
{
	fpComplex		xC;
	
	xC.re = x;
	xC.im = 0;
	
	mul(z, xC, y);
	
}/* mul */


// z = x * y
void mul(fpComplex& z, double x, const fpComplex& y)
{
	fpComplex		xC;
	
	xC.re = x;
	xC.im = 0;
	
	mul(z, xC, y);
	
}/* mul */


// z = x / y
void div(fpComplex& z, const fpComplex& x, const fpComplex& y)
{
	fp		ac, bd, den; // x=a+ib  y = c+id
	
	if(!y)
	{
		z.re = z.im = 0;
		return;
	}
	
	if(x==y)
	{
		z.re = 1.;
		z.im = 0;
		return;
	}
	
	if(!x.im.i.n && !y.im.i.n)
	{
		z.re = x.re/y.re;
		z.im = 0;
		return;
	}
	
	if(!x.re.i.n && !y.re.i.n)
	{
		z.re = x.im/y.im;
		z.im = 0;
		return;
	}
	
	if(!x.re.i.n && !y.im.i.n)
	{
		z.re = 0;
		z.im = x.im/y.re;
		return;
	}
	
	if(!x.im.i.n && !y.re.i.n)
	{
		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 */


// z = x / y
void div(fpComplex& z, const fpComplex& x, const fp& y)
{
	fpComplex		yC;
	
	yC.re = y;
	yC.im = 0;
	
	div(z, x, yC);
	
}/* div */


// z = x / y
void div(fpComplex& z, const fpComplex& x, const mb& y)
{
	fpComplex		yC;
	
	yC.re = y;
	yC.im = 0;
	
	div(z, x, yC);
	
}/* div */


// z = x / y
void div(fpComplex& z, const fpComplex& x, double y)
{
	fpComplex		yC;
	
	yC.re = y;
	yC.im = 0;
	
	div(z, x, yC);
	
}/* div */


// z = x / y
void div(fpComplex& z, const fp& x, const fpComplex& y)
{
	fpComplex		xC;
	
	xC.re = x;
	xC.im = 0;
	
	div(z, xC, y);
	
}/* div */


// z = x / y
void div(fpComplex& z, const mb& x, const fpComplex& y)
{
	fpComplex		xC;
	
	xC.re = x;
	xC.im = 0;
	
	div(z, xC, y);
	
}/* div */


// z = x / y
void div(fpComplex& z, double x, const fpComplex& y)
{
	fpComplex		xC;
	
	xC.re = x;
	xC.im = 0;
	
	div(z, xC, y);
	
}/* div */

