/**
 * THIS LIBRARY HAS SEVERAL BUGS AND DEPRECATED !!
 * PLEASE CHECK https://fenrir.naruoka.org/archives/000595.html
 * 
 */


/*extern "C"{
#include <DSPF_dp_mat_mul.h>
#include <DSPF_sp_mat_mul.h>
#include <DSPF_dp_mat_trans.h>
#include <DSPF_sp_mat_trans.h>
}*/

#include <cstdio>
#include <cstring>
#include <cassert>
#include <new>

//#include "matrix.hpp"

namespace fenrir{

using namespace std;

/**
 * 2次元配列の抽象クラス
 * 
 */
template<class T>
class Array2D{
  protected:
    unsigned int m_rows;
    unsigned int m_columns;
  
  public:
	/**
     * シャローコピーをします。
     * 
     * @return (Array2D *) 自分自身
     */
    virtual Array2D<T> *shallow_copy() const = 0;

    /**
     * Array2Dクラスのコンストラクタ。
     * 指定の行サイズ、指定の列サイズで仮想的な2次元配列を生成します。
     * 
     * @param rows 行数
     * @param columns 列数
     */
    Array2D(const unsigned int &rows, const unsigned int &columns)
        : m_rows(rows), m_columns(columns){set_new_handler(0);}
    
    /**
     * Array2Dクラスのデストラクタ。

     */
    virtual ~Array2D(){}
    
    /**
     * 2次元配列を複製(ディープコピー)します。
     * 
     * @return (Array2D<T>) コピー
     */
    virtual Array2D<T> *copy() const = 0;

    /**
	 * Tの配列に落とし込みます。
	 *
	 * @return (T*)
	 */
    virtual T* serialize() const = 0;
    
    /**
     * 要素のゼロクリアをします。
     * 
     */
    virtual void clear() = 0;
    
    /**
     * 行数を返します。
     * 
     * @return (int) 行数
     */
    virtual unsigned int rows() const{return m_rows;}
    /**
     * 列数を返します。
     * 
     * @return (int) 列数
     */
    virtual unsigned int columns() const{return m_columns;}
    
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    virtual T &operator()(const unsigned int &row, const unsigned int &column) = 0;
};

/**
 * 中身が詰まった2次元配列
 * 
 */
template <class T>
class Array2D_Dense : public Array2D<T> {
  protected:
    T *m_Values;
    int *ref; //参照カウンタ
  
  public:
    /**
     * シャローコピーをします。
     * 参照カウンタのインクリメントも同時に行います。
     * 
     * @return (Array2D_Dense *)自分自身
     */
    Array2D<T> *shallow_copy() const{return new Array2D_Dense<T>(*this);}
    
    /**
     * 2次元配列を複製(ディープコピー)します。
     * 
     * @return (Array2D<T>) コピー
     */
    Array2D<T> *copy() const {
      Array2D_Dense<T> *array = new Array2D_Dense<T>(Array2D<T>::rows(), Array2D<T>::columns());
      memcpy(array->m_Values, m_Values, sizeof(T) * Array2D<T>::rows() * Array2D<T>::columns());
      return array;
    }

    /**
     * Tの配列に落とし込みます。
     *
     * @return (T*)
     */
    T* serialize() const {return m_Values;}
    
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    T &operator()(const unsigned int &row, const unsigned int &column){
      assert((row < Array2D<T>::rows()) && (column < Array2D<T>::columns()));
      return *(m_Values + (row * Array2D<T>::columns()) + column);
    }
    
    /**
     * 要素のゼロクリアをします。
     * 
     */
    void clear(){
      for(unsigned int i = 0; i < Array2D<T>::rows() * Array2D<T>::columns(); i++){*(m_Values + i) = T(0);}
    }
    
    /**
     * Array2D_Denseクラスのコンストラクタ。
     * 指定の行サイズ、指定の列サイズで2次元配列を生成します。
     * また成分はすべてT(0)で初期化されます。
     * 
     * @param rows 行数
     * @param columns 列数
     */
    Array2D_Dense<T>(const unsigned int &rows, const unsigned int &columns) 
        : Array2D<T>(rows, columns), m_Values(new T[rows * columns]), ref(new int(0)) {
      assert(m_Values && ref);
      (*ref)++;
      clear();
    }
    
