/* ----------------------------------------------------------------- */
/* 'ZAlgebraClass.cc'									    */
/* Algebra- Class			 						         */
/* 													    */
/* 	last modified 07/04/2003								    */
/* 	provides functions and structures for the calculations		    */
/* 	of mathematic things needed in linear algebra			    */
/* 													    */
/* 													    */
/* 													    */
/* 	(C) Jochen Zimmermann								    */
/* ----------------------------------------------------------------- */

#include "ZHeaders.h"
/*
//Headerfile is included by ZHeaders.h
#include "ZAlgebraClass.h"
*/

//constructor
//initialises members
LinearAlgebra::LinearAlgebra() 
{ 
	int i,j;
	rows=columns=0;
	for(i=0;i<4;i++)
		for(j=0;j<4;j++)
			matrix[i][j]=0.0;
}


//non-standard constructor for initialisation-during-creation
LinearAlgebra::LinearAlgebra(double* data, int wrows, int wcolumns)
{
	int i,j;
	rows=wrows;
	columns=wcolumns;
	for(i=0;i<rows;i++)
		for(j=0;j<columns;j++)
			matrix[i][j]=data[i*columns+j];
}


//copy constructor
LinearAlgebra::LinearAlgebra(const LinearAlgebra& tocopy)
{
	int i,j;
	rows=tocopy.rows;
	columns=tocopy.columns;
	for(i=0;i<rows;i++)
		for(j=0;j<columns;j++)
			matrix[i][j]=tocopy.matrix[i][j];
}


//operator '*' overloaded for scalar multiplications
LinearAlgebra LinearAlgebra::operator * (double scalar)
{
	int i,j;
	LinearAlgebra temp=*this;
	for(i=0;i<rows;i++)
		for(j=0;j<columns;j++)
			temp.matrix[i][j]*=scalar;			
	return temp;			
}


//operator '*=' overloaded for scalar multiplications
LinearAlgebra LinearAlgebra::operator *=(double scalar)
{
	int i,j;
	for(i=0;i<rows;i++)
		for(j=0;j<columns;j++)
			matrix[i][j]*=scalar;			
	return *this;			
}


//operator '*' overloaded for matrix-matrix and matrix-vector multiplications
LinearAlgebra LinearAlgebra::operator * (const LinearAlgebra& lag)
{
	int i,j,k;
	LinearAlgebra temp;
	//matrices not matching
	if(columns!=lag.rows) return temp;
	
	//scalar product
	if(columns==1 && lag.columns==1 && rows==lag.rows)
	{
		temp.matrix[0][0]=0;
		temp.rows=rows;
		temp.columns=columns;
		for(i=0;i<rows;i++)
			temp.matrix[0][0]+=matrix[i][0]*lag.matrix[i][0];
		return temp;
	}
		
	//multiply matrices/vector
	for (i=0;i<rows;i++)
		for(j=0;j<lag.columns;j++)
			for (k=0;k<columns;k++)
				(temp.matrix[i][j])+=matrix[i][k]*(lag.matrix[k][j]);
	
	temp.rows=rows;
	temp.columns=lag.columns;
	return temp;
}


//operator '*=' overloaded for matrix-matrix and matrix-vector multiplications
LinearAlgebra LinearAlgebra::operator *=(const LinearAlgebra& lag)
{
	int i,j,k;
	LinearAlgebra temp;
	
	//matrices not matching
	if(columns!=lag.rows) return temp;
	
	//scalar product
	if(columns==1 && lag.columns==1 && rows==lag.rows)
	{
		temp.matrix[0][0]=0;
		temp.rows=rows;
		temp.columns=columns;
		for(i=0;i<rows;i++)
			temp.matrix[0][0]+=matrix[i][0]*lag.matrix[i][0];
		*this=temp;
		return *this;
	}

	//multiply matrices
	for (i=0;i<rows;i++)
		for(j=0;j<lag.columns;j++)
			for (k=0;k<columns;k++)
				(temp.matrix[i][j])+=matrix[i][k]*(lag.matrix[k][j]);
	
	temp.rows=rows;
	temp.columns=lag.columns;
	*this=temp;
	
	return *this;
}


//method for private-member-retrieving
void LinearAlgebra::GetData(double* data,int wrows,int wcolumns,int firstrow, int firstcolumn)
{
	int i,j;
	for(i=0;i<wrows;i++)
		for(j=0;j<wcolumns;j++)
			data[i*wcolumns+j]=matrix[i+firstrow-1][j+firstcolumn-1];
}


