/* ProjectionM.cp */

#include <iostream>
#include "ProjectionM.h"
#include "Complex.h"
#include "cmatinv.h"
#include <stdlib.h>

using namespace std;

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++)
		{
			z[i][j].Re = x[i][j].Re;
			z[i][j].Im = x[i][j].Im;
		}
	
}/* equateM_M */


// z = x*y
static inline void MulM_M(Complex **z, Complex **x, Complex **y, Complex **t, long 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 */


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


static inline void assignZeroM(Complex **z, long 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 */


// give this the k distinct root[], multiplicity m[], and have it calc k projection polys P
bool projectionPoly(Polynomial *P, Complex *root, INT32 *m, INT32 k, INT32 n)
{
	Polynomial		*a, *b;
	long			i, j, jj, rj;
	Polynomial		linear;
	Complex			**M, det;
	bool			theReturn;
	
	a = (Polynomial*)malloc(k*sizeof(Polynomial));
	b = (Polynomial*)malloc(k*sizeof(Polynomial));
	
	M = (Complex**)malloc(n*sizeof(Complex*));
	for(i=0;i<n;i++)
		M[i] = (Complex*)malloc(n*sizeof(Complex));
		
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
	{
		M[i][j].Re = 0;
		M[i][j].Im = 0;
	}
	
	CreateP(linear, 1);
	linear.coeff[1].Re = 1;
	linear.coeff[1].Im = 0;
	
	for(i=0;i<k;i++)
	{
		CreateP(a[i], m[i]-1);
		CreateP(b[i], n-m[i]);
		// init b[i] as 1
		b[i].coeff[0].Re = 1;
	}
	
	// calculate the b[]
	for(i=0;i<k;i++)
		for(j=0;j<k;j++)
		{
			if(j!=i)
			{
				linear.coeff[0].Re = -root[j].Re;
				linear.coeff[0].Im = -root[j].Im;
				
				for(jj=0;jj<m[j];jj++)
					MulP_P(b[i], b[i], linear);
			}
		}
		
	// print out the b[]
	/*
	for(i=0;i<k;i++)
	{
		cout << endl << "b[" << i << "]:" << endl;
		for(j=0;j<=b[i].n;j++)
			cout << "coeff[" << j << "] = " << b[i].coeff[j].Re << " , " << b[i].coeff[j].Im << endl;
			
	}
	*/
		
	// construct the matrix M
	j = 0;
	for(rj=0;rj<k;rj++)
	{
		for(jj=0;jj<m[rj];jj++)
		{
			for(i=jj;i<=n-m[rj]+jj;i++)
				Add(M[i][j+jj], M[i][j+jj], b[rj].coeff[i-jj]);
		}
		j+=m[rj];
	}
	
	// print out M
	/*
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			cout << M[i][j].Re << " , " << M[i][j].Im << endl;
	*/
	
	// invert the matrix
	det = cmatinv(M, n);
	if(isEqualZero(det))
	{
		theReturn = false;
		goto endgame;
	}
		
	// the unknowns in a[i] are the zeroth column of M, M[j][0]
	j = 0;
	for(rj=0;rj<k;rj++)
	{
		for(jj=0;jj<m[rj];jj++)
		{
			EquateC_C(a[rj].coeff[jj], M[j+jj][0]);
		}
		j+=m[rj];
	}
	
	// now form the projection polys
	for(i=0;i<k;i++)
		MulP_P(P[i], a[i], b[i]);
		
	theReturn = true;
		
	endgame:
	
	for(i=0;i<k;i++)
	{
		FreeP(a[i]);
		FreeP(b[i]);
	}
	
	free(a);
	free(b);
	
	FreeP(linear);
	
	for(i=0;i<n;i++)
		free(M[i]);
		
	free(M);
	
	return theReturn;

}/* projectionPoly */


void projMatrices(Complex ***A, Complex **M, Polynomial *P, INT32 k, INT32 n)
{
	long		i, j;
	Complex		**t;
	
	t = (Complex**)malloc(n*sizeof(Complex*));
	for(i=0;i<n;i++)
		t[i] = (Complex*)malloc(n*sizeof(Complex));
	
	// all P[i] have highest power n-1
	
	/*
	EquateC_C(P, c[n]);
	
	for(i=n-1;i>=0;i--)
	{
		MulC_C(P, P, x);
		AddC_C(P, P, c[i]);
	}
	*/
	
	for(i=0;i<k;i++)
	{
		/*
		equateM_M(A[i], M, n);
		MulM_C(A[i], P[i].coeff[n-1], n);
		*/
		
		assignZeroM(A[i], n);
		AddM_C(A[i], P[i].coeff[n-1], n);
		
		for(j=n-2;j>=0;j--)
		{
			MulM_M(A[i], A[i], M, t, n);
			AddM_C(A[i], P[i].coeff[j], n);
		}
	}
	
	for(i=0;i<n;i++)
		free(t[i]);
		
	free(t);

}/* projMatrices */

