C++14による正方行列クラスの作成（その５）

前回は除算以外の算術二項演算子と論理演算子を定義した。
今回で、いよいよ行列式と余因子行列を求める。
この実装は余因子展開を使うため、行列の一時オブジェクトをたくさん使う。しかし、要素型が整数の行列でもその行列式の値を浮動小数点型に拡張せずに、要素型で求められる。

// matrix.cpp
#include <iostream>
#include <initializer_list>


template <typename T, int N>
class Matrix {
  static_assert(N > 0, "");


public:
  //------------------------------------------
  //
  //  挿入演算子
  //
  //------------------------------------------
  friend std::ostream& operator <<(std::ostream& os, const Matrix& x) {
    os << "[";
    for (int iRow = 0; iRow < N; iRow++) {
      os << "[";
      for (int iColumn = 0; iColumn < N; iColumn++) {
        os << x.elements_[iRow][iColumn];
        if (iColumn != N - 1)
          os << ",";
      }
      os << "]";
      if (iRow != N - 1)
        os << ",";
    }
    os << "]";
    return os;
  }


public:
  //------------------------------------------
  //
  //  コンストラクタ、デストラクタ、コピー、ムーブ
  //
  //------------------------------------------
  Matrix() noexcept(noexcept(T())) {}

  Matrix(const T& a) {
    for (int i = 0; i < N; i++)
      elements_[i][i] = a;
  }

  Matrix(std::initializer_list<std::initializer_list<T>> ll) {
    auto itRow = ll.begin();
    for (int iRow = 0; iRow < N && itRow != ll.end(); iRow++, ++itRow) {
      auto l = *itRow;
      auto itColumn = l.begin();
      for (int iColumn = 0; iColumn < N && itColumn != l.end(); iColumn++, ++itColumn) {
        elements_[iRow][iColumn] = *itColumn;
      }
    }
  }

  Matrix(const Matrix&) noexcept = default;
  Matrix(Matrix&&) noexcept = default;
  ~Matrix() = default;

  Matrix& operator =(const Matrix&) & noexcept = default;
  Matrix& operator =(Matrix&&) & noexcept = default;

  Matrix&& operator =(const Matrix&) && noexcept = delete;
  Matrix&& operator =(Matrix&&) && noexcept = delete;


public:
  //------------------------------------------
  //
  //  要素の参照
  //
  //------------------------------------------
  T& operator ()(int row, int column) {
    return elements_[row][column];
  }

  const T& operator ()(int row, int column) const {
    return elements_[row][column];
  }


public:
  //------------------------------------------
  //
  //  算術単項演算子
  //
  //------------------------------------------
  Matrix operator +() const {
    return Matrix(*this);
  }

  Matrix operator -() const {
    Matrix ret;
    for (int iRow = 0; iRow < N; iRow++) {
      for (int iColumn = 0; iColumn < N; iColumn++) {
        ret.elements_[iRow][iColumn] = -elements_[iRow][iColumn];
      }
    }
    return ret;
  }


public:
  //------------------------------------------
  //
  //  算術代入演算子
  //
  //------------------------------------------
  Matrix& operator +=(const Matrix& other) & {
    for (int iRow = 0; iRow < N; iRow++) {
      for (int iColumn = 0; iColumn < N; iColumn++) {
        elements_[iRow][iColumn] += other.elements_[iRow][iColumn];
      }
    }
    return *this;
  }

  Matrix& operator -=(const Matrix& other) & {
    for (int iRow = 0; iRow < N; iRow++) {
      for (int iColumn = 0; iColumn < N; iColumn++) {
        elements_[iRow][iColumn] -= other.elements_[iRow][iColumn];
      }
    }
    return *this;
  }

  Matrix& operator *=(const T& a) & {
    for (int iRow = 0; iRow < N; iRow++) {
      for (int iColumn = 0; iColumn < N; iColumn++) {
        elements_[iRow][iColumn] *= a;
      }
    }
    return *this;
  }

  Matrix& operator *=(const Matrix& other) & {
    Matrix tmp = *this * other;
    using std::swap;
    swap(*this, tmp);
    return *this;
  }

  Matrix&& operator +=(const Matrix& other) && = delete;
  Matrix&& operator -=(const Matrix& other) && = delete;
  Matrix&& operator *=(const T& a) && = delete;
  Matrix&& operator *=(const Matrix& other) && = delete;


public:
  //------------------------------------------
  //
  //  算術二項演算子
  //
  //------------------------------------------
  friend Matrix operator +(const Matrix& lhs, const Matrix& rhs) {
    Matrix ret(lhs);
    return ret += rhs;
  }

