/* myComplexMatrixCalc.cp */

#include "myTypedef.h"
#include "myComplexMatrixCalc.h"
#include "matinv.h"
#include "cmatinv.h"
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "ComplexFunctions.h"
#include "HandleArrays.h"
#include "MatrixFunctions.h"
#include "Complex.h"
#include "fpMath.h"

using namespace std;

double  zeroDouble=1.e-10;

static long min(long x, long y)
{
    if(x<=y)
        return x;
    
    return y;
    
}/* min */

static double dabs(double x)
{
	if(x<0)
		x = -x;
	
	return x;
	
}/* dabs */

static double abs2(const Complex x)
{
    return x.Re*x.Re+x.Im*x.Im;
    
}/* abs2 */

static double abs2D(const double x)
{
    return x*x;
    
}/* abs2D */

// swap row i with row j
static void swapRows(Complex **z, long nr, long nc, long i, long j)
{
    Complex     temp;
    long        k;
    
    for(k=0;k<nc;++k)
        {
            temp.Re = z[i][k].Re;
            temp.Im = z[i][k].Im;
            z[i][k].Re = z[j][k].Re;
            z[i][k].Im = z[j][k].Im;
            z[j][k].Re = temp.Re;
            z[j][k].Im = temp.Im;
        }
    
}/* swapRows */

// swap row i with row j
static void swapRowsD(double **z, long nr, long nc, long i, long j)
{
    double     temp;
    long        k;
    
    for(k=0;k<nc;++k)
    {
        temp = z[i][k];
        z[i][k] = z[j][k];
        z[j][k] = temp;
    }
    
}/* swapRowsD */


static bool doGaussD(double **x, long nr, long nc, long ii, long jj)
{
    double          temp;
    long            i, j, k;
    bool            isZero;
    
    if(ii<0 || jj<0 || ii>=nr || jj>=nc)
        return false;
    
    k = 0; // to avoid incorrect warning
    // is pivot zero?
    if(!x[ii][jj])
    {
        // must swap rows
        // find non-zero element below pivot in pivot column
        isZero = true;
        for(i=ii+1;i<nr;++i)
        {
            if(x[i][jj])
            {
                isZero = false;
                k = i;
                break;
            }
        }
        
        if(isZero)
            return false;
        
        // have non-zero element in row k, so swap rows ii with k
        swapRowsD(x, nr, nc, ii, k);
    }
    // now we do the Gauss Elim. operation
    // pivot element xt[ii][jj] is not zero
   
    
    // now for every k not ii subtract row ii multiplied by xt[k][jj] from row k provided x[k][jj] != 0
    for(k=0;k<nr;++k)
    {
        if(k!=ii && x[k][jj])
        {
            temp = x[k][jj] / x[ii][jj];
            for(j=0;j<nc;++j)
            {
                x[k][j] = x[k][j] - (x[ii][j] * temp);
                if(dabs(x[k][j])<zeroDouble)
                {
                    x[k][j] = 0.;
                }
            }
        }
    }
    return true;
    
}/* doGaussD */

static bool doReduceD(double **x, long nr, long nc, long ii, long jj)
{
    double          temp;
    long            i, j, k;
    bool            isZero;
    
    if(ii<0 || jj<0 || ii>=nr || jj>=nc)
        return false;
    
    k = 0; // to avoid incorrect warning
    // is pivot zero?
    if(!x[ii][jj])
    {
        // must swap rows
        // find non-zero element below pivot in pivot column
        isZero = true;
        for(i=ii+1;i<nr;++i)
        {
            if(x[i][jj])
            {
                isZero = false;
                k = i;
                break;
            }
        }
        
        if(isZero)
            return false;
        
        // have non-zero element in row k, so swap rows ii with k
        swapRowsD(x, nr, nc, ii, k);
    }
    // now we do the reduce operation
    // pivot element xt[ii][jj] is not zero
    // divide each element of row ii by xt[ii][jj] making xt[ii][jj] = 1
    temp = x[ii][jj];
    for(j=0;j<nc;++j)
    {
        x[ii][j] = x[ii][j] / temp;
    }
    
    // now for every k not ii subtract row ii multiplied by xt[k][jj] from row k provided x[k][jj] != 0
    for(k=0;k<nr;++k)
    {
        if(k!=ii && x[k][jj])
        {
            temp = x[k][jj];
            for(j=0;j<nc;++j)
            {
                x[k][j] = x[k][j] - (x[ii][j] * temp);
                if(dabs(x[k][j])<zeroDouble)
                {
                    x[k][j] = 0.;
                }
            }
        }
    }
    return true;
    
}/* doReduceD */

// z must be inited by caller
static void EquateArrays(Complex **z, long nr, long nc, Complex **x)
{
    long    i, j;
    
    for(i=0;i<nr;++i)
        for(j=0;j<nc;++j)
        {
            z[i][j].Re = x[i][j].Re;
            z[i][j].Im = x[i][j].Im;
        }
    
}/* EquateArrays */

// z must be inited by caller
static void EquateArraysD(double **z, long nr, long nc, double **x)
{
    long    i, j;
    
    for(i=0;i<nr;++i)
        for(j=0;j<nc;++j)
            z[i][j] = x[i][j];
    
}/* EquateArraysD */


static void PrettyFormat(Complex **z, int n)
{
	double		averageSize;
	long		i, j;
	
	if(n<=0)
		return;
		
	averageSize = 0;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			averageSize = averageSize + dabs(z[i][j].Re) + dabs(z[i][j].Im);
			
	averageSize = 1.e-16*averageSize/n;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(dabs(z[i][j].Re)<averageSize)
				z[i][j].Re = 0;
						
			if(dabs(z[i][j].Im)<averageSize)
				z[i][j].Im = 0;
		}

}/* PrettyFormat */


static inline bool isEqualCC(Complex& z, const Complex& x)
{
	if(z.Re==x.Re && z.Im==x.Im)
		return true;
	else
		return false;
	
}/* isEqualCC */

/*
static inline void EquateCC(Complex& z, const Complex& x)
{
	z.Re = x.Re;
	z.Im = x.Im;

}*//* EquateCC */


static inline void NegateCC(Complex& z)
{
	z.Re = -z.Re;
	z.Im = -z.Im;

}/* NegateCC */


static inline void AddCC(Complex& z, const Complex& x, const Complex& y)
{
	z.Re = x.Re + y.Re;
	z.Im = x.Im + y.Im;

}/* AddCC */


static inline void SubCC(Complex& z, const Complex& x, const Complex& y)
{
	z.Re = x.Re - y.Re;
	z.Im = x.Im - y.Im;

}/* SubCC */


