#ifndef __MATRIX_H__
#define __MATRIX_H__

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <string>
#include <exception>

#define  printf  ((int (*)(const char *,...))0x00002bb0)
#define  scanf   ((int (*)(const char *,...))0x00002c02)

template <class T>
class Matrix{
  private:
    T *m_Values;
    int *ref;				//参照カウンタ
    unsigned int m_Rows;
    unsigned int m_Columns;
  public:
    /**
     * Matrixクラスのコンストラクタ。
     * 指定の行数、指定の列数で行列を生成します。
     * また成分はすべてT(0)で初期化されます。
     * 
     * @param rows 行数
     * @param columns 列数
     */
    Matrix(const unsigned int &rows, const unsigned int &columns) : m_Rows(rows), m_Columns(columns){
    	m_Values = new T[rows * columns];
    	ref = new int(0);
      if((!m_Values) || (!ref)){}
      (*ref)++;
      for(unsigned int i = 0; i < rows * columns; i++){*(m_Values + i) = T(0);}
    }
    /**
     * コピーコンストラクタ。
     * シャローコピーを生成します。
     * 
     * @param matrix コピー元
     */
    Matrix(const Matrix &matrix){
      if(m_Values = matrix.m_Values){
        (*(ref = matrix.ref))++;
        m_Rows = matrix.m_Rows;
        m_Columns = matrix.m_Columns;
      }
    }
    /**
     * デストラクタ。
     */
    ~Matrix(){
      if(m_Values && ((--(*ref)) <= 0)){delete [] m_Values; delete ref;}
     }
    
    /**
     * 指定の単位行列を生成します。
     * 
     * @parma size 指定の行数(列数)
     */
    static Matrix<T> getI(const int &size){
      Matrix<T> result(size, size);
      for(int i = 0; i < size; i++){result(i, i) = T(1);}
      return result;
    }
    
    /**
     * 行数を返します。
     * 
     * @return (int) 行数
     */
    unsigned int rows() const{return m_Rows;}
    /**
     * 列数を返します。
     * 
     * @return (int) 列数
     */
    unsigned int columns() const{return m_Columns;}
    
    /**
     * 指定した行の行ベクトルを生成します。
     * 
     * @param row 行インデックス
     * @return (Matrix<T>) 行ベクトル
     */
    Matrix<T> rowVector(const unsigned int &row) const {
      if(row >= m_Rows){}
      Matrix<T> vector(1, m_Columns);
      for(int j = 0; j < m_Columns; j++){vector(0, j) = (*const_cast<Matrix<T> *>(this))(row, j);}
      return vector;
    }
    /**
     * 指定した列の列ベクトルを生成します。
     * 
     * @param column 列インデックス
     * @return (Matrix<T>) 列ベクトル
     */
    Matrix<T> columnVector(const unsigned int &column) const {
      if(column >= m_Columns){}
      Matrix<T> vector(m_Rows, 1);
      for(int i = 0; i < m_Rows; i++){vector(i, 0) = (*const_cast<Matrix<T> *>(this))(i, column);}
      return vector;
    }
    
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    T &operator()(const unsigned int &row, const unsigned int &column) {
      if((row < 0 || row >= m_Rows) || (column < 0 || column >= m_Columns)){}
      return *(m_Values + (row * m_Columns) + column);
    }
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行番号(開始番号は1〜)
     * @param column 列番号(開始番号は1〜)
     * @return (T) 成分
     */
    T &index(const unsigned int &row, const unsigned int &column) {
      return (*this)(row - 1, column - 1);
    }
    
    /**
     * 行列を複製(ディープコピー)します。
     * 
     * @return (Matrix<T>) コピー
     */
    Matrix<T> copy() const{
      Matrix<T> matrix(m_Rows, m_Columns);
      if(matrix.m_Values){
      	memcpy(matrix.m_Values, m_Values, sizeof(T) * m_Rows * m_Columns);
      }
      return matrix;
    }
    
    /**
     * 行を入れ替えます。破壊的メソッドです。
     * 
     * @param row1 行インデックス1
     * @param row2 行インデックス2
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &exchangeRows(const unsigned int &row1, const unsigned int &row2) {
      if(row1 >= m_Rows || row2 >= m_Rows){}
      T temp;
      for(unsigned int j = 0; j < m_Columns; j++){
      	temp = (*this)(row1, j);
      	(*this)(row1, j) = (*this)(row2, j);
      	(*this)(row2, j) = temp;
      }
      return *this;
    }
    /**
     * 列を入れ替えます。破壊的メソッドです。
     * 
     * @param column1 列インデックス1
     * @param column2 列インデックス2
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &exchangeColumns(const unsigned int &column1, const unsigned int &column2){
      if(column1 >= m_Columns || column2 >= m_Columns){}
      T temp;
      for(unsigned int i = 0; i < m_Rows; i++){
				temp = (*this)(i, column1);
      	(*this)(i, column1) = (*this)(i, column2);
      	(*this)(i, column2) = temp;
      }
      return *this;
    }
    
    /**
     * 代入演算子。
     * 
     * @return (Matrix<T>) 自分自身
     */
    Matrix &operator=(const Matrix &matrix){
      if((this == &matrix) || (m_Values == matrix.m_Values)){return *this;}
      if(m_Values && ((--(*ref)) <= 0)){delete [] m_Values; delete ref;}
      if(m_Values = matrix.m_Values){
        (*(ref = matrix.ref))++;
        m_Rows = matrix.rows();
        m_Columns = matrix.columns();
      }
      return *this;
    }
    
