#ifndef MATRIX_H
#define MATRIX_H
/*
* @Author:  xuyong;
* mail: xyzw557@sohu.com
*/
#include <iostream>
using namespace std;
namespace _salt
{
	template <class T>
		class matrix{
		public :
			matrix(int row,int colum);
			matrix(const matrix& omatrix);
            ~matrix();
			friend istream& operator >> (istream& ins, matrix& omatrix);
		    friend ostream& operator << (ostream& outs,const matrix& omatrix);
			friend matrix operator+ (const matrix& matrix1,const matrix& matrix2);
			friend matrix operator- (const matrix& matrix1,const matrix& matrix2);
			friend matrix operator* (const matrix& matrix1,const matrix& matrix2);
			//friend matrix operator= ( const matrix& matrix2);
			friend int symbol(int m,int n);
			T get(int i,int j);
			void set(int i,int j,T elem);
			matrix mrest(int i,int j); 
            T hls();
            matrix trans();
			matrix adjoint();// adjoint matrix;
			matrix inverse_matrix();
		private:
			int mrow,mcolum;
			T **p;
	};
	template <class T>
		matrix<T>::matrix(int row,int colum)//create a matrix by given row and colum;
	{
		int i;
		mcolum=colum;
		mrow=row;
		p=new T* [row];
		for(i=0;i<row;i++)
           p[i]=new T[colum];
	}
	template <class T>
		matrix<T>::matrix(const matrix<T>& omatrix)//create a matrix by copying another matrix;
	{
		int i,j;
		mrow=omatrix.mrow;
		mcolum=omatrix.mcolum;
		p=new T* [mrow];
		for(i=0;i<mrow;i++)
			p[i]=new T[mcolum];
		for(i=0;i<mrow;i++)
			for(j=0;j<mcolum;j++)
				p[i][j]=omatrix.p[i][j];
	}
	template <class T>
		matrix<T>::~matrix()//destroy the matrix when the program exit;
	{
		int i;
		for(i=0;i<mrow;i++)
			delete [] p[i];
	}
	template <class T>
		T matrix<T>::get(int i,int j)//get element in matrix by given its row and colum;
	{
		try{
			if(i<0||i>mrow||j<0||j>mcolum) throw 0;
			return p[i][j];
		}
		catch(int n){
			cout<<"i,j are not reasonable!"<<endl;
		}
	}
	template <class T>
		void matrix<T>::set(int i,int j, T elem)//set element whose row and colum is given  in matrix 
	{
		try{
			if(i<0||i>mrow||j<0||j>mcolum) throw 0;
			p[i][j]=elem;
		}
		catch(int n){
			cout<<"i,j are not reasonable!"<<endl;
		}
	}
	/*template <class T>
		matrix<T> operator=(const matrix<T>& omatrix)
	{
		int i,j;
		mrow=omatrix.mrow;
		mcolum=omatrix.mcolum;
		p=new T* [mrow];
		for(i=0;i<mrow;i++)
			p[i]=new T[mcolum];
		for(i=0;i<mrow;i++)
			for(j=0;j<mcolum;j++)
				p[i][j]=omatrix.p[i][j];
		return *this;
	}*/

	template <class T>
		matrix<T> operator+(const matrix<T>& matrix1,const matrix<T>& matrix2)
	{
		//two matrixes add when they are resemble;
		try{
		    int i,j,c,r;
		    c=matrix1.mrow; r=matrix1.mcolum;
	     	matrix<T> temp(c,r);
			if(matrix1.mrow!=matrix2.mrow||matrix1.mcolum!=matrix2.mcolum) throw 0;
		    for(i=0;i<c;i++)
			    for(j=0;j<r;j++)
				   temp.p[i][j]=matrix1.p[i][j]+matrix2.p[i][j];
			return temp;
		}
		catch(int n){
			cout<<"Matrix cannot match!"<<endl;
		}
	}
	template <class T>
		matrix<T> operator-(const matrix<T>& matrix1,const matrix<T>& matrix2)
	{
		//two matrix minus ;
		try{
		    int i,j,c,r;
		    c=matrix1.mrow; r=matrix1.mcolum;
	     	matrix<T> temp(c,r);
			if(matrix1.mcolum!=matrix2.mcolum||matrix1.mrow!=matrix2.mrow) throw 0;
		    for(i=0;i<c;i++)
			    for(j=0;j<r;j++)
				   temp.p[i][j]=matrix1.p[i][j]-matrix2.p[i][j];
			return temp;
		}
		catch(int n){
			cout<<"Matrix cannot match!"<<endl;
		}
	}				
	template <class T>
		matrix<T> operator*(const matrix<T>& matrix1,const matrix<T>& matrix2)
	{
		//two matrixes multiple;
		try{
		    int i,j,k,c,r;
		    c=matrix1.mrow; r=matrix2.mcolum;
	     	matrix<T> temp(c,r);
			if(matrix1.mcolum!=matrix2.mrow) throw 0;
		    for(i=0;i<c;i++)
			    for(j=0;j<r;j++)
					for(k=0;k<matrix1.mcolum;k++){
						if(k==0){
							temp.p[i][j]=matrix1.p[i][0]*matrix2.p[0][j];
							continue;
						}
					   temp.p[i][j]+=matrix1.p[i][k]*matrix2.p[k][j];
					}
			return temp;
		}
		catch(int n){
			cout<<"Matrix cannot match!"<<endl;
		}
	}	