static inline void MulCC(Complex& z, const Complex& x, const Complex& y)
{
	double		ac, bd;
	
	if(x.Re!=0 && x.Im!=0 && y.Re!=0 && y.Im!=0)
	{
		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==0)
	{
		ac = -x.Im*y.Im;
		z.Im = x.Im*y.Re;
		z.Re = ac;
		return;
	}
	
	if(x.Im==0)
	{
		ac = x.Re*y.Re;
		z.Im = x.Re*y.Im;
		z.Re = ac;
		return;
	}
	
	if(y.Re==0)
	{
		ac = -x.Im*y.Im;
		z.Im = x.Re*y.Im;
		z.Re = ac;
		return;
	}
	
	if(y.Im==0)
	{
		ac = x.Re*y.Re;
		z.Im = x.Im*y.Re;
		z.Re = ac;
		return;
	}
	
}/* MulCC */


static inline void DivCC(Complex& z, const Complex& x, const Complex& y)
{
    double		ac, bd, den;
    
    if(y.Re==0 && y.Im==0)
    {
        z.Re = NAN;
        z.Im = NAN;
        return;
    }
    
    if(x.Im==0 && y.Im==0)
    {
        z.Re = x.Re/y.Re;
        z.Im = 0;
        return;
    }
    
    if(x.Re==0 && y.Re==0)
    {
        z.Re = x.Im/y.Im;
        z.Im = 0;
        return;
    }
    
    if(x.Re==0 && y.Im==0)
    {
        z.Im = x.Im/y.Re;
        z.Re = 0;
        return;
    }
    
    if(x.Im==0 && y.Re==0)
    {
        z.Im = -x.Re/y.Im;
        z.Re = 0;
        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;
    
}/* DivCC */


static inline void DivCC_D(Complex& z, double x)
{
	z.Re = z.Re/x;
	z.Im = z.Im/x;

}/* DivCC_D */


// z = x*y; z<>x, z<>y
static inline void MulCM_CM(Complex **z, Complex **x, Complex **y, int n)
{
	long		i, j, k;
	Complex		temp;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = 0;
			z[i][j].Im = 0;
			for(k=0;k<n;k++)
			{
				MulCC(temp, x[i][k], y[k][j]);
				AddCC(z[i][j], z[i][j], temp);
			}
		}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
			
}/* MulCM_CM */


// z = x*y
static inline void MulCM_M(Complex **z, Complex **x, double **y, int n)
{
	long			i, j, k;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = 0;
			z[i][j].Im = 0;
			for(k=0;k<n;k++)
			{
				z[i][j].Re = z[i][j].Re + x[i][k].Re*y[k][j];
				z[i][j].Im = z[i][j].Im + x[i][k].Im*y[k][j];
			}
		}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
			
}/* MulCM_M */


// z = x*y
static inline void MulM_CM(Complex **z, double **x, Complex **y, int n)
{
	long			i, j, k;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = 0;
			z[i][j].Im = 0;
			for(k=0;k<n;k++)
			{
				z[i][j].Re = z[i][j].Re + x[i][k]*y[k][j].Re;
				z[i][j].Im = z[i][j].Im + x[i][k]*y[k][j].Im;
			}
		}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
			
}/* MulM_CM */


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

}/* EquateCM_CM */

// z = z/k
static inline void DivCM_D(Complex **z, int k, int n)
{
	long	i, j;
	
	if(k==0)
		return;
		
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			DivCC_D(z[i][j], (double)k);

}/* DivCM_D */

// z = -z
static inline void NegateCM(Complex **z, int n)
{
	long	i, j;
		
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			NegateCC(z[i][j]);

}/* NegateCM */

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


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


// does z==x?
static inline bool isEqualCM_CM(Complex **z, Complex **x, int n)
{
	long		i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(!isEqualCC(z[i][j], x[i][j]))
				return false;
		}
		
	return true;
	
}/* isEqualCM_CM */

// z = x*x
void mySqrCM(Complex **z, Complex **x, int n)
{
	if(n<=0)
		return;
		
	if(n==1)
	{
		MulCC(z[0][0], x[0][0], x[0][0]);
		return;
	}
	
	MulCM_CM(z, x, x, n);
		
}/* mySqrCM */


// z = sqrt(x)
void mySqrtCM(Complex **z, Complex **x, int n)
{
	Complex			*root;
	int			*multiplicity;
	int			k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	Complex			t;
	bool			isSingular;
	bool			isFillArraysGood;

	if(n<=0)
		return;
		
	if(n==1)
	{
		mySqrt(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	if(!isFillArraysGood)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		EquateCM_CM(z, x, n);
		return;
	}
	
	t.Re = .5;
	t.Im = 0;
	M_t(R, P, t, Omega, root, multiplicity, k, n, tempM, tempM1, sum, current);
	
	EquateCM_CM(z, R, n);
	
	PrettyFormat(z, n);
	
	killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		
}/* mySqrtCM */


// z = sqrt(x) for a real matrix x and complex matrix z
void mySqrtM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	mySqrtCM(z, z, n);
		
}/* mySqrtM */