//operator '+' overloaded for matrix-matrix and vector-vector additions
LinearAlgebra LinearAlgebra::operator + (const LinearAlgebra& lag)
{
	int i,j;
	LinearAlgebra temp;
	
	if(lag.rows!=rows && lag.columns!=columns) return temp;
	else temp=*this;
	for(i=0;i<temp.rows;i++)
		for(j=0;j<temp.columns;j++)
			temp.matrix[i][j]+=lag.matrix[i][j];
	return temp;
}


//operator '+=' overloaded for matrix-matrix and vector-vector additions
LinearAlgebra LinearAlgebra::operator +=(const LinearAlgebra& lag)
{
	int i,j;
	LinearAlgebra temp;
	if(lag.rows!=rows && lag.columns!=columns) return temp;
	else temp=*this;
	for(i=0;i<temp.rows;i++)
		for(j=0;j<temp.columns;j++)
			temp.matrix[i][j]+=lag.matrix[i][j];
	
	*this=temp;			
	return *this;
}


//operator '-' overloaded for matrix-matrix and vector-vector subtractions
LinearAlgebra LinearAlgebra::operator - (const LinearAlgebra& lag)
{
	int i,j;
	LinearAlgebra temp;
	
	if(lag.rows!=rows && lag.columns!=columns) return temp;
	else temp=*this;
	for(i=0;i<temp.rows;i++)
		for(j=0;j<temp.columns;j++)
			temp.matrix[i][j]=temp.matrix[i][j]-lag.matrix[i][j];
	return temp;
}


//operator '-=' overloaded for matrix-matrix and vector-vector subtractions
LinearAlgebra LinearAlgebra::operator -=(const LinearAlgebra& lag)
{
	int i,j;
	LinearAlgebra temp;
	if(lag.rows!=rows && lag.columns!=columns) return temp;
	else temp=*this;
	for(i=0;i<temp.rows;i++)
		for(j=0;j<temp.columns;j++)
			temp.matrix[i][j]=temp.matrix[i][j]-lag.matrix[i][j];	
	*this=temp;			
	return *this;
}


//method for assigning member values (similar operator '=')
LinearAlgebra LinearAlgebra::Equals(double* data, int wrows, int wcolumns)
{
	int i,j;
	for(i=0;i<wrows;i++)
		for(j=0;j<wcolumns;j++)
			matrix[i][j]=data[i*wcolumns+j];
	rows=wrows;
	columns=wcolumns;
	return *this;
}


//method for addition (similar Operator '+')
LinearAlgebra LinearAlgebra::Plus(double* data, int wrows, int wcolumns)
{
	int i,j;
	LinearAlgebra temp;
	
	if(wrows!=rows && wcolumns!=columns) return temp;
	else temp.Equals(data,wrows,wcolumns);
	for(i=0;i<temp.rows;i++)
		for(j=0;i<temp.columns;j++)
			temp.matrix[i][j]+=matrix[i][j];
	return temp;
}


//method to determine if actual instance is a matrix or not
bool  LinearAlgebra::IsMatrix()
{
	if(columns>1) return true;
	else return false;
}


//method to determine if actual instance is a vector or not
bool  LinearAlgebra::IsVector()
{
	if(columns==1) return true;
	else return false;
}


//method for retrieving member 'rows'
int   LinearAlgebra::GetRows()
{
	return rows;
}


//method for retrieving member 'columns'
int   LinearAlgebra::GetColumns()
{
	return columns;
}


//method to calculate absolute if instance is a vector or determinant if instance is a square-matrix
double LinearAlgebra::Absolute()
{
	double val=0;
	int i;
	if(columns==1)
	{
		for(i=0;i<rows;i++)
			val+=(matrix[i][0]*matrix[i][0]);
		val=sqrt(val);
	}
	else val=Det();	
	return val;
}