    /**
     * コピーコンストラクタ
     * 
     */
    Array2D_Dense(const Array2D_Dense &array) : Array2D<T>(array.m_rows, array.m_columns){
      if(m_Values = array.m_Values){(*(ref = array.ref))++;}
    }
    /**
     * デストラクタ。
     */
    ~Array2D_Dense(){
      if(ref && ((--(*ref)) <= 0)){
        delete [] m_Values;
        delete ref;
      }
    }
    
    /**
     * 代入演算子。
     * 
     * @return (Array2D_Dense<T>) 自分自身
     */
    Array2D_Dense &operator=(const Array2D_Dense &array){
      if(this != &array){
        if(ref && ((--(*ref)) <= 0)){delete [] m_Values;}
        if(m_Values = array.m_Values){
          Array2D<T>::m_rows = array.m_rows;
          Array2D<T>::m_columns = array.m_columns;
          (*(ref = array.ref))++;
        }        
      }
      return *this;
    }  
};

/**
 * 別の2次元配列に委譲を行う2次元配列クラス
 * 
 */
template<class T>
class Array2D_Delegate : public Array2D<T>{
  private:
    Array2D<T> *m_target;
    
  protected:
    /**
     * 委譲先を返します。
     * 
     */
    Array2D<T> &getTarget() const{return *m_target;}
  
  public:
    
    /**
     * Array2D_Partial(部分2次元配列)クラスのコンストラクタ。
     * 
     * @param rows 行数
     * @param columns 列数
     * @param array 元の配列
     */
    Array2D_Delegate(const unsigned int &rows,
                     const unsigned int &columns,
                     const Array2D<T> &array)
                        : Array2D<T>(rows, columns), m_target(array.shallow_copy()){}
    
    /**
     * コピーコンストラクタ
     * 
     */
    Array2D_Delegate(const Array2D_Delegate &array)
      : Array2D<T>(array.rows(), array.columns()), m_target(array.m_target->shallow_copy()){}
        
    /**
     * デストラクタ。
     */
    ~Array2D_Delegate(){delete m_target;}
        
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    T &operator()(const unsigned int &row, const unsigned int &column){
      return m_target->operator()(row, column);
    }
    
    /**
     * 部分2次元配列を複製(ディープコピー)します。
     * 複製後は単純な2次元配列になります。
     * 
     * @return (Array2D<T>) コピー
     */
    Array2D<T> *copy() const{
      Array2D_Dense<T> *array = new Array2D_Dense<T>(Array2D<T>::rows(), Array2D<T>::columns());
      for(unsigned int i = 0; i < Array2D<T>::rows(); i++){
        for(unsigned int j = 0; j < Array2D<T>::columns(); j++){
          array->operator()(i, j) = (const_cast<Array2D_Delegate *>(this))->operator()(i, j);
        }
      }
      return array;
    }

    /**
     * Tの配列に落とし込みます。
     *
     * @return (T*)
     */
    T* serialize() const {return m_target->serialize();}
    
    /**
     * 要素のゼロクリアを行います。
     * 
     */
    void clear(){
      for(unsigned int i = 0; i < Array2D<T>::rows(); i++){
        for(unsigned int j = 0; j < Array2D<T>::columns(); j++){
          this->operator()(i, j) = T(0);
        }
      } 
    }
};

/**
 * 転置2次元配列をあらわすクラス
 * 
 */
template<class T>
class Array2D_Transpose : public Array2D_Delegate<T>{
  
  public:
    /**
     * シャローコピーをします。
     * 
     * @return (Array2D *)自分自身
     */
    Array2D<T> *shallow_copy() const{return new Array2D_Transpose(*this);}
    
    /**
     * Array2D_Transpose(部分2次元配列)クラスのコンストラクタ。
     * 
     * @param array 元の配列
     */
    Array2D_Transpose(const Array2D<T> &array)
                       : Array2D_Delegate<T>(array.columns(), array.rows(), array){}
    
    /**
     * コピーコンストラクタ
     * 
     */
    Array2D_Transpose(const Array2D_Transpose &array)
      : Array2D_Delegate<T>(array.rows(), array.columns(), array.getTarget()){}
        
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    T &operator()(const unsigned int &row, const unsigned int &column){
      return Array2D_Delegate<T>::operator()(column, row);
    }