// z = exp(x) x and z must be square of order n
void myExpCM(Complex **z, Complex **x, int n)
{
	long			i, j, k;
	Complex			**temp, **currentTerm, **zCurrent;
	double			temp1;
	double			maxEl;
	int             m;
	Complex			**xt;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		temp1 = exp(x[0][0].Re);
		z[0][0].Re = temp1*cos(x[0][0].Im);
		z[0][0].Im = temp1*sin(x[0][0].Im);
		return;
	}
	
	// find maximum element of x in absolute value
	maxEl = 0;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if(AbsC(x[i][j])>maxEl)
				maxEl = AbsC(x[i][j]);
	
	if(maxEl==0)
	{
		// z = unit matrix
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
				if(i==j)
				{
					z[i][j].Re = 1.;
					z[i][j].Im = 0;
				}
				else
				{
					z[i][j].Re = 0;
					z[i][j].Im = 0;
				}
		return;
	}
	
	xt = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		xt[i] = (Complex *)malloc(n*sizeof(Complex));
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			xt[i][j].Re = x[i][j].Re;
			xt[i][j].Im = x[i][j].Im;
		}
	
	m = 1;
	if(maxEl<(double)0x7FFFFFFF)
	{
		m = maxEl;
		m++;
		if(m>1)
		{
			// divide matrix xt by m
			for(i=0;i<n;i++)
				for(j=0;j<n;j++)
				{
					xt[i][j].Re = xt[i][j].Re/m;
					xt[i][j].Im = xt[i][j].Im/m;
				}
		}
	}

	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	currentTerm = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		currentTerm[i] = (Complex *)malloc(n*sizeof(Complex));
		
	zCurrent = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		zCurrent[i] = (Complex *)malloc(n*sizeof(Complex));
	
	
	
	// make z, currentTerm, and zCurrent Unit matrices
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if(i==j)
			{
				z[i][j].Re = currentTerm[i][j].Re = zCurrent[i][j].Re = 1;
				z[i][j].Im = currentTerm[i][j].Im = zCurrent[i][j].Im = 0;
			}
			else
			{
				z[i][j].Re = currentTerm[i][j].Re = zCurrent[i][j].Re = 0;
				z[i][j].Im = currentTerm[i][j].Im = zCurrent[i][j].Im = 0;
			}
	
	
	for(k=1; ;k++)
	{
		MulCM_CM(temp, xt, currentTerm, n);
		EquateCM_CM(currentTerm, temp, n);
		DivCM_D(currentTerm, k, n);
		
		AddCM_CM(z, z, currentTerm, n);
		
		if(isEqualCM_CM(z, zCurrent, n))
			break;
		
		EquateCM_CM(zCurrent, z, n);
	}
	
	if(m>1)
	{
		myPowerCM(z, zCurrent, m, n);
	}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	for(i=0;i<n;i++)
	{
		free(xt[i]);
		free(temp[i]);
		free(currentTerm[i]);
		free(zCurrent[i]);
	}
	
	free(xt);
	free(temp);
	free(currentTerm);
	free(zCurrent);
	
}/* myExpCM */

// z=log(x)
void myLogCM(Complex **z, Complex **x, int n)
{
	Complex			*root;
	int             *multiplicity;
	int             i, j, k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	bool			isSingular;
	bool			isFillArraysGood;

	if(n<=0)
		return;
		
	if(n==1)
	{
		myLog(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	if(!isFillArraysGood || isSingular)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		EquateCM_CM(z, x, n);
		return;
	}
	
	
	log_M(R, P, Omega, root, multiplicity, k, n, tempM, tempM1, sum, current);
	
	EquateCM_CM(z, R, n);
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	PrettyFormat(z, n);
	
	killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
    
}/* myLogCM */


// z=log(x)
void myLogM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myLogCM(z, z, n);

}/* myLogM */


// z = sin(x) x and z must be square of order n
void mySinCM(Complex **z, Complex **x, int n)
{
	long			i, j, k;
	Complex		**temp, **currentTerm, **zCurrent, **xSqr;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		z[0][0].Re = sin(x[0][0].Re)*cosh(x[0][0].Im);
		z[0][0].Im = cos(x[0][0].Re)*sinh(x[0][0].Im);
		return;
	}
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	currentTerm = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		currentTerm[i] = (Complex *)malloc(n*sizeof(Complex));
		
	zCurrent = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		zCurrent[i] = (Complex *)malloc(n*sizeof(Complex));
		
	xSqr = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		xSqr[i] = (Complex *)malloc(n*sizeof(Complex));
	
	// Equate z, currentTerm, zCurrent, and temp to x
	EquateCM_CM(z, x, n);
	EquateCM_CM(currentTerm, x, n);
	EquateCM_CM(zCurrent, x, n);
	
	MulCM_CM(xSqr, x, x, n);
	
	for(k=2; ;k+=2)
	{
		MulCM_CM(temp, xSqr, currentTerm, n);
		EquateCM_CM(currentTerm, temp, n);
		DivCM_D(currentTerm, -k*(k+1), n);
		
		AddCM_CM(z, z, currentTerm, n);
		
		if(isEqualCM_CM(z, zCurrent, n))
			break;
		
		EquateCM_CM(zCurrent, z, n);
	}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	for(i=0;i<n;i++)
	{
		free(temp[i]);
		free(currentTerm[i]);
		free(zCurrent[i]);
		free(xSqr[i]);
	}
	
	free(temp);
	free(currentTerm);
	free(zCurrent);
	free(xSqr);
	
}/* mySinCM */


// z = cos(x) x and z must be square of order n
void myCosCM(Complex **z, Complex **x, int n)
{
	long			i, j, k;
	Complex		**temp, **currentTerm, **zCurrent, **xSqr;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		z[0][0].Re = cos(x[0][0].Re)*cosh(x[0][0].Im);
		z[0][0].Im = -sin(x[0][0].Re)*sinh(x[0][0].Im);
		return;
	}
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	currentTerm = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		currentTerm[i] = (Complex *)malloc(n*sizeof(Complex));
		
	zCurrent = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		zCurrent[i] = (Complex *)malloc(n*sizeof(Complex));
		
	xSqr = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		xSqr[i] = (Complex *)malloc(n*sizeof(Complex));
		
	// make z, currentTerm, and zCurrent Unit matrices
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if(i==j)
			{
				z[i][j].Re = currentTerm[i][j].Re = zCurrent[i][j].Re = 1;
				z[i][j].Im = currentTerm[i][j].Im = zCurrent[i][j].Im = 0;
			}
			else
			{
				z[i][j].Re = currentTerm[i][j].Re = zCurrent[i][j].Re = 0;
				z[i][j].Im = currentTerm[i][j].Im = zCurrent[i][j].Im = 0;
			}
	
	MulCM_CM(xSqr, x, x, n);
	
	for(k=1; ;k+=2)
	{
		MulCM_CM(temp, xSqr, currentTerm, n);
		EquateCM_CM(currentTerm, temp, n);
		DivCM_D(currentTerm, -k*(k+1), n);
		
		AddCM_CM(z, z, currentTerm, n);
		
		if(isEqualCM_CM(z, zCurrent, n))
			break;
		
		EquateCM_CM(zCurrent, z, n);
	}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	for(i=0;i<n;i++)
	{
		free(temp[i]);
		free(currentTerm[i]);
		free(zCurrent[i]);
		free(xSqr[i]);
	}
	
	free(temp);
	free(currentTerm);
	free(zCurrent);
	free(xSqr);
	
}/* myCosCM */