    /**
     * 正方行列かどうか調べます。
     * 
     * @return (bool) 正方行列である場合true、それ以外の場合false
     */
    bool isSquare() const{return m_Rows == m_Columns;}
    
    /**
     * 対称行列かどうか調べます。
     * 
     * @return (bool) 対称行列である場合true、それ以外の場合false
     */
    bool isSymmetric() const{
    	if(isSquare()){
    		for(unsigned int i = 0; i < m_Rows; i++){
    			for(unsigned int j = i + 1; j < m_Columns; j++){
    				if((*const_cast<Matrix *>(this))(i, j) != (*const_cast<Matrix *>(this))(j, i)){return false;}
    			}	
    		}
    		return true;
    	}else{return false;}
    } 
    
    /**
     * 行列の成分全てを指定倍します。破壊的メソッドです。
     * 
     * @param scalar 倍数
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator*=(const T &scalar){
      for(unsigned int i = 0; i < m_Rows; i++){
        for(unsigned int j = 0; j < m_Columns; j++){
          (*this)(i, j) *= scalar;
        }
      }
      return *this;
    }
    /**
     * 行列の成分全てを指定倍します。
     * 
     * @param scalar 倍数
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator*(const T &scalar) const{return (copy() *= scalar);}
    /**
     * 行列の成分全てを指定倍します。
     * 
     * @param scalar 倍数
     * @return (Matrix<T>) 結果
     */
    friend Matrix<T> operator*(const T &scalar, const Matrix<T> &matrix){return matrix * scalar;}
    /**
     * 行列の成分全てを除算します。破壊的メソッドです。
     * 
     * @param scalar 倍数
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator/=(const T &scalar){return (*this) *= (1 / scalar);}
    /**
     * 行列の成分全てを除算します。
     * 
     * @param scalar 倍数
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator/(const T &scalar) const{return (copy() /= scalar);}
    /**
     * 行列の成分全てを除算します。
     * 
     * @param scalar 倍数
     * @return (Matrix<T>) 結果
     */
    friend Matrix<T> operator/(const T &scalar, const Matrix<T> &matrix){return matrix / scalar;}
    
    /**
     * 行列を転置します。
     * 
     * @return (Matrix<T>) 転置行列
     */
    Matrix<T> transpose() const{
      Matrix<T> result(m_Columns, m_Rows);
      for(unsigned int i = 0; i < result.rows(); i++){
        for(unsigned int j = 0; j < result.columns(); j++){
          result(i, j) = (*const_cast<Matrix<T> *>(this))(j, i);
        }
      }
      return result;
    }
    
		/**
     * 行列を成分ごとに加算します。破壊的メソッドです。
     * 
     * @param matrix 加える行列
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator+=(const Matrix<T> &matrix){
    	if(m_Rows != matrix.rows() || m_Columns != matrix.columns()){}
		  for(unsigned int i = 0; i < m_Rows; i++){
    		for(unsigned int j = 0; j < m_Columns; j++){
      		(*this)(i, j) += (*const_cast<Matrix<T> *>(&matrix))(i, j);
    		}
  		}
  		return *this;
    }

    /**
     * 行列を成分ごとに加算します。
     * 
     * @param matrix 加える行列
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator+(const Matrix<T> &matrix) const{return (copy() += matrix);}
    
    /**
     * 行列を成分ごとに減算します。
     * 
     * @param matrix 引く行列
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator-=(const Matrix<T> &matrix){
    	if(m_Rows != matrix.rows() || m_Columns != matrix.columns()){}
		  for(int i = 0; i < m_Rows; i++){
    		for(int j = 0; j < m_Columns; j++){
      		(*this)(i, j) -= (*const_cast<Matrix<T> *>(&matrix))(i, j);
    		}
  		}
  		return *this;
    }
    
    /**
     * 行列を成分ごとに減算します。
     * 
     * @param matrix 引く行列
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator-(const Matrix<T> &matrix) const{return (copy() -= matrix);}
    
    /**
     * 行列を乗算します。
     * 
     * @param matrix かける行列
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator*(const Matrix<T> &matrix) const {
      if(m_Columns != matrix.rows()){}
      Matrix<T> result(m_Rows, matrix.columns());
      for(unsigned int i = 0; i < result.rows(); i++){
        for(unsigned int j = 0; j < result.columns(); j++){
          for(unsigned int k = 0; k < m_Columns; k++){
            result(i, j) += (*const_cast<Matrix<T> *>(this))(i, k) * (*const_cast<Matrix<T> *>(&matrix))(k, j);
          }
        }
      }
      return result;
    }
    /**
     * 行列を乗算します。破壊的メソッドです。
     * 
     * @param matrix かける行列
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator*=(const Matrix<T> &matrix) {
      return (*this = (*this * matrix));
    }
    
    /**
     * 単項演算子-。
     * 効果は matrix * -1と同じです。
     * 
     * @return (Matrix<T>) -matrix
     */
     Matrix<T> operator-() const{return (copy() *= -1);}
    