    /**
     * Tの配列に落とし込みます。
     *
     * @return (T*)
     */
    T* serialize() const {
      T *serialized = new T[sizeof(T) * Array2D_Delegate<T>::rows() * Array2D_Delegate<T>::columns()];
	  int index(0);
      for(int i = 0; i < Array2D_Delegate<T>::rows(); i++){
        for(int j = 0; j < Array2D_Delegate<T>::columns(); j++){
          serialized[index++] = const_cast<Array2D_Transpose *>(this)->operator()(i, j);
        }
      }
      return serialized;
    }
};


/**
 * 部分2次元配列をあらわすクラス
 * 
 */
template<class T>
class Array2D_Partial : public Array2D_Delegate<T>{
  
  private:
    unsigned int m_RowOffset;
    unsigned int m_ColumnOffset;
    
  protected:
    /**
     * オフセット行を返します。
     * 
     * @return (int) 行数
     */
    unsigned int row_offset() const{return m_RowOffset;}
    /**
     * オフセット列を返します。
     * 
     * @return (int) 列数
     */
    unsigned int column_offset() const{return m_ColumnOffset;}
  
  public:
    /**
     * シャローコピーをします。
     * 
     * @return (Array2D *)自分自身
     */
    Array2D<T> *shallow_copy() const{return new Array2D_Partial(*this);}
    
    /**
     * Array2D_Partial(部分2次元配列)クラスのコンストラクタ。
     * 
     * @param array 元の配列
     * @param offsetRow 部分2次元配列で(0,0)となる元の2次元配列のピボットの行インデックス
     * @param offsetColumn 同じく列インデックス
     * @param rows 行数
     * @param columns 列数
     */
    Array2D_Partial(const unsigned int &rows,
                    const unsigned int &columns,
                    const Array2D<T> &array,
                    const unsigned int &rowOffset,
                    const unsigned int &columnOffset)
                       : Array2D_Delegate<T>(rows, columns, array),
                         m_RowOffset(rowOffset), m_ColumnOffset(columnOffset){}
    
    /**
     * コピーコンストラクタ
     * 
     */
    Array2D_Partial(const Array2D_Partial &array)
      : Array2D_Delegate<T>(array.rows(), array.columns(), array.getTarget()),
        m_RowOffset(array.row_offset()), m_ColumnOffset(array.column_offset()){}
        
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    T &operator()(const unsigned int &row, const unsigned int &column){
      return Array2D_Delegate<T>::operator()(row + row_offset(), column + column_offset());
    }

    /**
     * Tの配列に落とし込みます。
     *
     * @return (T*)
     */
    T* serialize() const {
      T *serialized = new T[sizeof(T) * Array2D_Delegate<T>::rows() * Array2D_Delegate<T>::columns()];
	  int index(0);
      for(int i = 0; i < Array2D_Delegate<T>::rows(); i++){
        for(int j = 0; j < Array2D_Delegate<T>::columns(); j++){
          serialized[index++] = const_cast<Array2D_Partial *>(this)->operator()(i, j);
        }
      }
      return serialized;
    }
};

template<class T>
class Matrix;

/**
 * 転置行列をあらわすクラス
 * 
 */
template <class T>
class TransposeMatrix : public Matrix<T>{
  
  public:
    /**
     * TransposeMatrix(転置行列)クラスのコンストラクタ。
     * 
     * @param matrix 元の行列
     */
    TransposeMatrix(const Matrix<T> &matrix)
                      : Matrix<T>(new Array2D_Transpose<T>(*(matrix.strage()))){}

    /**
     * デストラクタ。
     */
    ~TransposeMatrix(){}
};

/**
 * 部分行列をあらわすクラス
 * 
 */
template <class T>
class PartialMatrix : public Matrix<T>{
  
  public:
    /**
     * PartialMatrix(部分行列)クラスのコンストラクタ。
     * 
     * @param matrix 元の行列
     * @param offsetRow 部分行列で(0,0)となる元の行列のピボットの行インデックス
     * @param offsetColumn 同じく列インデックス
     * @param rows 行数
     * @param columns 列数
     */
    PartialMatrix(const Matrix<T> &matrix,
                  const unsigned int &rows,
                  const unsigned int &columns,
                  const unsigned int &rowOffset,
                  const unsigned int &columnOffset)
                    : Matrix<T>(new Array2D_Partial<T>(rows, columns, *(matrix.strage()), rowOffset, columnOffset)){}

    /**
     * デストラクタ。
     */
    ~PartialMatrix(){}
};

/**
 * 行列のクラス
 * 
 */
template <class T = double>
class Matrix{
  protected:
    Array2D<T> *m_Strage;
    