// z = tan(x) x and z must be square of order n
// note from power series that [sin(x), cos(x)] = 0
// so tan(x) = Inverse(cos(x))*sin(x) = sin(x)*Inverse(cos(x))
void myTanCM(Complex **z, Complex **x, int n)
{
	Complex		**theSin, **theCos;
	Complex		determinant, temp1, temp2;
	long			i, j;
	double			sinx, cosx, sinhy, coshy;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		sinx = sin(x[0][0].Re);
		cosx = cos(x[0][0].Re);
		sinhy = sinh(x[0][0].Im);
		coshy = cosh(x[0][0].Im);
		temp1.Re = sinx*coshy;
		temp1.Im = cosx*sinhy;
		temp2.Re = cosx*coshy;
		temp2.Im = -sinx*sinhy;
		DivCC(z[0][0], temp1, temp2);
		return;
	}
	
	theSin = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		theSin[i] = (Complex *)malloc(n*sizeof(Complex));
		
	theCos = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		theCos[i] = (Complex *)malloc(n*sizeof(Complex));
		
	mySinCM(theSin, x, n);
	myCosCM(theCos, x, n);
	
	determinant = cmatinv(theCos, n); //theCos is now Inverted
	
	if(!(determinant.Re==0 && determinant.Im==0))
		MulCM_CM(z, theCos, theSin, n);
	
	for(i=0;i<n;i++)
	{
		free(theSin[i]);
		free(theCos[i]);
	}
	
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
    
	free(theSin);
	free(theCos);
	
}/* myTanCM */


// z = aSin(x) = -i*log(sqrt(1-x^2) + i*x)
void myASinCM(Complex **z, Complex **x, int n)
{
	Complex			*root;
	int			*multiplicity;
	int			i, j, k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	bool			isSingular;
	bool			isFillArraysGood;

	if(n<=0)
		return;
		
	if(n==1)
	{
		myASin(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	if(!isFillArraysGood)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		EquateCM_CM(z, x, n);
		return;
	}
	
	MulCM_CM(tempM, x, x, n);
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(i==j)
			{
				tempM[i][j].Re = 1 - tempM[i][j].Re;
				tempM[i][j].Im = -tempM[i][j].Im;
			}
			else
			{
				tempM[i][j].Re = -tempM[i][j].Re;
				tempM[i][j].Im = -tempM[i][j].Im;
			}
		}
		
	mySqrtCM(R, tempM, n);
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			tempM[i][j].Re = R[i][j].Re - x[i][j].Im;
			tempM[i][j].Im = R[i][j].Im + x[i][j].Re;
		}
		
	myLogCM(R, tempM, n);
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = R[i][j].Im;
			z[i][j].Im = -R[i][j].Re;
		}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	PrettyFormat(z, n);
	
	killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);

}/* myASinCM */


// z = aCos(x) = -i*log(sqrt(x^2-1) + x)
void myACosCM(Complex **z, Complex **x, int n)
{
	Complex			*root;
	int			*multiplicity;
	int			i, j, k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	bool			isSingular;
	bool			isFillArraysGood;

	if(n<=0)
		return;
		
	if(n==1)
	{
		myACos(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	if(!isFillArraysGood)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		EquateCM_CM(z, x, n);
		return;
	}
	
	MulCM_CM(tempM, x, x, n);
	for(i=0;i<n;i++)
		tempM[i][i].Re = tempM[i][i].Re - 1;
	
		
	mySqrtCM(R, tempM, n);
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			tempM[i][j].Re = R[i][j].Re + x[i][j].Re;
			tempM[i][j].Im = R[i][j].Im + x[i][j].Im;
		}
		
	myLogCM(R, tempM, n);
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = R[i][j].Im;
			z[i][j].Im = -R[i][j].Re;
		}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	PrettyFormat(z, n);
	
	killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);

}/* myACosCM */


// z = aTan(x) = -(i/2)*log((1+ix)/(1-ix))
void myATanCM(Complex **z, Complex **x, int n)
{
	Complex			*root;
	int			*multiplicity;
	int			i, j, k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	bool			isSingular;
	bool			isFillArraysGood;
	Complex			determinant;

	if(n<=0)
		return;
		
	if(n==1)
	{
		myATan(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	if(!isFillArraysGood)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		EquateCM_CM(z, x, n);
		return;
	}
	
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			if(i==j)
			{
				tempM[i][j].Re = 1 - x[i][j].Im;
				tempM[i][j].Im = x[i][j].Re;
				
				R[i][j].Re = 1 + x[i][j].Im;
				R[i][j].Im = -x[i][j].Re;
				
			}
			else
			{
				tempM[i][j].Re = -x[i][j].Im;
				tempM[i][j].Im = x[i][j].Re;
				
				R[i][j].Re = x[i][j].Im;
				R[i][j].Im = -x[i][j].Re;
			}
		}
		
	determinant = cmatinv(R, n);
	
	if(determinant.Re==0 && determinant.Im==0)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		EquateCM_CM(z, x, n);
		return;
	}
	
	MulCM_CM(tempM1, tempM, R, n);
		
	myLogCM(R, tempM1, n);
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = R[i][j].Im/2;
			z[i][j].Im = -R[i][j].Re/2;
		}
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	PrettyFormat(z, n);
	
	killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);

}/* myATanCM */


void myASinM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myASinCM(z, z, n);

}/* myASinM */


void myACosM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myACosCM(z, z, n);

}/* myACosM */


void myATanM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myATanCM(z, z, n);

}/* myATanM */


void mySinhCM(Complex **z, Complex **x, int n)
{
	Complex		**temp;
	long			i, j;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		z[0][0].Re = sinh(x[0][0].Re)*cos(x[0][0].Im);
		z[0][0].Im = cosh(x[0][0].Re)*sin(x[0][0].Im);
		return;
	}
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	EquateCM_CM(temp, x, n); // temp = x
	NegateCM(temp, n); // temp = -x
	myExpCM(z, temp, n); // z = exp(-x)
	EquateCM_CM(temp, z, n); // temp = exp(-x)
	
	myExpCM(z, x, n); // z = exp(x)
	SubCM_CM(z, z, temp, n); // z = exp(x) - exp(-x)
	DivCM_D(z, 2, n); // z = sinh(x)
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
		
	for(i=0;i<n;i++)
		free(temp[i]);
	
	free(temp);
	
}/* mySinhCM */


void myCoshCM(Complex **z, Complex **x, int n)
{
	Complex		**temp;
	long			i, j;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		z[0][0].Re = cosh(x[0][0].Re)*cos(x[0][0].Im);
		z[0][0].Im = sinh(x[0][0].Re)*sin(x[0][0].Im);
		return;
	}
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	EquateCM_CM(temp, x, n); // temp = x
	NegateCM(temp, n); // temp = -x
	myExpCM(z, temp, n); // z = exp(-x)
	EquateCM_CM(temp, z, n); // temp = exp(-x)
	
	myExpCM(z, x, n); // z = exp(x)
	AddCM_CM(z, z, temp, n); // z = exp(x) + exp(-x)
	DivCM_D(z, 2, n); // z = sinh(x)
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
		
	for(i=0;i<n;i++)
		free(temp[i]);
	
	free(temp);
	
}/* myCoshCM */


