/* charPoly.cp */

#include "Complex.h"
#include "charPolyOld.h"
#include <stdlib.h>
#include <math.h>

/*
	Finds the characteristic polynomial of a square matrix using a method described by
	Prof. Bumby of Rutgers Univ.
	
	Let M be the matrix of size nxn, and the characteristic polynomial be:
	
	x^n - p[1]*x^n-1 - p[2]*x^n-2 - ... - p[n], and let p[0] = 1
	
	Then define matrices F such that
	
	F_1(M) = M
	
	and for k = 1, 2,...,n
	
	p[k] = (1/k)*Tr(F_k)
	
	F_k+1 = M*(F_k - p[k]*I)
	
*/

/*
	Calling program should do:
	
	M = (Complex**)malloc(n*sizeof(Complex*));
	for(i=0;i<n;i++)
		M[i] = (Complex*)malloc(n*sizeof(Complex));
		
	p = (Complex*)malloc((n+1)*sizeof(Complex));
	
	and free memory when done
*/

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

}/* MulM_M */


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


static inline Complex TraceM(Complex **x, long n)
{
	long		i;
	Complex		result;
	
	assignZeroC(result);
	
	for(i=0;i<n;i++)
		Add(result, result, x[i][i]);
	
	return result;
	
}/* TraceM */



void CharPolyOld(Complex **M, Complex *p, long n)
{
	Complex		**temp, **F, temp1;
	long		i, k;
	
	temp = (Complex**)malloc(n*sizeof(Complex*));
	for(i=0;i<n;i++)
		temp[i] = (Complex*)malloc(n*sizeof(Complex));
		
	F = (Complex**)malloc(n*sizeof(Complex*));
	for(i=0;i<n;i++)
		F[i] = (Complex*)malloc(n*sizeof(Complex));
		
	p[0].Re = 1;
	p[0].Im = 0;
	
	EquateM_M(F, M, n);
	
	for(k=1;k<n;k++)
	{
		EquateC_C(temp1, TraceM(F, n));
		p[k].Re = temp1.Re/k;
		p[k].Im = temp1.Im/k;
		
		for(i=0;i<n;i++)
			Sub(F[i][i], F[i][i], p[k]);
			
		MulM_M(temp, M, F, n);
		EquateM_M(F, temp, n);
	}
	
	EquateC_C(temp1, TraceM(F, n));
	p[n].Re = temp1.Re/n;
	p[n].Im = temp1.Im/n;
	
	// change signs except for p[0]
	for(i=1;i<=n;i++)
	{
		p[i].Re = -p[i].Re;
		p[i].Im = -p[i].Im;
	}
	
	// now reverse p[] to get coefficients in standard form
	for(i=0;i<=n/2;i++)
	{
		EquateC_C(temp1, p[i]);
		EquateC_C(p[i], p[n-i]);
		EquateC_C(p[n-i], temp1);
	}
	
	for(i=0;i<n;i++)
		free(temp[i]);
		
	free(temp);
	
	for(i=0;i<n;i++)
		free(F[i]);
		
	free(F);
	
}/* CharPoly */
