/* MatrixFunctions.cp */

#include "MatrixFunctions.h"
#include "ComplexFunctions.h"
#include "Complex.h"
#include <stdlib.h>


static inline void assignZeroM(Complex **z, int n)
{
	long	i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = 0;
			z[i][j].Im = 0;
		}

}/* assignZeroM */


static inline void assignUnitM(Complex **z, int n)
{
	long	i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(i==j)
				z[i][j].Re = 1;
			else
				z[i][j].Re = 0;
				
			z[i][j].Im = 0;
		}

}/* assignUnitM */


static inline void equateM_M(Complex **z, Complex **x, int n)
{
	long	i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = x[i][j].Re;
			z[i][j].Im = x[i][j].Im;
		}
	
}/* equateM_M */


static inline void AddM_M(Complex **z, Complex **x, Complex **y, int n)
{
	long	i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			Add(z[i][j], x[i][j], y[i][j]);
		
}/* AddM_M */


// adds a complex X unit matrix to z
static inline void AddM_C(Complex **z, Complex x, int n)
{
	long	i;
	
	for(i=0;i<n;i++)
		Add(z[i][i], z[i][i], x);
	
}/* AddM_C */


// subtracts a complex X unit matrix from z
static inline void SubM_C(Complex **z, Complex x, int n)
{
	long	i;
	
	for(i=0;i<n;i++)
		Sub(z[i][i], z[i][i], x);
	
}/* SubM_C */


// z = x*y
static inline void MulM_M(Complex **z, Complex **x, Complex **y, Complex **t, int n)
{
	long		i, j, k;
	Complex		temp;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			assignZeroC(t[i][j]);
			for(k=0;k<n;k++)
			{
				Mul(temp, x[i][k], y[k][j]);
				Add(t[i][j], t[i][j], temp);
			}
		}
		
	equateM_M(z, t, n);

}/* MulM_M */


static inline void MulM_C(Complex **z, Complex x, int n)
{
	long		i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			Mul(z[i][j], z[i][j], x);
		
}/* MulM_C */

static inline void DivM_D(Complex **z, double x, int n)
{
	long	i, j;
	double	tempD;
	
	if(x==0)
		return;
		
	tempD = 1/x;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = tempD*z[i][j].Re;
			z[i][j].Im = tempD*z[i][j].Im;
		}
	
}/* DivM_D */


void exp_tM(Complex **R, Complex ***P, Complex t, Complex **M, Complex *root, int *m, int k, int n,
			Complex **tempM, Complex **tempM1, Complex **sum, Complex **current)
{
	long		i, j;
	Complex		tempC;
	
	assignZeroM(R, n);
	
	for(i=0;i<k;i++)
	{
		assignUnitM(current, n);
		assignUnitM(sum, n);
		
		for(j=0;j<m[i]-1;j++)
		{
			equateM_M(tempM, M, n);
			SubM_C(tempM, root[i], n);
			MulM_M(current, current, tempM, tempM1, n);
			MulM_C(current, t, n);
			DivM_D(current, j+1, n);
			AddM_M(sum, sum, current, n);
		}
		
		MulM_M(tempM, sum, P[i], tempM1, n);
		Mul(tempC, root[i], t);
		myExp(tempC.Re, tempC.Im, tempC.Re, tempC.Im);
		MulM_C(tempM, tempC, n);
		AddM_M(R, R, tempM, n);
	}
	
}/* exp_tM */


void log_M(Complex **R, Complex ***P, Complex ***Omega, Complex *root, int *m, int k, int n,
			Complex **tempM, Complex **tempM1, Complex **sum, Complex **current)
{
	long		i, j;
	Complex		tempC;
		
	assignZeroM(R, n);
	
	for(i=0;i<k;i++)
	{
		assignZeroM(tempM, n);
		myLog(tempC.Re, tempC.Im, root[i].Re, root[i].Im);
		AddM_C(tempM, tempC, n);
		
		if(m[i]>1)
		{
			equateM_M(current, Omega[i], n);
			equateM_M(sum, Omega[i], n);
		}
		
		for(j=0;j<m[i]-2;j++)
		{
			MulM_M(current, current, Omega[i], tempM1, n);
			DivM_D(current, -(j+2), n);
			AddM_M(sum, sum, current, n);
		}
		
		if(m[i]>1)
			AddM_M(tempM, tempM, sum, n);
		
		MulM_M(tempM, tempM, P[i], tempM1, n);
		
		AddM_M(R, R, tempM, n);
	}

}/* log_M */
	
// this will work for a zero root also of any multiplicity
void M_t(Complex **R, Complex ***P, Complex t, Complex ***Omega, Complex *root, int *m, int k, int n,
		Complex **tempM, Complex **tempM1, Complex **sum, Complex **current)
{
	int		i, j;
	Complex		tempC;
		
	assignZeroM(R, n);
	
	for(i=0;i<k;i++)
	{
		assignZeroM(tempM, n);
		myPowCC(tempC.Re, tempC.Im, root[i].Re, root[i].Im, t.Re, t.Im);
		assignUnitM(tempM, n);
		
		if(!isEqualZero(root[i]))
		{
			if(m[i]>1)
			{
				equateM_M(current, Omega[i], n);
				MulM_C(current, t, n);
				equateM_M(sum, current, n);
			}
			
			for(j=0;j<m[i]-2;j++)
			{
				MulM_M(current, current, Omega[i], tempM1, n);
				DivM_D(current, j+2, n);
				tempC.Re = t.Re - 1 - j;
				tempC.Im = t.Im;
				MulM_C(current, tempC, n);
				AddM_M(sum, sum, current, n);
			}
			
			if(m[i]>1)
				AddM_M(tempM, tempM, sum, n);
		}
			
		MulM_C(tempM, tempC, n);
		
		MulM_M(tempM, tempM, P[i], tempM1, n);
		
		AddM_M(R, R, tempM, n);
	}
	
}/* M_t */