// z = tanh(x) x and z must be square of order n
// note from power series that [sinh(x), cosh(x)] = 0
// so tanh(x) = Inverse(cosh(x))*sinh(x) = sinh(x)*Inverse(cosh(x))
void myTanhCM(Complex **z, Complex **x, int n)
{
	Complex		**theSinh, **theCosh;
	Complex		determinant, temp1, temp2;
	long			i, j;
	double			sinhx, coshx, siny, cosy;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		sinhx = sinh(x[0][0].Re);
		coshx = cosh(x[0][0].Re);
		siny = sin(x[0][0].Im);
		cosy = cos(x[0][0].Im);
		temp1.Re = sinhx*cosy;
		temp1.Im = coshx*siny;
		temp2.Re = coshx*cosy;
		temp2.Im = sinhx*siny;
		DivCC(z[0][0], temp1, temp2);
		return;
	}
	
	theSinh = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		theSinh[i] = (Complex *)malloc(n*sizeof(Complex));
		
	theCosh = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		theCosh[i] = (Complex *)malloc(n*sizeof(Complex));
		
	mySinhCM(theSinh, x, n);
	myCoshCM(theCosh, x, n);
	
	determinant = cmatinv(theCosh, n); //theCosh is now Inverted
	
	if(!(determinant.Re==0 && determinant.Im==0))
		MulCM_CM(z, theCosh, theSinh, n);
	
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
    
	for(i=0;i<n;i++)
	{
		free(theSinh[i]);
		free(theCosh[i]);
	}
	
	free(theSinh);
	free(theCosh);
	
}/* myTanhCM */


// z = aSinh(x) = -i*aSin(i*x)
void myASinhCM(Complex **z, Complex **x, int n)
{
	long		i, j;
	double		temp;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		myASinhC(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			temp = x[i][j].Re;  // in case z=x
			z[i][j].Re = -x[i][j].Im;
			z[i][j].Im = temp;
		}
		
	myASinCM(z, z, n);
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			temp = z[i][j].Re;
			z[i][j].Re = z[i][j].Im;
			z[i][j].Im = -temp;
		}
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
		
}/* myASinhCM */


// z = aCosh(x) = -i*aCos(x)
void myACoshCM(Complex **z, Complex **x, int n)
{
	long		i, j;
	double		temp;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		myACoshC(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	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;
		}
		
	myACosCM(z, z, n);
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			temp = z[i][j].Re;
			z[i][j].Re = z[i][j].Im;
			z[i][j].Im = -temp;
		}
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
		
}/* myACoshCM */


// z = aTanh(x) = -i*aTan(i*x)
void myATanhCM(Complex **z, Complex **x, int n)
{
	long		i, j;
	double		temp;
	
	if(n<=0)
		return;
		
	if(n==1)
	{
		myATanhC(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im);
		return;
	}
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			temp = x[i][j].Re;  // in case z=x
			z[i][j].Re = -x[i][j].Im;
			z[i][j].Im = temp;
		}
		
	myATanCM(z, z, n);
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			temp = z[i][j].Re;
			z[i][j].Re = z[i][j].Im;
			z[i][j].Im = -temp;
		}
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
		
}/* myATanhCM */


void myASinhM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myASinhCM(z, z, n);

}/* myASinhM */


void myACoshM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myACoshCM(z, z, n);

}/* myACoshM */


void myATanhM(Complex **z, double **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];
			z[i][j].Im = 0;
		}
		
	myATanhCM(z, z, n);

}/* myATanhM */


// z=Pow(x, t)
void myPowCM_C(Complex **z, Complex **x, Complex t, int n)
{
	Complex			*root;
	int             *multiplicity;
	int             i, j, k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	bool			isSingular;
	bool			isFillArraysGood;

	if(n<=0)
		return;
		
	if(n==1)
	{
		myPowCC(z[0][0].Re, z[0][0].Im, x[0][0].Re, x[0][0].Im, t.Re, t.Im);
		return;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	if(!isFillArraysGood)
	{
		EquateCM_CM(z, x, n);
		return;
	}
	
	M_t(R, P, t, Omega, root, multiplicity, k, n, tempM, tempM1, sum, current);
	
	EquateCM_CM(z, R, n);
    
    // we cannot have any extremely small but non-zero elements
    for(i=0;i<n;++i)
    {
        for(j=0;j<n;++j)
        {
            if(dabs(z[i][j].Re)<zeroDouble)
                z[i][j].Re = 0;
            
            if(dabs(z[i][j].Im)<zeroDouble)
                z[i][j].Im = 0;
            
        }
    }
	
	PrettyFormat(z, n);
	
	killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
	
}/* myPowCM_C */


// z=Pow(x, t)
void myPowCM_D(Complex **z, Complex **x, double t, int n)
{
	Complex		y;
	
	y.Re = t;
	y.Im = 0;
	
	myPowCM_C(z, x, y, n);

}/* myPowCM_D */


// z=Pow(x, t)
void myPowM_C(Complex **z, double **x, Complex t, int n)
{
	long		i, j;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = x[i][j];
			z[i][j].Im = 0;
		}
	
	myPowCM_C(z, z, t, n);

}/* myPowM_C */


// z=Pow(x, t)
void myPowM_D(Complex **z, double **x, double t, int n)
{
	long		i, j;
	Complex		y;
	
	y.Re = t;
	y.Im = 0;
	
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			z[i][j].Re = x[i][j];
			z[i][j].Im = 0;
		}
	
	myPowCM_C(z, z, y, n);

}/* myPowM_D */


// z = x*y - y*x
void myCommutatorCM_CM(Complex **z, Complex **x, Complex **y, int n)
{
	Complex		**temp;
	long			i;
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	MulCM_CM(temp, y, x, n);  // temp = y*x
	MulCM_CM(z, x, y, n);  // z = x*y
	SubCM_CM(z, z, temp, n);  // z = x*y - y*x
		
	for(i=0;i<n;i++)
		free(temp[i]);
		
	free(temp);
	
}/* myCommutatorCM_CM */