//method to calculate determinant of a square-matrix
double LinearAlgebra::Det()
{
	double val1=0,val2=1.0,val3=1.0;
	int i,j,k,l;
	LinearAlgebra temp;
	if(columns!=rows) return 0;
	else switch(columns)
	{
	case 1:	return matrix[0][0];
	case 2:	val1=(matrix[0][0]*matrix[1][1]-matrix[1][0]*matrix[0][1]);
			return val1;
			break;
	case 3:	//algorithm of Sarrus
			for(i=0;i<3;i++)
			{
				val2=val3=1.0;
				for(j=0;j<3;j++)
				{
					val2*=matrix[j][(i+j)%3];
					val3*=matrix[j][(i+2-j)%3];
				}
				val1+=(val2-val3);
			}
			return val1;
			break;
	case 4:	//development algorithm of Laplace
			for(i=0;i<4;i++)
			{
				val2=1.0;
				for(j=0;j<=i;j++)
					val2*=(-1.0);
					
				for(k=0;k<3;k++)
					for(l=0;l<3;l++)
						if(l>=i)	temp.matrix[k][l]=matrix[k][l+1];
						else 	temp.matrix[k][l]=matrix[k][l];
				temp.rows=temp.columns=3;
				
				if(matrix[3][i]!=0)
					val1+=(val2*matrix[3][i]*temp.Det());			
			}
			return val1;
			break;
	default:	return 0;
	}
}


/*
//Square-root lookup-method (not used, because of bad performance)
//disconnect this OwnSQRT- function and the lookup- table to increase speed (and size of .hex- file)
int   LinearAlgebra::OwnSqrt(double x)
{
	int i=0;
	if(x>=64620 || x<0) return -1;
	while(x>=(sqrtlookup[i]*2.0+0.25) )i++;
	return i;
}
//Square-root lookup-table (not used, because of bad performance)
//disconnect this lookup- table and the OwnSQRT- function to increase speed (and size of .hex- file)
const unsigned short LinearAlgebra::sqrtlookup[360]={
    0,    1,    3,    6,   10,   15,   21,   28,   36,   45,   55,   66,   78,   91,  105,  120,
  136,  153,  171,  190,  210,  231,  253,  276,  300,  325,  351,  378,  406,  435,  465,  496,
  528,  561,  595,  630,  666,  703,  741,  780,  820,  861,  903,  946,  990, 1035, 1081, 1128,
 1176, 1225, 1275, 1326, 1378, 1431, 1485, 1540, 1596, 1653, 1711, 1770, 1830, 1891, 1953, 2016,
 2080, 2145, 2211, 2278, 2346, 2415, 2485, 2556, 2628, 2701, 2775, 2850, 2926, 3003, 3081, 3160,
 3240, 3321, 3403, 3486, 3570, 3655, 3741, 3828, 3916, 4005, 4095, 4186, 4278, 4371, 4465, 4560,
 4656, 4753, 4851, 4950, 5050, 5151, 5253, 5356, 5460, 5565, 5671, 5778, 5886, 5995, 6105, 6216,
 6328, 6441, 6555, 6670, 6786, 6903, 7021, 7140, 7260, 7381, 7503, 7626, 7750, 7875, 8001, 8128,
 8256, 8385, 8515, 8646, 8778, 8911, 9045, 9180, 9316, 9453, 9591, 9730, 9870,10011,10153,10296,
10440,10585,10731,10878,11026,11175,11325,11476,11628,11781,11935,12090,12246,12403,12561,12720,
12880,13041,13203,13366,13530,13695,13861,14028,14196,14365,14535,14706,14878,15051,15225,15400,
15576,15753,15931,16110,16290,16471,16653,16836,17020,17205,17391,17578,17766,17955,18145,18336,
18528,18721,18915,19110,19306,19503,19701,19900,20100,20301,20503,20706,20910,21115,21321,21528,
21736,21945,22155,22366,22578,22791,23005,23220,23436,23653,23871,24090,24310,24531,24753,24976,
25200,25425,25651,25878,26106,26335,26565,26796,27028,27261,27495,27730,27966,28203,28441,28680,
28920,29161,29403,29646,29890,30135,30381,30628,30876,31125,31375,31626,31878,32131,32385,32640,
32896,33153,33411,33670,33930,34191,34453,34716,34980,35245,35511,35778,36046,36315,36585,36856,
37128,37401,37675,37950,38226,38503,38781,39060,39340,39621,39903,40186,40470,40755,41041,41328,
41616,41905,42195,42486,42778,43071,43365,43660,43956,44253,44551,44850,45150,45451,45753,46056,
46360,46665,46971,47278,47586,47895,48205,48516,48828,49141,49455,49770,50086,50403,50721,51040,
51360,51681,52003,52326,52650,52975,53301,53628,53956,54285,54615,54946,55278,55611,55945,56280,
56616,56953,57291,57630,57970,58311,58653,58996,59340,59685,60031,60378,60726,61075,61425,61776,
62128,62481,62835,63190,63546,63903,64261,64620};
*/