    /**
     * Matrixクラスのコンストラクタ。
     * ストレージを指定して新しく行列を作ります。
     * 
     * @param array ストレージ
     */
    Matrix(const Array2D<T> *strage) : m_Strage(const_cast<Array2D<T> *>(strage)){}
  public:
    /**
     * 内部的な保存形式を返します。
     * 
     * @return (const Array2D<T> *) ストレージ
     */
    const Array2D<T> *strage() const{return m_Strage;}

    /**
     * Matrixクラスのコンストラクタ。
     * 
     */
    Matrix() : m_Strage(NULL){}
    
    /**
     * Matrixクラスのコンストラクタ。
     * 指定の行数、指定の列数で行列を生成します。
     * また成分はすべてT(0)で初期化されます。
     * 
     *     
     * @param rows 行数
     * @param columns 列数
     */
    Matrix(const unsigned int &rows, const unsigned int &columns)
      : m_Strage(new Array2D_Dense<T>(rows, columns)){}
    
    /**
     * コピーコンストラクタ。
     * シャローコピーを生成します。
     * 
     * @param matrix コピー元
     */
    Matrix(const Matrix &matrix) : m_Strage(matrix.m_Strage->shallow_copy()){}
    /**
     * デストラクタ。
     */
    virtual ~Matrix(){delete m_Strage;}
    
    /**
     * 代入演算子。
     * 
     * @return (Matrix<T>) 自分自身
     */
    Matrix &operator=(const Matrix &matrix){
      if(this != &matrix){
        delete m_Strage;
        if(matrix.m_Strage){
          m_Strage = matrix.m_Strage->shallow_copy();
        }
      }
      return *this;
    }
    
    /**
     * 行列を複製(ディープコピー)します。
     * 
     * @return (Matrix<T>) コピー
     */
    Matrix<T> copy() const{
      return Matrix<T>(m_Strage->copy());
    }
    
    /**
     * 行数を返します。
     * 
     * @return (int) 行数
     */
    unsigned int rows() const{return m_Strage->rows();}
    /**
     * 列数を返します。
     * 
     * @return (int) 列数
     */
    unsigned int columns() const{return m_Strage->columns();}
    
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行インデックス(開始番号は0〜)
     * @param column 列インデックス(開始番号は0〜)
     * @return (T) 成分
     */
    T &operator()(const unsigned int &row, const unsigned int &column){
      return m_Strage->operator()(row, column);
    }
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行番号(開始番号は1〜)
     * @param column 列番号(開始番号は1〜)
     * @return (T) 成分
     */
    T &index(const int &row, const int &column){
      return (*this)(row - 1, column - 1);
    }
    
    /**
     * 要素をゼロクリアします
     * 
     */
    Matrix<T> clear(){
      m_Strage->clear();
      return *this;
    }
    
    /**
     * 指定の単位行列を生成します。
     * 
     * @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;
    }
    
    /**
     * 行列を転置します。
     * 転置された行列はもとの行列とリンクしています。
     * もとの行列との切り離しを行うにはtranspose().copy()としてください。
     * 
     * @return (Matrix<T>) 転置行列
     */
    Matrix<T> transpose() const{return TransposeMatrix<T>(*this);}
    
    /**
     * 指定した行の行ベクトルを返します。
     * 
     * @param row 行インデックス
     * @return (Matrix<T>) 行ベクトル
     */
    Matrix<T> rowVector(const unsigned int &row) const{
      assert(row < rows());
      return PartialMatrix<T>(*this, 1, columns(), row, 0);
    }
    /**
     * 指定した列の列ベクトルを返します。
     * 
     * @param column 列インデックス
     * @return (Matrix<T>) 列ベクトル
     */
    Matrix<T> columnVector(const unsigned int &column) const{
      assert(column < columns());
      return PartialMatrix<T>(*this, rows(), 1, 0, column);
    }
    
    /**
     * 行を入れ替えます。破壊的メソッドです。
     * 
     * @param row1 行インデックス1
     * @param row2 行インデックス2
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &exchangeRows(const unsigned int &row1, const unsigned int &row2){
      assert(row1 < rows() && row2 < rows());
      T temp;
      for(unsigned int j = 0; j < 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){
      assert(column1 < columns() && column2 < columns());
      T temp;
      for(unsigned int i = 0; i < rows(); i++){
				temp = (*this)(i, column1);
      	(*this)(i, column1) = (*this)(i, column2);
      	(*this)(i, column2) = temp;
      }
      return *this;
    }
    
    /**
     * 正方行列かどうか調べます。
     * 
     * @return (bool) 正方行列である場合true、それ以外の場合false
     */
    bool isSquare() const{return rows() == columns();}
    