// z = x*y - y*x
void myCommutatorCM_M(Complex **z, Complex **x, double **y, int n)
{
	Complex		**temp;
	long			i;
	
	temp = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		temp[i] = (Complex *)malloc(n*sizeof(Complex));
		
	MulM_CM(temp, y, x, n);  // temp = y*x
	MulCM_M(z, x, y, n);  // z = x*y
	SubCM_CM(z, z, temp, n);  // z = x*y - y*x
		
	for(i=0;i<n;i++)
		free(temp[i]);
		
	free(temp);
	
}/* myCommutatorCM_M */


// z = x*y - y*x
void myCommutatorM_CM(Complex **z, double **x, Complex **y, int n)
{
	myCommutatorCM_M(z, y, x, n);
	NegateCM(z, n);
	
}/* myCommutatorM_CM */


bool myEigenvalues(Complex **x, Complex *eigenvalue, int n)
{
	Complex			*root;
	int			*multiplicity;
	int			k;
	Complex			***P; // projection matrices
	Complex			***Omega; // need for many matrix functions
	Complex			**R, **tempM, **tempM1, **sum, **current;
	bool			isSingular;
	bool			isFillArraysGood;
	int             i, j, m;
	
	if(n<=0)
		return false;
		
	if(n==1)
	{
		EquateC_C(eigenvalue[0], x[0][0]);
		return true;
	}
	
	isFillArraysGood = FillArrays(n, x, isSingular, root, multiplicity, k, P, Omega,
						R, tempM, tempM1, sum, current);
						
	// k is the number of distinct roots
						
	if(!isFillArraysGood)
	{
		killArrays(n, root, multiplicity, P, Omega, R, tempM, tempM1, sum, current);
		return false;
	}
	/*
	if(k==n)
	{
		for(i=0;i<n;i++)
			EquateC_C(eigenvalue[i], root[i]);
	}
	*/
	
	
	m = 0;
	for(i=0;i<k;i++)
	{
		for(j=0;j<multiplicity[i];j++)
		{
			EquateC_C(eigenvalue[m], root[i]);
			m++;
		}
	}
    
    for(i=0;i<n;++i)
    {
        if(dabs(eigenvalue[i].Re)<zeroDouble)
            eigenvalue[i].Re = 0;
        
        if(dabs(eigenvalue[i].Im)<zeroDouble)
            eigenvalue[i].Im = 0;
    }
    
    // bubble sort
    bool    isDirty;
    double  temp;
    
    do
    {
        isDirty = false;
        for(i=0;i<n-1;++i)
        {
            if(eigenvalue[i].Re>eigenvalue[i+1].Re)
            {
                temp = eigenvalue[i+1].Re;
                eigenvalue[i+1].Re = eigenvalue[i].Re;
                eigenvalue[i].Re = temp;
                
                temp = eigenvalue[i+1].Im;
                eigenvalue[i+1].Im = eigenvalue[i].Im;
                eigenvalue[i].Im = temp;
                isDirty = true;
            }
        }
    }
    while (isDirty);
    
	return true;
	
}/* myEigenvalues */

// z=x^m, m an integer. x must be square
// z<>x, plugin code makes z<>x
Complex myPowerCM(Complex **z, Complex **x, int m, int n)
{
	int	i, j, ii, jj, nb;
	static unsigned int	mask; // 0x80000000
	bool	isNegative;
	Complex	det;
	Complex	**zt;
	static bool	initGood=false;
	
	if(!initGood)
	{
		mask = 1;
		mask = (mask<<31); // 0x80000000
		initGood = true;
	}
	
	det.Re = 1.;
	det.Im = 0;
	
	if (m==0)  // z = unit matrix
	{
		for(i=0;i<n;i++)
			for(j=0;j<n;j++)
				if(i==j)
				{
					z[i][j].Re = 1.;
					z[i][j].Im = 0;
				}
				else
				{
					z[i][j].Re = 0;
					z[i][j].Im = 0;
				}
		
		return det;
	}
	
	if (m==1)  // z = x
	{
		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;
			}
		
		return det;
	}
	
	if (m==-1)  // z = 1/x
	{
		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;
			}
		
		det = cmatinv(z, n);
		
		return det;
	}
	
	if (m<0)
	{
		isNegative = true;
		m = -m;
	}
	else
	{
		isNegative = false;
	}
	
	// initiate zt
	zt = (Complex **)malloc(n*sizeof(Complex *));
	for (i=0;i<n;i++)
		zt[i] = (Complex *)malloc(n*sizeof(Complex));
	
	// now m>0, make z a unit matrix
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			if(i==j)
			{
				z[i][j].Re = 1.;
				z[i][j].Im = 0;
			}
			else
			{
				z[i][j].Re = 0;
				z[i][j].Im = 0;
			}
	
	nb = NumBits(m);
	
	m = (m<<(32-nb)); // go to leading 1 bit
	
	for(i=0;i<nb;i++)
	{
		MulCM_CM(zt, z, z, n);
		for(ii=0;ii<n;ii++)
			for(jj=0;jj<n;jj++)
			{
				z[ii][jj].Re = zt[ii][jj].Re;
				z[ii][jj].Im = zt[ii][jj].Im;
			}
				
		if(m & mask)
		{
			MulCM_CM(zt, z, x, n);
			for(ii=0;ii<n;ii++)
				for(jj=0;jj<n;jj++)
				{
					z[ii][jj].Re = zt[ii][jj].Re;
					z[ii][jj].Im = zt[ii][jj].Im;
				}
		}
		
		m = (m<<1);
		
	}
	
	for(i=0;i<n;i++)
		free(zt[i]);
	
	free(zt);
	
	if(isNegative)
		det = cmatinv(z, n);
	
	return det;
	
}/* myPowerCM */


void GaussElimD(double **z, long nr, long nc, double **x)
{
    double          **xt;
    long            i, ii, jj;
    
    xt = (double**)malloc(nr*sizeof(double*));
    for(i=0;i<nr;++i)
        xt[i] = (double*)malloc(nc*sizeof(double));
    
    EquateArraysD(xt, nr, nc, x);
    
    ii = jj = 0;
    while(ii<nr && jj<nc)
    {
        if(doGaussD(xt, nr, nc, ii, jj))
        {
            ii++;
            jj++;
        }
        else
            jj++;
    }
    
    EquateArraysD(z, nr, nc, xt);
    for(i=0;i<nr;++i)
        free(xt[i]);
    free(xt);
    
}/* GaussElimD */