    /**
     * 補行列(余因子行列)を求めます。
     * 
     * @param row 行インデックス
     * @param column 列インデックス
     * @return (Matrix<T>) 補行列
     */
    Matrix<T> coMatrix(const int &row, const int &column) const {
      if(row < 0 && row >= m_Rows && column < 0 && column >= m_Rows){}
      Matrix<T> result(m_Rows - 1, m_Columns - 1);
      for(unsigned int i = 0; i < m_Rows - 1; i++){
        for(unsigned int j = 0; j < m_Columns - 1; j++){
          result(i, j) = (*const_cast<Matrix<T> *>(this))((i < row ? i : i + 1), (j < column ? j : j + 1));
        }  
      }
      return result;
    }
    
    /**
     * 行列式を計算します。
     * 
     * @return (T) 結果
     */
    T determinant() const {
      if(!isSquare()){}
      if(m_Rows == 1){
        return (*const_cast<Matrix<T> *>(this))(0, 0);
      }else{
        T sum(0);
        for(unsigned int i = 0; i < m_Rows; i++){
          if((*const_cast<Matrix<T> *>(this))(i, 0) != 0){
            sum += (*const_cast<Matrix<T> *>(this))(i, 0) * (coMatrix(i, 0).determinant()) * (i % 2 == 0 ? 1 : -1);
          }
        }
        return sum;
      }
    }
    
    /**
     * LU分解をします。
     * (0, 0)〜(n-1, n-1): L行列
     * (0, n)〜(n-1, 2n-1):U行列
     * 
     * @return (Matrix<T>) LU分解
     */
   	Matrix<T> divideLU() const {
   		if(!isSquare()){}
   		Matrix<T> LU(m_Rows, m_Columns * 2);
   		for(unsigned int i = 0; i < m_Rows; i++){
   			for(unsigned int j = 0; j < m_Columns; j++){
   				if(i >= j){
   					LU(i, j) = (*const_cast<Matrix *>(this))(i, j);
   					for(unsigned int k = 0; k < j; k++){
   						LU(i, j) -= (LU(i, k) * LU(k, j + m_Columns));
   					}
   				}else{
   					LU(i, j + m_Columns) = (*const_cast<Matrix *>(this))(i, j);
   					for(unsigned int k = 0; k < i; k++){
   						LU(i, j + m_Columns) -= (LU(i, k) * LU(k, j + m_Columns));
   					}
   					LU(i, j + m_Columns) /= LU(i, i);
   				}
   			}
   			LU(i, i + m_Columns) = 1;
   		}
   		return LU;
   	}
    