  friend Matrix operator -(const Matrix& lhs, const Matrix& rhs) {
    Matrix ret(lhs);
    return ret -= rhs;
  }

  friend Matrix operator *(const Matrix& lhs, const T& rhs) {
    Matrix ret(lhs);
    return ret *= rhs;
  }

  friend Matrix operator *(const T& lhs, const Matrix& rhs) {
    return rhs *= lhs;
  }

  friend Matrix operator *(const Matrix& lhs, const Matrix& rhs) {
    Matrix ret;
    for (int iRow = 0; iRow < N; iRow++) {
      for (int iColumn = 0; iColumn < N; iColumn++) {
        for (int i = 0; i < N; i++) {
          ret.elements_[iRow][iColumn] += lhs.elements_[iRow][i] * rhs.elements_[i][iColumn];
        }
      }
    }
    return ret;
  }


public:
  //------------------------------------------
  //
  //  論理演算子
  //
  //------------------------------------------
  friend bool operator ==(const Matrix& lhs, const Matrix& rhs) {
    for (int iRow = 0; iRow < N; iRow++) {
      for (int iColumn = 0; iColumn < N; iColumn++) {
        if (lhs.elements_[iRow][iColumn] != rhs.elements_[iRow][iColumn])
          return false;
      }
    }
    return true;
  }

  friend bool operator !=(const Matrix& lhs, const Matrix& rhs) {
    return !(lhs == rhs);
  }


  //------------------------------------------
  //
  // iRow行とiCokumn列を除いた成分を行列にして返す。
  //
  //------------------------------------------
  Matrix<T, N - 1> minor_matrix(int iRow, int iColumn) const {
    static_assert(N > 1, "");

    Matrix<T, N - 1> ret;
    for (int iRR = 0, i = -1; iRR < N; iRR++) {
      if (iRR == iRow) continue;
      i++;
      for (int iCC = 0, j = -1; iCC < N; iCC++) {
        if (iCC == iColumn) continue;
        j++;
        ret(i, j) = elements_[iRR][iCC];
      }
    }
    return ret;
  }


protected:
  T elements_[N][N] = {};

};  // class Matrix


//------------------------------------------
//
// (i, j)成分の余因子を求める。
//
//------------------------------------------
template <typename T, int N> inline
T cofactor(const Matrix<T, N>& m, int iRow, int iColumn) {
  T sign = (iRow + iColumn) % 2 == 0 ? 1 : -1;
  return sign * det(m.minor_matrix(iRow, iColumn));
}


//------------------------------------------
//
// 行列式を求めるときのヘルパー。
// 余因子展開で再帰的(detとcofactorが相互に呼ぶ)に計算する。
// 部分的特殊化で N == 1 の動きを変えるために必要。
// (関数テンプレートは部分的特殊化ができないが、クラステンプレートならできる)
//
//------------------------------------------
template <typename T, int N>
struct Det {
  inline static T det(const Matrix<T, N>& matrix) {
    T ret = 0;
    for (int iColumn = 0; iColumn < N; iColumn++) {
      ret += matrix(0, iColumn) * cofactor(matrix, 0, iColumn);
    }
    return ret;
  }
};


template<typename T>
struct Det<T, 1> {
  inline static T det(const Matrix<T, 1>& matrix) {
    return matrix(0, 0);
  }
};


//------------------------------------------
//
// 行列式を求める
//
//------------------------------------------
template <typename T, int N> inline
T det(const Matrix<T, N>& matrix) {
  return Det<T, N>::det(matrix);
}


//------------------------------------------
//
// 余因子行列を求める
//
//------------------------------------------
template <typename T, int N> inline
Matrix<T, N> cofactor_matrix(const Matrix<T, N>& m) {
  Matrix<T, N> ret;
  for (int iRow = 0; iRow < N; iRow++) {
    for(int iColumn = 0; iColumn < N; iColumn++) {
      //ここで行と列を入れ替えているため、転置したことになっている。
      ret(iRow, iColumn) = cofactor(m, iColumn, iRow);
    }
  }
  return ret;
}


int main()
{
  Matrix<int, 3> x{{1,2,3},{4,5,6},{7,8,9}};
  Matrix<int, 3> y = cofactor_matrix(x);

  std::cout << det(x) << std::endl;
  std::cout << y << std::endl;
  std::cout << x * y << std::endl;

  return 0;
}

実行結果は

0
[[-3,6,-3],[6,-12,6],[-3,6,-3]]
[[0,0,0],[0,0,0],[0,0,0]]

のようになる。

det(x)が0なので、余因子行列が零因子になっていることがわかる。

次回は、転置行列、要素型での除算、逆行列、いくつかのis関数を定義してみたい。
ずいぶんと、ソースが長くなってしまった。