static bool doGauss(Complex **x, long nr, long nc, long ii, long jj)
{
    Complex         temp, temp1;
    long            i, j, k;
    bool            isZero;
    
    if(ii<0 || jj<0 || ii>=nr || jj>=nc)
        return false;
    
    k = 0; // to avoid incorrect warning
    // is pivot zero?
    if(!x[ii][jj].Re && !x[ii][jj].Im)
    {
        // must swap rows
        // find non-zero element below pivot in pivot column
        isZero = true;
        for(i=ii+1;i<nr;++i)
        {
            if(x[i][jj].Re || x[i][jj].Im)
            {
                isZero = false;
                k = i;
                break;
            }
        }
        
        if(isZero)
            return false;
        
        // have non-zero element in row k, so swap rows ii with k
        swapRows(x, nr, nc, ii, k);
    }
    // now we do the Gauss Elim. operation
    // pivot element xt[ii][jj] is not zero
    
    
    // now for every k not ii subtract row ii multiplied by xt[k][jj] from row k provided x[k][jj] != 0
    for(k=0;k<nr;++k)
    {
        if(k!=ii && (x[k][jj].Re || x[k][jj].Im))
        {
            //temp = x[k][jj] / x[ii][jj];
            DivCC(temp, x[k][jj], x[ii][jj]);
            for(j=0;j<nc;++j)
            {
                //x[k][j] = x[k][j] - (x[ii][j] * temp);
                MulCC(temp1, x[ii][j], temp);
                SubCC(x[k][j], x[k][j], temp1);
                if(dabs(x[k][j].Re)<zeroDouble)
                {
                    x[k][j].Re = 0.;
                }
                
                if(dabs(x[k][j].Im)<zeroDouble)
                {
                    x[k][j].Im = 0.;
                }
            }
        }
    }
    return true;
    
}/* doGauss */

void GaussElimC(Complex **z, long nr, long nc, Complex **x)
{
    Complex         **xt;
    long            i, j, ii, jj;
    
    xt = (Complex**)malloc(nr*sizeof(Complex*));
    for(i=0;i<nr;++i)
        xt[i] = (Complex*)malloc(nc*sizeof(Complex));
    
    EquateArrays(xt, nr, nc, x);
    
    ii = jj = 0;
    while(ii<nr && jj<nc)
    {
        if(doGauss(xt, nr, nc, ii, jj))
        {
            ii++;
            jj++;
        }
        else
            jj++;
    }
    
    for(i=0;i<nr ;++i)
        for(j=0;j<nc;++j)
        {
            if(!xt[i][j].Re)
                xt[i][j].Re = xt[i][j].Re + 0; // corrects -0
            
            if(!xt[i][j].Im)
                xt[i][j].Im = xt[i][j].Im + 0;
        }
    
    EquateArrays(z, nr, nc, xt);
    for(i=0;i<nr;++i)
        free(xt[i]);
    free(xt);
    
}/* GaussElimC */


// 1 when i=j, 0 when i!=j
// so, when nr=nc have unit matrix
// z must be inited by caller
void CreateUnit(Complex **z, long nr, long nc)
{
    long    i, j;
    
    for(i=0;i<nr;++i)
        for(j=0;j<nc;++j)
        {
            if(i==j)
            {
                z[i][j].Re = 1;
                z[i][j].Im = 0;
            }
            else
            {
                z[i][j].Re = 0;
                z[i][j].Im = 0;
            }
            
        }
    
}/* CreateUnit */


// 1 when i=j, 0 when i!=j
// so, when nr=nc have unit matrix
// z must be inited by caller
void CreateUnitD(double **z, long nr, long nc)
{
    long    i, j;
    
    for(i=0;i<nr;++i)
        for(j=0;j<nc;++j)
        {
            if(i==j)
            {
                z[i][j] = 1;
            }
            else
                z[i][j] = 0;
        }
    
}/* CreateUnitD */


void RealToComplex(Complex **z, long nr, long nc, double **x)
{
    long        i, j;
    
    for(i=0;i<nr;++i)
        for(j=0;j<nc;++j)
        {
            z[i][j].Re = x[i][j];
            z[i][j].Im = 0;
        }
    
}/* RealToComplex */

static bool doReduce(Complex **x, long nr, long nc, long ii, long jj)
{
    Complex         temp, temp1;
    long            i, j, k;
    bool            isZero;
    
    if(ii<0 || jj<0 || ii>=nr || jj>=nc)
        return false;
    
    k = 0; // to avoid incorrect warning
    // is pivot zero?
    if(!x[ii][jj].Re && !x[ii][jj].Im)
    {
        // must swap rows
        // find non-zero element below pivot in pivot column
        isZero = true;
        for(i=ii+1;i<nr;++i)
        {
            if(x[i][jj].Re || x[i][jj].Im)
            {
                isZero = false;
                k = i;
                break;
            }
        }
        
        if(isZero)
            return false;
        
        // have non-zero element in row k, so swap rows ii with k
        swapRows(x, nr, nc, ii, k);
    }
    // now we do the reduce operation
    // pivot element xt[ii][jj] is not zero
    // divide each element of row ii by xt[ii][jj] making xt[ii][jj] = 1
    temp.Re = x[ii][jj].Re;
    temp.Im = x[ii][jj].Im;
    for(j=0;j<nc;++j)
    {
        //x[ii][j] = x[ii][j] / temp;
        DivCC(x[ii][j], x[ii][j], temp);
    }
    
    // now for every k not ii subtract row ii multiplied by xt[k][jj] from row k provided x[k][jj] != 0
    for(k=0;k<nr;++k)
    {
        if(k!=ii && (x[k][jj].Re || x[k][jj].Im))
        {
            temp.Re = x[k][jj].Re;
            temp.Im = x[k][jj].Im;
            for(j=0;j<nc;++j)
            {
                //x[k][j] = x[k][j] - (x[ii][j] * temp);
                MulCC(temp1, x[ii][j], temp);
                SubCC(x[k][j], x[k][j], temp1);
                if(dabs(x[k][j].Re)<zeroDouble)
                {
                    x[k][j].Re = 0.;
                }
                
                if(dabs(x[k][j].Im)<zeroDouble)
                {
                    x[k][j].Im = 0.;
                }
            }
        }
    }
    return true;
    
}/* doReduce */

// z must be inited by caller
void Reduce(Complex **z, long nr, long nc, Complex **x)
{
    Complex          **xt;
    long            i, j, ii, jj;
    
    xt = (Complex**)malloc(nr*sizeof(Complex*));
    for(i=0;i<nr;++i)
        xt[i] = (Complex*)malloc(nc*sizeof(Complex));
    
    EquateArrays(xt, nr, nc, x);
    
    ii = jj = 0;
    while(ii<nr && jj<nc)
    {
        if(doReduce(xt, nr, nc, ii, jj))
        {
            ii++;
            jj++;
        }
        else
            jj++;
    }
    
    for(i=0;i<nr ;++i)
        for(j=0;j<nc;++j)
        {
            if(!xt[i][j].Re)
                xt[i][j].Re = xt[i][j].Re + 0; // corrects -0
            
            if(!xt[i][j].Im)
                xt[i][j].Im = xt[i][j].Im + 0;
        }
    
    EquateArrays(z, nr, nc, xt);
    for(j=0;j<nr;++j)
        free(xt[j]);
    free(xt);
    
}/* Reduce */