    /**
     * 逆行列を求めます。
     * 
     * @return (Matrix<T>) 逆行列
     */
    Matrix<T> inverse() const {
      if(!isSquare()){}
      
      //クラメール(遅い)
      /*
      Matrix<T> result(m_Rows, m_Columns);
      T det;
      if((det = determinant()) == 0){throw MatrixException("Operatiorn void!!");}
      for(int i = 0; i < m_Rows; i++){
        for(int j = 0; j < m_Columns; j++){
          result(i, j) = coMatrix(i, j).determinant() * ((i + j) % 2 == 0 ? 1 : -1);
        }
      }
      return result.transpose() / det;
      */
      
      //ガウス消去法
      Matrix<T> left = copy();
      Matrix<T> right = Matrix<T>::getI(m_Rows);
      if(left(0, 0) == T(0)){       //(0, 0)が存在するように並べ替え
        for(unsigned int i = 1; i <= m_Rows; i++){
          if(i == m_Rows){}
          if(left(i, 0) != T(0)){
            left.exchangeRows(0, i);
            right.exchangeRows(0, i);
            break;
          }
        }
      }
      for(unsigned int i = 0; i < m_Rows; i++){
        if(left(i, i) == T(0)){}
        if(left(i, i) != T(1)){
          for(unsigned int j = 0; j < m_Columns; j++){right(i, j) /= left(i, i);}
          for(unsigned int j = i+1; j < m_Columns; j++){left(i, j) /= left(i, i);}
          left(i, i) = T(1);
        }
        for(unsigned int k = 0; k < m_Rows; k++){
          if(k == i){continue;}
          if(left(k, i) != T(0)){
            for(unsigned int j = 0; j < m_Columns; j++){right(k, j) -= right(i, j) * left(k, i);}
            for(unsigned int j = i+1; j < m_Columns; j++){left(k, j) -= left(i, j) * left(k, i);}
            left(k, i) = T(0);
          }
        }
      }
      return right;
      
      //LU分解
      /*
      */
    }
    /**
     * 逆行列をかけます。破壊的メソッドです。
     * 
     * @param matrix 行列
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator/=(const Matrix<T> &matrix) {return (*this) *= matrix.inverse();}
    /**
     * 逆行列をかけます。
     * 
     * @param matrix 行列
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator/(const Matrix<T> &matrix) const {return (copy() /= matrix);}
    
    /**
     * ピボットを指定して、加算します。
     * 破壊的です。
     * 
     * @param row 行インデックス
     * @param column 列インデックス
     * @param matrix 足す行列
     */
    Matrix<T> pivotMerge(const int &row, const int &column, const Matrix &matrix){
      for(unsigned int i = 0; i < matrix.rows(); i++){
      	if(row + i < 0){continue;}
        else if(row + i >= m_Rows){break;}
        for(unsigned int j = 0; j < matrix.columns(); j++){
        	if(column + j < 0){continue;}
          else if(column + j >= m_Columns){break;}
          (*this)(row + i, column + j) += (*const_cast<Matrix *>(&matrix))(i, j);
        }
      }
      return *this;
    }
    
    /**
     * ピボットを指定して、加算します。
     * 
     * @param row 行インデックス
     * @param column 列インデックス
     * @param matrix 足す行列
     */
    Matrix<T> pivotAdd(const int &row, const int &column, const Matrix &matrix) const{
      return copy().pivotMerge(row, column, matrix);
    }
    
    /**
     * ハウスホルダー変換をしてヘッセンベルク行列を得ます。
     * 
     * @param transform 変換に用いた行列の積を格納するポインタ
     * @return (Matrix<T>) ヘッセンベルク行列
     */
    Matrix<T> hessenberg(Matrix *transform) const {
    	if(!isSquare()){}
    	
    	Matrix<T> result = copy();
    	for(unsigned int j = 0; j < m_Columns - 2; j++){
    		T t(0);
    		for(unsigned int i = j + 1; i < m_Rows; i++){
    			t += pow(result(i, j), 2);
    		}
        T s = sqrt(t);
        if(result(j + 1, j) < 0){s *= -1;}
        
    		Matrix<T> omega(m_Rows - (j+1), 1);
    		{
    			for(unsigned int i = 0; i < omega.rows(); i++){
    				omega(i, 0) = result(j+i+1, j);	
    			}
    			omega(0, 0) += s;
    		}
    		Matrix<T> P = Matrix<T>::getI(m_Rows).pivotMerge(j+1, j+1, -(omega * omega.transpose() / (t + result(j + 1, j) * s)));
    		
    		result = P * result * P;
        if(transform){(*transform) *= P;}
    	}
    	
    	//ゼロ処理
    	bool sym = isSymmetric();
    	for(unsigned int j = 0; j < m_Columns - 2; j++){
    		for(unsigned int i = j + 2; i < m_Rows; i++){
    			result(i, j) = T(0);
    			if(sym){result(j, i) = T(0);}
    		}	
    	}
    	
    	return result;	
    }
    
    /**
     * ハウスホルダー変換をしてヘッセンベルク行列を得ます。
     * 
     * @return (Matrix<T>) ヘッセンベルク行列
     */
    Matrix<T> hessenberg() const {return hessenberg(NULL);}
    
    /**
     * 行列を見やすい形で出力します。
     * 
     */
    void print() const{
    	if(m_Values){
    		printf("{");
	      for(unsigned int i = 0; i < m_Rows; i++){
  	      if(i > 0){printf(",");}
  	      printf("\n{");
    	    for(unsigned int j = 0; j < m_Columns; j++){
      	    if(j > 0){printf(",");}
      	    char c[16];
      	    sprintf(c, "%f", (*const_cast<Matrix<T> *>(this))(i, j));
      	    printf("%s", c);
        	}
	        printf("}");
  	    }
  	    printf("\n");
    	  printf("}\n");
    	}
    }
};
#endif /* __MATRIX_H__ */