
/* complex.c */

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

Complex zadd(Complex z1, Complex z2)
{
        Complex temp;

	temp.r = z1.r + z2.r;
	temp.i = z1.i + z2.i;
	return temp;
}

Complex zsub(Complex z1, Complex z2)
{
        Complex temp;

	temp.r = z1.r - z2.r;
	temp.i = z1.i - z2.i;
	return temp;
}

Complex zconj(Complex z)
{
	z.i = -z.i;
	return z;
}

Complex zmult(Complex z1, Complex z2)
{
	Complex temp;

	temp.r = z1.r*z2.r - z1.i*z2.i;
	temp.i = z1.r*z2.i + z1.i*z2.r;
	return temp;
}

Complex zdiv(Complex z1, Complex z2)
{
	Complex temp;
	double div;

	div = z2.r*z2.r + z2.i*z2.i;
	temp.r = (z1.r*z2.r + z1.i*z2.i) / div;
	temp.i = (z1.i*z2.r - z1.r*z2.i) / div;
	return temp;
}

Complex zlog(Complex z)
{
	Complex temp;
	
	temp.r = log(sqrt(z.r*z.r + z.i*z.i));
	temp.i = atan2(z.i,z.r); /* Between -PI and PI */
	return temp;
}

Zarray4 zaddm(Zarray4 zm1, Zarray4 zm2)
{
        Complex zadd(Complex, Complex);
        Zarray4 temp;

	temp.e[0] = zadd(zm1.e[0],zm2.e[0]);
	temp.e[1] = zadd(zm1.e[1],zm2.e[1]);
	temp.e[2] = zadd(zm1.e[2],zm2.e[2]);
	temp.e[3] = zadd(zm1.e[3],zm2.e[3]);
	return temp;
}

Zarray4 zmultsm(Zarray4 zm1, Zarray4 zm2)
{
        Complex zmult(Complex, Complex);
        Zarray4 temp;

	temp.e[0] = zmult(zm1.e[0],zm2.e[0]);
	temp.e[1] = zmult(zm1.e[1],zm2.e[1]);
	temp.e[2] = zmult(zm1.e[2],zm2.e[2]);
	temp.e[3] = zmult(zm1.e[3],zm2.e[3]);
	return temp;
}

Zarray4 zmultm(Zarray4 zm1, Zarray4 zm2)
{
	Zarray4 temp;
	Complex u,v;
	Complex zadd(Complex, Complex);
	Complex zmult(Complex, Complex);

	u = zmult(zm1.e[0],zm2.e[0]); v = zmult(zm1.e[1],zm2.e[2]);
	temp.e[0] = zadd(u,v);
	u = zmult(zm1.e[0],zm2.e[1]); v = zmult(zm1.e[1],zm2.e[3]);
	temp.e[1] = zadd(u,v);
	u = zmult(zm1.e[2],zm2.e[0]); v = zmult(zm1.e[3],zm2.e[2]);
	temp.e[2] = zadd(u,v);
	u = zmult(zm1.e[2],zm2.e[1]); v = zmult(zm1.e[3],zm2.e[3]);
	temp.e[3] = zadd(u,v);
	return temp;
}

Zarray2 zmultv(Zarray4 zm, Zarray2 zv)
{
	Zarray2 temp;
	Complex u,v;
	Complex zadd(Complex, Complex);
	Complex zmult(Complex, Complex);

	u = zmult(zm.e[0],zv.e[0]); v = zmult(zm.e[1],zv.e[1]);
	temp.e[0] = zadd(u,v);
	u = zmult(zm.e[2],zv.e[0]); v = zmult(zm.e[3],zv.e[1]);
	temp.e[1] = zadd(u,v);
	return temp;
}

Complex zdetm(Zarray4 zm)
{
	Complex u,v;
	Complex zsub(Complex, Complex);
	Complex zmult(Complex, Complex);

	u = zmult(zm.e[0],zm.e[3]); v = zmult(zm.e[1],zm.e[2]);

	return zsub(u,v);
}

Zarray4 zinvm(Zarray4 zm)
{
	Complex d, u;
	Complex zdiv(Complex, Complex);
	Complex zdetm(Zarray4);

	d = zdetm(zm);

	u = zm.e[0]; 
	zm.e[0].r = zm.e[3].r; zm.e[0].i = zm.e[3].i;
	zm.e[1].r = -zm.e[1].r; zm.e[1].i = -zm.e[1].i;
	zm.e[2].r = -zm.e[2].r; zm.e[2].i = -zm.e[2].i;
	zm.e[3].r = u.r; zm.e[3].i = u.i;

	zm.e[0] = zdiv(zm.e[0],d);
	zm.e[1] = zdiv(zm.e[1],d);
	zm.e[2] = zdiv(zm.e[2],d);
	zm.e[3] = zdiv(zm.e[3],d);

	return zm;
}