// z must be inited by caller
void ReduceD(double **z, long nr, long nc, double **x)
{
    double          **xt;
    long            i, ii, jj;
    
    xt = (double**)malloc(nr*sizeof(double*));
    for(i=0;i<nr;++i)
        xt[i] = (double*)malloc(nc*sizeof(double));
    
    EquateArraysD(xt, nr, nc, x);
    
    ii = jj = 0;
    while(ii<nr && jj<nc)
    {
        if(doReduceD(xt, nr, nc, ii, jj))
        {
            ii++;
            jj++;
        }
        else
            jj++;
    }
    
    EquateArraysD(z, nr, nc, xt);
    for(i=0;i<nr;++i)
        free(xt[i]);
    free(xt);
    
}/* ReduceD */

// augments y on right of x if xnr==ynr
void Augment(double **z, double **x, long xnr, long xnc, double **y, long ynr, long ync)
{
    double      **zt;
    long        i, j, znr, znc;
    
    if(xnr!=ynr)
        return;
    
    znr = xnr;
    znc = xnc + ync;
    
    zt = (double**)malloc(znr*sizeof(double*));
    for(i=0;i<znr;++i)
        zt[i] = (double*)malloc(znc*sizeof(double));
    
    for(i=0;i<znr;++i)
        for(j=0;j<znc;++j)
        {
            if(j<xnc)
            {
                zt[i][j] = x[i][j];
            }
            else
            {
                zt[i][j] = y[i][j-xnc];
            }
        }
    
    EquateArraysD(z, znr, znc, zt);
    for(j=0;j<znr;++j)
        free(zt[j]);
    free(zt);
    
}/* Augment */

// augments y on right of x if xnr==ynr
void Augment(Complex **z, Complex **x, long xnr, long xnc, double **y, long ynr, long ync)
{
    Complex      **zt;
    long         i, j, znr, znc;
    
    if(xnr!=ynr)
        return;
    
    znr = xnr;
    znc = xnc + ync;
    
    zt = (Complex**)malloc(znr*sizeof(Complex*));
    for(i=0;i<znr;++i)
        zt[i] = (Complex*)malloc(znc*sizeof(Complex));
    
    for(i=0;i<znr;++i)
        for(j=0;j<znc;++j)
        {
            if(j<xnc)
            {
                zt[i][j].Re = x[i][j].Re;
                zt[i][j].Im = x[i][j].Im;
            }
            else
            {
                zt[i][j].Re = y[i][j-xnc];
                zt[i][j].Im = 0;
            }
        }
    
    EquateArrays(z, znr, znc, zt);
    for(j=0;j<znr;++j)
        free(zt[j]);
    free(zt);
    
}/* Augment */

void Augment(Complex **z, double **x, long xnr, long xnc, Complex **y, long ynr, long ync)
{
    Complex      **zt;
    long         i, j, znr, znc;
    
    if(xnr!=ynr)
        return;
    
    znr = xnr;
    znc = xnc + ync;
    
    zt = (Complex**)malloc(znr*sizeof(Complex*));
    for(i=0;i<znr;++i)
        zt[i] = (Complex*)malloc(znc*sizeof(Complex));
    
    for(i=0;i<znr;++i)
        for(j=0;j<znc;++j)
        {
            if(j<xnc)
            {
                zt[i][j].Re = x[i][j];
                zt[i][j].Im = 0;
            }
            else
            {
                zt[i][j].Re = y[i][j-xnc].Re;
                zt[i][j].Im = y[i][j-xnc].Im;
            }
        }
    
    EquateArrays(z, znr, znc, zt);
    for(j=0;j<znr;++j)
        free(zt[j]);
    free(zt);
    
}/* Augment */


void Augment(Complex **z, Complex **x, long xnr, long xnc, Complex **y, long ynr, long ync)
{
    Complex      **zt;
    long         i, j, znr, znc;
    
    if(xnr!=ynr)
        return;
    
    znr = xnr;
    znc = xnc + ync;
    
    zt = (Complex**)malloc(znr*sizeof(Complex*));
    for(i=0;i<znr;++i)
        zt[i] = (Complex*)malloc(znc*sizeof(Complex));
    
    for(i=0;i<znr;++i)
        for(j=0;j<znc;++j)
        {
            if(j<xnc)
            {
                zt[i][j].Re = x[i][j].Re;
                zt[i][j].Im = x[i][j].Im;
            }
            else
            {
                zt[i][j].Re = y[i][j-xnc].Re;
                zt[i][j].Im = y[i][j-xnc].Im;
            }
        }
    
    EquateArrays(z, znr, znc, zt);
    for(j=0;j<znr;++j)
        free(zt[j]);
    free(zt);
    
}/* Augment */


// unAugment splits column-wise a matrix z into matrices x and y (z=xy) using xnc for split, xnc<znc
// neither x or y can be z
void unAugmentD(double **x, double **y, long xnc, long znr, long znc, double **z)
{
    long        i, j;
    
    if(xnc>=znc)
        return;
    
    for(i=0;i<znr;++i)
        for(j=0;j<xnc;++j)
            x[i][j] = z[i][j];
    
    for(i=0;i<znr;++i)
        for(j=xnc;j<znc;++j)
            y[i][j-xnc] = z[i][j];
    
}/* unAugmentD */

// unAugment splits column-wise a matrix z into matrices x and y (z=xy) using xnc for split, xnc<znc
// neither x or y can be z
void unAugmentC(Complex **x, Complex **y, long xnc, long znr, long znc, Complex **z)
{
    long        i, j;
    
    if(xnc>=znc)
        return;
    
    for(i=0;i<znr;++i)
        for(j=0;j<xnc;++j)
        {
            x[i][j].Re = z[i][j].Re;
            x[i][j].Im = z[i][j].Im;
        }
    
    
    for(i=0;i<znr;++i)
        for(j=xnc;j<znc;++j)
        {
            y[i][j-xnc].Re = z[i][j].Re;
            y[i][j-xnc].Im = z[i][j].Im;
        }
    
    
}/* unAugmentC */