    /**
     * 対称行列かどうか調べます。
     * 
     * @return (bool) 対称行�である場合true、それ以外の場合false
     */
    bool isSymmetric() const{
      if(isSquare()){
        for(unsigned int i = 0; i < rows(); i++){
          for(unsigned int j = i + 1; j < 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 < rows(); i++){
        for(unsigned int j = 0; j < 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;}
    
    /**
     * 行列を成分ごとに加算します。破壊的メソッドです。
     * 
     * @param matrix 加える行列
     * @return (Matrix<T>) 自分自身
     */
    Matrix<T> &operator+=(const Matrix<T> &matrix){
    	assert(rows() == matrix.rows() && columns() == matrix.columns());
		for(unsigned int i = 0; i < rows(); i++){
    	  for(unsigned int j = 0; j < 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){
    	assert(rows() == matrix.rows() && columns() == matrix.columns());
		for(unsigned int i = 0; i < rows(); i++){
    	  for(unsigned int j = 0; j < 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{
      assert(columns() == matrix.rows());
      Matrix<T> result(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 < 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 unsigned int &row, const unsigned int &column) const{
      assert(row < rows() && column < columns());
      Matrix<T> result(rows() - 1, columns() - 1);
      for(int i = 0; i < rows() - 1; i++){
        for(int j = 0; j < 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(rows() == 1){
        return (*const_cast<Matrix<T> *>(this))(0, 0);
      }else{
        T sum(0);
        for(int i = 0; i < 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> decompressLU() const{
   		Matrix<T> LU(rows(), columns() * 2);
#define L(i, j) LU(i, j)
#define U(i, j) LU(i, j + columns())
   		for(int i = 0; i < rows(); i++){
   			for(int j = 0; j < columns(); j++){
   				if(i >= j){
   					L(i, j) = (*const_cast<Matrix *>(this))(i, j);
   					for(int k = 0; k < j; k++){
   						L(i, j) -= (L(i, k) * U(k, j));
   					}
   				}else{
   					U(i, j) = (*const_cast<Matrix *>(this))(i, j);
   					for(int k = 0; k < i; k++){
   						U(i, j) -= (L(i, k) * U(k, j));
   					}
   					U(i, j) /= L(i, i);
   				}
   			}
   			U(i, i) = 1;
   		}
#undef L
#undef U
   		return LU;
   	}
     
    /**
     * UD分解をします。
     * (0, 0)〜(n-1,n-1):  U行列
     * (0, n)〜(n-1,2n-1): D行列
     * 
     * @return (Matrix<T>) UD分解
     */
    Matrix<T> decompressUD() const{
      Matrix<T> P(copy());
      Matrix<T> UD(rows(), columns() * 2);
#define U(i, j) UD(i, j)
#define D(i, j) UD(i, j + columns())
      for(int i = rows() - 1; i >= 0; i--){
        D(i, i) = P(i, i);
        U(i, i) = T(1);
        for(int j = 0; j < i; j++){
          U(j, i) = P(j, i) / D(i, i);
          for(int k = 0; k <= j; k++){
            P(k, j) -= U(k, i) * D(i, i) * U(j, i);
          }
        }
      }
#undef U
#undef D
      return UD;
    }
    
    /**
     * 逆行列を求めます。
     * 
     * @return (Matrix<T>) 逆行列
     */
    Matrix<T> inverse() const{
      
      assert(isSquare());
      
      //ガウス消去法
      Matrix<T> left(copy());
      Matrix<T> right(Matrix<T>::getI(rows()));
      if(left(0, 0) == T(0)){       //(0, 0)が存在するように並べ替え
        for(unsigned int i = 1; i <= rows(); i++){
          assert(i != rows());
          if(left(i, 0) != T(0)){
            left.exchangeRows(0, i);
            right.exchangeRows(0, i);
            break;
          }
        }
      }
      for(unsigned int i = 0; i < rows(); i++){
        assert(left(i, i) != T(0));
        if(left(i, i) != T(1)){
          for(unsigned int j = 0; j < columns(); j++){right(i, j) /= left(i, i);}
          for(unsigned int j = i+1; j < columns(); j++){left(i, j) /= left(i, i);}
          left(i, i) = T(1);
        }
        for(unsigned int k = 0; k < rows(); k++){
          if(k == i){continue;}
          if(left(k, i) != T(0)){
            for(unsigned int j = 0; j < columns(); j++){right(k, j) -= right(i, j) * left(k, i);}
            for(unsigned int j = i+1; j < columns(); j++){left(k, j) -= left(i, j) * left(k, i);}
            left(k, i) = T(0);
          }
        }
      }
      //cout << "L:" << left << endl;
      //cout << "R:" << right << endl;
      
      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(int i = 0; i < matrix.rows(); i++){
      	if(row + i < 0){continue;}
        else if(row + i >= rows()){break;}
        for(int j = 0; j < matrix.columns(); j++){
        	if(column + j < 0){continue;}
          else if(column + j >= 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);
    }
    
    /**
     * 行列を見やすい形で出力します。
     * 
     */
    void inspect(char *buffer, int buffer_size) const{
      if(m_Strage){
        int printed;
        printed = snprintf(buffer, buffer_size, "{");
		buffer += printed;
		buffer_size -= printed;
	    for(unsigned int i = 0; i < rows(); i++){
  	      printed = snprintf(buffer, buffer_size, (i == 0 ? "\n{" : ",\n{"));
          buffer += printed;
          buffer_size -= printed;
    	  for(unsigned int j = 0; j < columns(); j++){
            /*printed = printf((j == 0 ? "%f" : ",%f"), 
      	                       const_cast<Matrix<T> *>(this)->operator()(i, j));*/
      	    printed = snprintf(buffer, buffer_size, (j == 0 ? "%f" : ",%f"), 
      	                       const_cast<Matrix<T> *>(this)->operator()(i, j));
            buffer += printed;
            buffer_size -= printed;
          }
	      printed = snprintf(buffer, buffer_size, "}");
          buffer += printed;
          buffer_size -= printed;
  	    }
        snprintf(buffer, buffer_size, "\n}");
      }
    }
};

#ifdef DSPF_DP_MAT_MUL_H_
/**
 * 行列を乗算します。
 * doubleの場合の高速化
 * 
 * @param matrix かける行列
 * @return (Matrix<double>) 結果
 */
template<>
Matrix<double> Matrix<double>::operator*(const Matrix<double> &matrix) const{
  assert(columns() == matrix.rows());
  Matrix<double> result(rows(), matrix.columns());
  DSPF_dp_mat_mul(strage()->serialize(), rows(), columns(), 
                  matrix.strage()->serialize(), matrix.columns(), 
                  result.strage()->serialize());  
  return result;
}
#endif

#ifdef DSPF_DP_MAT_TRANS_ASM_H_
/**
 * Tの配列に落とし込みます。
 * doubleの場合の高速化
 *
 * @return (T*)
 */
template<>
double* Array2D_Transpose<double>::serialize() const {
  double *serialized(new double[sizeof(double) * Array2D_Delegate<double>::rows() * Array2D_Delegate<double>::columns()]);
  DSPF_dp_mat_trans(Array2D_Delegate<double>::serialize(), Array2D_Delegate<double>::rows(), Array2D_Delegate<double>::columns(), serialized);
  return serialized;
}
#endif

#ifdef DSPF_SP_MAT_MUL_ASM_H_
/**
 * 行列を乗算します。
 * floatの場合の高速化
 * 
 * @param matrix かける行列
 * @return (Matrix<float>) 結果
 */
template<>
Matrix<float> Matrix<float>::operator*(const Matrix<float> &matrix) const{
  assert(columns() == matrix.rows());
  Matrix<float> result(rows(), matrix.columns());
  DSPF_sp_mat_mul(strage()->serialize(), rows(), columns(), 
                  matrix.strage()->serialize(), matrix.columns(), 
                  result.strage()->serialize());
  return result;
}
#endif

#ifdef DSPF_SP_MAT_TRANS_H_
/**
 * Tの配列に落とし込みます。
 * floatの場合の高速化
 *
 * @return (T*)
 */
template<>
float* Array2D_Transpose<float>::serialize() const {
  float *serialized(new float[sizeof(float) * Array2D_Delegate<float>::rows() * Array2D_Delegate<float>::columns()]);
  DSPF_sp_mat_trans(Array2D_Delegate<float>::serialize(), Array2D_Delegate<float>::rows(), Array2D_Delegate<float>::columns(), serialized);
  return serialized;
}
#endif

}