#ifndef __MATRIX_H
#define __MATRIX_H

#include <string>
#include <exception>

/**
 * Matrixクラスの例外クラス。
 * 例として、演算が成立しない場合など。
 * 
 */
class MatrixException: public exception{
  private:
    string what_str;
  public:
    MatrixException(const string &what_arg): what_str(what_arg){}
    ~MatrixException() throw() {}
    /**
     * エラー内容を取得します。
     * 
     * @return (chsr *) エラー内容
     */
    const char *what() const throw(){
      return what_str.c_str(); 
    }
};

#include <string.h>
#include <ostream>
#include "complex.h"

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) throw(MatrixException) : m_Rows(rows), m_Columns(columns){
        m_Values = new T[rows * columns];
        ref = new int(0);
      if((!m_Values) || (!ref)){throw MatrixException("Lack of memory!!");}
      (*ref)++;
      for(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 throw(MatrixException){
      if(row >= m_Rows){throw MatrixException("Index Incorrect!!");}
      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 throw(MatrixException){
      if(column >= m_Columns){throw MatrixException("Index Incorrect!!");}
      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 int &row, const int &column) throw(MatrixException){
      if((row < 0 || row >= m_Rows) || (column < 0 || column >= m_Columns)){
        throw MatrixException("Index incorrect");
      }
      return *(m_Values + (row * m_Columns) + column);
    }
    /**
     * 指定した行列成分を返します。
     * 
     * @param row 行番号(開始番号は1〜)
     * @param column 列番号(開始番号は1〜)
     * @return (T) 成分
     */
    T &index(const int &row, const int &column) throw(MatrixException){
      return (*this)(row -1, column -1);
    }
    
    /**
     * 行列を複製(ディープコピー)します。
     * 
     * @return (Matrix<T>) コピー
     */
    Matrix<T> copy() const throw(MatrixException){
      Matrix<T> matrix(m_Rows, m_Columns);
      /*
      for(int i = 0; i < m_Rows; i++){
        for(int j = 0; j < m_Columns; j++){
          matrix(i, j) = (*const_cast<Matrix<T> *>(this))(i, j);
        }
      }
      */
            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) throw(MatrixException){
      if(row1 >= m_Rows || row2 >= m_Rows){throw MatrixException("Index incorrect");}
      T temp;
      for(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){throw MatrixException("Index incorrect");}
      T temp;
      for(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(int i = 0; i < m_Rows; i++){
                for(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(int i = 0; i < m_Rows; i++){
        for(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(int i = 0; i < result.rows(); i++){
        for(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()){throw MatrixException("Operatiorn void!!");}
          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){
        if(m_Rows != matrix.rows() || m_Columns != matrix.columns()){throw MatrixException("Operatiorn void!!");}
          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 throw(MatrixException){
      if(m_Columns != matrix.rows()){throw MatrixException("Operatiorn void!!");}
      Matrix<T> result(m_Rows, matrix.columns());
      for(int i = 0; i < result.rows(); i++){
        for(int j = 0; j < result.columns(); j++){
          for(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) throw(MatrixException){
      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 throw(MatrixException){
      if(row < 0 && row >= m_Rows && column < 0 && column >= m_Rows){
        throw MatrixException("Index incorrect");
      }
      Matrix<T> result(m_Rows -1, m_Columns -1);
      for(int i = 0; i < m_Rows -1; i++){
        for(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 throw(MatrixException){
      if(!isSquare()){throw MatrixException("Operatiorn void!!");}
      if(m_Rows == 1){
        return (*const_cast<Matrix<T> *>(this))(0, 0);
      }else{
        T sum(0);
        for(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 throw(MatrixException){
        if(!isSquare()){throw MatrixException("Operatiorn void!!");}
        Matrix<T> LU(m_Rows, m_Columns * 2);
        for(int i = 0; i < m_Rows; i++){
            for(int j = 0; j < m_Columns; j++){
                if(i >= j){
                    LU(i, j) = (*const_cast<Matrix *>(this))(i, j);
                    for(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(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 throw(MatrixException){
      if(!isSquare()){throw MatrixException("Operatiorn void!!");}
      
      //クラメール(遅い)
      /*
      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(int i = 1; i <= m_Rows; i++){
          if(i == m_Rows){throw MatrixException("Operatiorn void!! ; Invert matrix not exist");}
          if(left(i, 0) != T(0)){
            left.exchangeRows(0, i);
            right.exchangeRows(0, i);
            break;
          }
        }
      }
      for(int i = 0; i < m_Rows; i++){
        if(left(i, i) == T(0)){throw MatrixException("Operatiorn void!! ; Invert matrix not exist");}
        if(left(i, i) != T(1)){
          for(int j = 0; j < m_Columns; j++){right(i, j) /= left(i, i);}
          for(int j = i+1; j < m_Columns; j++){left(i, j) /= left(i, i);}
          left(i, i) = T(1);
        }
        for(int k = 0; k < m_Rows; k++){
          if(k == i){continue;}
          if(left(k, i) != T(0)){
            for(int j = 0; j < m_Columns; j++){right(k, j) -= right(i, j) * left(k, i);}
            for(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) throw(MatrixException){return (*this) *= matrix.inverse();}
    /**
     * 逆行列をかけます。
     * 
     * @param matrix 行列
     * @return (Matrix<T>) 結果
     */
    Matrix<T> operator/(const Matrix<T> &matrix) const throw(MatrixException){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 >= m_Rows){break;}
        for(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 throw(MatrixException){
        if(!isSquare()){throw MatrixException("Operation void!!");}
        
        Matrix<T> result = copy();
        for(int j = 0; j < m_Columns -2; j++){
            T t(0);
            for(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(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(int j = 0; j < m_Columns -2; j++){
            for(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 throw(MatrixException){return hessenberg(NULL);}
    
    /**
     * 2次小行列の固有値を求めます。
     * 
     * @param row 2次小行列の左上項の行インデックス
     * @param column 2次小行列の左上項の列インデックス
     * @param upper 結果(固有値1)
     * @param lower 結果(固有値2)
     */
    void eigen22(const int &row, const int &column, Complex<T> * &upper, Complex<T> * &lower) const throw(MatrixException){
      T a = (*const_cast<Matrix *>(this))(row, column),
        b = (*const_cast<Matrix *>(this))(row, column + 1),
        c = (*const_cast<Matrix *>(this))(row + 1, column),
        d = (*const_cast<Matrix *>(this))(row + 1, column + 1);
      T root = pow((a -d), 2) + b * c * 4;
      if(root >= T(0)){
        root = sqrt(root);
        upper = new Complex<T>((a + d + root) / 2);
        lower = new Complex<T>((a + d -root) / 2);
      }else{
        root = sqrt(root * -1);
        upper = new Complex<T>((a + d) / 2, root / 2);
        lower = new Complex<T>((a + d) / 2, root / 2 * -1);
      }
    }
    
    /**
     * 固有値、固有ベクトルを求めます。
     * 返却値はMatrix<Complex<T> >型で、
     * (0,0)〜(n-1,n-1)要素が固有ベクトルの行列
     * (0,n)〜(n-1,n)要素が対応する固有値の列ベクトル
     * になっています。
     * (固有ベクトル1,固有ベクトル2,…,固有ベクトルn,固有値)のn×(n+1)行列
     * 
     * @param threshold 収束判定に用いる閾値
     * @return (Matrix<T>) 固有値、固有ベクトル
     */
    Matrix<Complex<T> > eigen(const T &threshold) const throw(MatrixException){
        
        if(!isSquare()){throw MatrixException("Operatiorn void!!");}
        
      //パワー法(べき乗法)
      /*Matrix<T> result(m_Rows, m_Rows + 1);
      Matrix<T> source = copy();
      for(int i = 0; i < m_Columns; i++){result(0, i) = T(1);}
      for(int i = 0; i < m_Columns; i++){
        while(true){
          Matrix<T> approxVec = source * result.columnVector(i);
          T approxVal(0);
          for(int j = 0; j < approxVec.rows(); j++){approxVal += pow(approxVec(j, 0), 2);}
          approxVal = sqrt(approxVal);
          for(int j = 0; j < approxVec.rows(); j++){result(j, i) = approxVec(j, 0) / approxVal;}
          T before = result(i, m_Rows);
          if(abs(before - (result(i, m_Rows) = approxVal)) < threshold){break;}
        }
        for(int j = 0; (i < m_Rows - 1) && (j < m_Rows); j++){
          for(int k = 0; k < m_Rows; k++){
            source(j, k) -= result(i, m_Rows) * result(j, i) * result(k, i);
          }
        }
      }
      return result;*/
      
      //ダブルQR法
      /* <手順>
       * ハウスホルダー法を適用して、上ヘッセンベルク行列に置換後、
       * ダブルQR法を適用。
       * 結果、固有値が得られるので、固有ベクトルを計算。
       */
      
      //固有値の計算
      Matrix<Complex<T> > lambda(m_Rows, 1);   //固有値
      
      T mu_sum = T(0), mu_multi = T(0);
      Complex<T> p1, p2;
      int m = m_Rows;
      bool first = true;
      
      Matrix transform = getI(m_Rows);
      Matrix A = hessenberg(&transform);
      Matrix A_ = A;

      while(true){
              
        //m = 1 or m = 2
        if(m == 1){
          lambda(0, 0) = A(0, 0);
          break;
        }else if(m == 2){
          Complex<T> *lambda1, *lambda2;
          A.eigen22(0, 0, lambda1, lambda2);
          lambda(0, 0) = *lambda1;
          lambda(1, 0) = *lambda2;
          break;
        }
        
        //μ、μ*の更新(4.143)
        {
          Complex<T> *p1_new, *p2_new;
          A.eigen22(m-2, m-2, p1_new, p2_new);
          if(first ? (first = false) : true){
            if((*p1_new -p1).abs() > p1_new->abs() / 2){
              if((*p2_new -p2).abs() > p2_new->abs() / 2){
                mu_sum = (p1 + p2).real();
                mu_multi = (p1 * p2).real();
              }else{
                mu_sum = p2_new->real() * 2;
                mu_multi = pow(p2_new->real(), 2);
              }
            }else{
              if((*p2_new -p2).abs() > p2_new->abs() / 2){
                mu_sum = p1_new->real() * 2;
                mu_multi = pow(p1_new->real(), 2);
              }else{
                mu_sum = (*p1_new + *p2_new).real();
                mu_multi = (*p1_new * *p2_new).real();
              }
            }
          }
          p1 = *p1_new, p2 = *p2_new;
        }
      
        //ハウスホルダー変換を繰り返す
        T b1, b2, b3, r;
        for(int i = 0; i < m -1; i++){
            if(i == 0){
                b1 = pow(A(0, 0), 2) -mu_sum * A(0, 0) + mu_multi + A(0, 1) * A(1, 0);
            b2 = A(1, 0) * (A(0, 0) + A(1, 1) -mu_sum);
                b3 = A(2, 1) * A(1, 0);
            }else{
                b1 = A(i, i -1);
            b2 = A(i + 1, i -1);
            b3 = (i == m -2 ? T(0) : A(i + 2, i -1));
            }
        
            r = sqrt(pow(b1, 2) + pow(b2, 2) + pow(b3, 2));
        
          Matrix omega(3, 1);
          {
              omega(0, 0) = b1 + r * (b1 >= T(0) ? 1 : -1);
                omega(1, 0) = b2;
                if(b3 != T(0)){omega(2, 0) = b3;}
          }
          Matrix P = getI(m_Rows).pivotMerge(i, i, omega * omega.transpose() * -2 / (omega.transpose() * omega)(0, 0));
       
          A = P * A * P;
        }
        
        //収束判定
        T epsilon = threshold
          + threshold * abs(abs(A(m-2, m-2)) < abs(A(m-1, m-1)) ? A(m-2, m-2) : A(m-1, m-1));
        if(abs(A(m-1, m-2)) < epsilon){
          lambda(--m, 0) = A(m, m);
        }else if(abs(A(m-2, m-3)) < epsilon){
          Complex<T> *lambda1, *lambda2;
          A.eigen22(m-2, m-2, lambda1, lambda2);
          lambda(--m, 0) = *lambda1;
          lambda(--m, 0) = *lambda2;
        }
      }
      
      //固有ベクトルの計算
      Matrix<Complex<T> > x(m_Rows, m_Rows);  //固有ベクトル
      A = A_;
      
      for(int j = 0; j < m_Rows; j++){
        int n = m_Rows;
        for(int i = 0; i < j; i++){if(lambda(j, 0) == lambda(i, 0)){--n;}}
        x(--n, j) = 1;
        while(n--> 0){
          x(n, j) = x(n+1, j) * (lambda(j, 0) -A(n+1, n+1));
          for(int i = n+2; i < m_Rows; i++){
            x(n, j) -= x(i, j) * A(n+1, i);
          }
          x(n, j) /= A(n+1, n);
        }
      }
      
      /*Matrix<Complex<T> > lambda2(m_Rows, m_Rows);
      for(int i = 0; i < m_Rows; i++){
        lambda2(i, i) = lambda(i, 0);
      }
      
      cout << "A:" << A << endl;
      //cout << "x * x^-1" << x * x.inverse() << endl;
      cout << "x * lambda * x^-1:" << x * lambda2 * x.inverse() << endl;*/
            
      //結果の格納
      Matrix<Complex<T> > result(m_Rows, m_Rows + 1);
      for(int j = 0; j < x.columns(); j++){
        result(j, m_Rows) = lambda(j, 0);
        for(int i = 0; i < x.rows(); i++){
          for(int k = 0; k < transform.columns(); k++){
            result(i, j) += transform(i, k) * x(k, j);
          }
        }
        
        //正規化
        Complex<T> _norm;
        for(int i = 0; i < m_Rows; i++){
            _norm += result(i, j).abs2();
        }
        T norm = sqrt(_norm.real());        
        for(int i = 0; i < m_Rows; i++){
            result(i, j) /= norm;   
        }  
      }
      
      return result;
    }
    
    /**
     * 固有値、固有ベクトルを求めます。
     * 返却値はMatrix<T>型で、
     * (0,0)〜(n-1,n-1)要素が固有ベクトルの行列
     * (0,n)〜(n-1,n)要素が対応する固有値の列ベクトル
     * になっています。
     * (固有ベクトル1,固有ベクトル2,…,固有ベクトルn,固有値)のn×(n+1)行列
     * 
     * @return (Matrix<T>) 固有値、固有ベクトル
     */
    Matrix<Complex<T> > eigen() const throw(MatrixException){return eigen(10E-10);}
    
    /**
     * 行列を見やすい形で出力します。
     * 
     */
    friend ostream &operator<<(ostream &out, const Matrix<T> &matrix){
        if(matrix.m_Values){
            out << "{";
          for(int i = 0; i < matrix.rows(); i++){
          out << (i == 0 ? "" : ",") << endl << "{";
            for(int j = 0; j < matrix.columns(); j++){
            out << (j == 0 ? "" : ",") << (*const_cast<Matrix<T> *>(&matrix))(i, j);
            }
            out << "}";
        }
          out << endl << "}";
        }
      return out;  
    }
};

#endif /* __MATRIX_H */