	template <class T>
		istream& operator>>(istream& ins,matrix<T>& mmatrix)
	{
		//reload opretator >> for input elements to matrix easily;
		int i,j;
		for(i=0;i<mmatrix.mrow;i++)
		{
			for(j=0;j<mmatrix.mcolum;j++)
				ins>>mmatrix.p[i][j];
		}

		return ins;
	}
	template <class T>
		ostream& operator<<(ostream& outs,const matrix<T>& mmatrix)
	{
		//reload operator << for output matrix;
		int i,j;
		for(i=0;i<mmatrix.mrow;i++)
			for(j=0;j<mmatrix.mcolum;j++){
				if(j%mmatrix.mcolum==0) outs<<endl;
				outs<<mmatrix.p[i][j]<<"    ";
			}
			outs<<"\n";
			return outs;
	}

	template<class T>
		T matrix<T>::hls()
	{
		//the determinant of matrix;
		int i;
		T val;
		try{
			if(mcolum!=mrow) throw 0;
			if(mrow==1) return p[0][0];
			for(i=0;i<mrow;i++)
			{
				matrix<T> temp1(mrest(0,i));
				if(i==0) {
					
					val=temp1.hls();
					val*=p[0][i];
				}
				else val+=symbol(0,i)*p[0][i]*temp1.hls();
			}
			return val;
		}
		catch(int n){
			cout<<"Matrix cannot match!"<<endl;
			exit(1);
		}
	}
	int symbol(int m,int n)
	{// get value of (-1)^(m+n ) ;
		int val,i;
		val=1;
		for(i=0;i<m+n;i++)
			val*=-1;
		return val;
	}
	template <class T>
		matrix<T> matrix<T>::mrest(int i,int j)
	{  //cofactor of determinant;
		int r,s,k,t;
		matrix<T> temp(mrow-1,mcolum-1);
		for(k=0, r=0;r<mrow;r++)
		{
			if(r==i) continue;
			for(t=0,s=0;s<mcolum;s++)
			{
				if(s==j) continue;
				temp.p[k][t]=p[r][s];
				t++;
			}
			k++;
		}
		return temp;
	}
	template <class T>
		matrix<T> matrix<T>::trans()
	{
		//transpose matrix;
		int i,j;
		matrix<T> temp(mcolum,mrow);
		for(i=0;i<mrow;i++)
		{
			for(j=0;j<mcolum;j++)
				temp.p[j][i]=p[i][j];
		}
		return temp;
	}
		template <class T>
		matrix<T> matrix<T>::adjoint()
	{
		//adjoint matrix;
			int i,j;
			T det;
			matrix<T> temp(mcolum,mrow);
			try{
				if(mcolum!=mrow) throw 0;
				for(i=0;i<mrow;i++)
					for(j=0;j<mcolum;j++)
						temp.p[i][j]=symbol(i,j)*mrest(i,j).hls();
				return temp.trans();
			}
			catch(int n){
				cout<<"Matrix cannot match!"<<endl;
				exit(1);
			}
	}
		template <class T>
		matrix<T> matrix<T>::inverse_matrix()
	{
			int i,j;
			T det;
			det=hls();
			matrix<T> temp(mcolum,mrow);
			try{
				if (det==0) throw 0;
				for (i=0;i<mrow;i++)
					for(j=0;j<mcolum;j++)
						temp.p[i][j]=adjoint().p[i][j]/det;
				return temp;
			}
			catch(int n){
				cout<<"Determinant of Matrix equals zero!"<<endl;
				exit(1);
			}
	}

}

#endif