#include <stdio.h>
#include <math.h>
#include <memory.h>
#include "Matrix.h"


matrix::matrix(int m, int n, double* d ) :
 M(m), N(n)
{
	X = new double [M*N];
	if(d){
 		memcpy(X, d, M*N*sizeof(double) );
   }else
		memset(X,  0, M*N*sizeof(double) );
};


matrix::matrix(int m, int n, unsigned int * data ):
	M(m),N(n)
{
	X=new double [M*N];
	double* u = X;
	unsigned int* v=data;
	if(data){
		for (int i = 0; i < M*N; i++)	*u++ = (double) *v++;
	}else{
		memset(X,  0, M*N*sizeof(double) );
	}
};

matrix::matrix(int m, int n, unsigned char * data ):
	M(m),N(n)
{
	X=new double [M*N];
	double* u = X;
	unsigned char* v=data;
	if(data){
		for (int i = 0; i < M*N; i++)	*u++ = (double) *v++;
	}else{
		memset(X,  0, M*N*sizeof(double) );
	}
};


matrix operator* (const matrix& a, const double& b)
{
	int N = a.n(); int M = a.m();
	double *d = new double [N*M];
	double *u = d; double *v = a.x();
	for (int y = 0; y < M; y++)  {
			for (int x = 0; x < N; x++) *u++ = *v++ * b;
	};
	return matrix(M, N, d);
};

matrix operator* (const matrix& a, const matrix& b)
{
	if (a.n() != b.m())  {
		fprintf(stderr, "illegal matrix * operation.\n");
		return matrix();
	};
	int M = a.m(); int N = b.n(); int K = a.n();
	double *d = new double [M*N]; double *u = d;
	for (int y = 0; y < M; y++)  {
			for (int x = 0; x < N; x++) {
			register double sum = 0.0;
			register double *ai = a.x()+y*K;
			register double *bi = b.x()+x;
			for (int k = 0; k < K; k++, bi += N) sum += *ai++ * *bi;
			*u++ = sum;
		};
	};
	return matrix(M, N, d);
}

matrix operator/ (const matrix& a, const double& b)
{
	if (b == 0.0)  {
		fprintf(stderr, "matrix / with 0 operand.");
		return matrix();
	};
	int M = a.m(); int N = a.n();
	double *d = new double [M*N];
	double *u = d; double *v = a.x();
	for (int i = 0; i < M; i++)
			for (int j = 0; j < N; j++) *u++ = *v++ / b;
	return matrix(M, N, d);
}

matrix& matrix::operator= (matrix const &a)
{
	if(X) delete [] X;
	N = a.n(); M = a.m();
	X = new double [N*M];
	memcpy(X, a.x(), N*M*sizeof(double));
	return *this;
}

matrix& matrix::operator= (int const d)
{
	double *u = X;
	for (int i = 0; i < M*N; i++) *u++ = (double)d;
	return *this;
}

matrix& matrix::operator= (int* const d)
{
	double *u = X;
	int *v=d;
	for (int i = 0; i < M*N; i++) *u++ = (double) *v++;
	return *this;
}


matrix& matrix::operator= (double const d)
{
	double *u = X;
	for (int i = 0; i < M*N; i++) *u++ = d;
	return *this;
}


int operator== (const matrix& a, const matrix& b)
{
	if (a.n() != b.n()) return 0;
	if (a.m() != b.m()) return 0;
	double *u = a.x(); double *v = b.x();
	for (int i = 0; i < a.m(); i++)  {
			for (int j = 0; j < a.n(); j++) if (*u++ != *v++) return 0;
	};
	return 1;
}

matrix operator+ (const matrix& a, const double& b)
{
	if (b == 0.0)  {
		fprintf(stderr, "matrix / with 0 operand.");
		return matrix();
	};
	int M = a.m(); int N = a.n();
	double *d = new double [M*N];
	double *u = d; double *v = a.x();
	for (int i = 0; i < M; i++)
			for (int j = 0; j < N; j++) *u++ = *v++ + b;
	return matrix(M, N, d);
}

matrix operator+ (const matrix& a, const matrix& b)
{
	if ((a.n() != b.n()) || (a.m() != b.m()))  {
		fprintf(stderr, "illegal matrix + operation.\n");
		return matrix();
	};
	int N = a.n(); int M = a.m();
	double *d = new double [N*M];
	double *u = d; double *v = a.x(); double *w = b.x();
	for (int y = 0; y < M; y++)  {
			for (int x = 0; x < N; x++) *u++ = *v++ + *w++;
	};
	return matrix(M, N, d);
}

matrix operator- (const matrix& a, const double& b)
{
	if (b == 0.0)  {
		fprintf(stderr, "matrix / with 0 operand.");
		return matrix();
	};
	int M = a.m(); int N = a.n();
	double *d = new double [M*N];
	double *u = d; double *v = a.x();
	for (int i = 0; i < M; i++)
			for (int j = 0; j < N; j++) *u++ = *v++ - b;
	return matrix(M, N, d);
}

matrix operator- (const matrix& a, const matrix& b)
{
	if ((a.n() != b.n()) || (a.m() != b.m()))  {
		fprintf(stderr, "illegal matrix + operation.\n");
		return matrix();
	};
	int N = a.n(); int M = a.m();
	double *d = new double [N*M];
	double *u = d; double *v = a.x(); double *w = b.x();
	for (int y = 0; y < M; y++)  {
			for (int x = 0; x < N; x++) *u++ = *v++ - *w++;
	};
	return matrix(M, N, d);
}

matrix matrix::c(int n) const			// extract column
{
	if (n >= N)  {
		fprintf(stderr, "illegal column extraction.\n");
		return matrix();
	};
	double *d = new double [M]; double *u = d; double *v = X+n;
	for (int y = 0; y < M; y++, v += N) *u++ = *v;
	return matrix(M, 1, d);
}

matrix matrix::sum() const
{
	double* d = new double [M];
	for (int i = 0; i < M; i++)  {
		double sum = 0; double* u = X+i*N;
		for (int j = 0; j < N; j++) sum += *u++;
		d[i] = sum;
	};
	return matrix(1, M, d);
}

matrix matrix::sumsq() const
{
	double* d = new double [M];
	for (int i = 0; i < M; i++)  {
		double sum = 0; double* u = X+i*N;
		for (int j = 0; j < N; j++, u++) sum += *u * *u;
		d[i] = sum;
	};
	return matrix(1, M, d);
}


matrix matrix::min() const
{
	double* d = new double [M];
	for (int i = 0; i < M; i++)  {
		double* u = X+i*N; double Min = *u;
		for (int j = 0; j < N; j++, u++) if (Min > *u) Min = *u;
		d[i] = Min;
	};
	return matrix(1, M, d);
}

double matrix::max() const
{
	double* d = new double [M];
        double Max;
	for (int i = 0; i < M; i++)  {
		double* u = X+i*N;
                Max = *u;
		for (int j = 0; j < N; j++, u++) if (Max < *u) Max = *u;
		d[i] = Max;
	};
	return Max;
}


void matrix::identity() const
{
	for (int i = 0; i < M; i++)
		for (int j = 0; j < N; j++)
			X[i*N+j] = (i == j)? 1.0: 0.0;
}


matrix::matrix(int m, int n, const double& a)
{
	M = m; N = n; X = new double [m*n]; double* t = X;
	for (int i = 0; i < m*n; i++) *t++ = a;
}