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

前回は行列式と余因子行列を定義した。
最後の今回は、転置行列、逆行列、要素型との除算を実装する。また、いくつかのis関数を実装する。
今回のプログラムの実行にはBoostライブラリを必要とする。

// matrix.cpp
#include <iostream>
#include <initializer_list>
#include <stdexcept>
#include <boost/rational.hpp>



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 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) && = delete;
  Matrix&& operator -=(const Matrix& other) && = delete;
  Matrix&& operator *=(const T& a) && = delete;
  Matrix&& operator *=(const Matrix& other) && = delete;
  Matrix&& operator /=(const T& a) && = 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;
  }

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


public:
  //------------------------------------------
  //
  // その他のメンバ関数
  //
  //------------------------------------------

  // 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;
  }

  // 零行列ならtrue
  bool is_zero() const {
    for (int iRow = 0; iRow < N; iRow++)
      for(int iColumn = 0; iColumn < N; iColumn++) {
        if (elements_[iRow][iColumn] != 0)
          return false;
      }
    return true;
  }

  // 対角行列ならtrue
  bool is_diag() const {
    for (int iRow = 0; iRow < N; iRow++)
      for(int iColumn = 0; iColumn < N; iColumn++) {
        if (iRow == iColumn)
          continue;
        if (elements_[iRow][iColumn] != 0)
          return false;
      }
    return true;
  }

  // 零因子ならtrue
  bool is_zerodevisor() const {
    return det(*this) == static_cast(0);
  }


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;
}


//------------------------------------------
//
// 転置行列を求める
//
//------------------------------------------
template <typename T, int N> inline
Matrix<T, N> t(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) = m(iColumn, iRow);
    }
  }
  return ret;
}


//------------------------------------------
//
// 逆行列が存在しない場合の例外
//
//------------------------------------------
class ZeroDivisorException : public std::exception {
public:
  const char* what() const noexcept {
    return "Zero divisor cannot create inverse matrix.";
  }
};


//------------------------------------------
//
// 逆行列を求める
//
//------------------------------------------
template <typename T, int N> inline
Matrix<T, N> inverse(const Matrix<T, N>& m) {
  if (det(m) == 0) {
    throw ZeroDivisorException();
  }
  T d = det(m);
  Matrix<T, N> co_mat = cofactor_matrix(m);
  return co_mat / d;
}


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

  try {
    Matrix<int, 3> z = inverse(x);
  } catch (ZeroDivisorException& e) {
    std::cout << e.what() << std::endl;
  }

  std::cout << x / 9 << std::endl << std::endl;


  using boost::rational;
  using M = Matrix<rational<int>, 3>;

  M a{{1,0,4},{3,-1,0},{-2,1,-1}};

  try {
    M b = inverse(a);
    std::cout << b << std::endl;
    std::cout << a * b << std::endl;
  } catch (ZeroDivisorException& e) {
    std::cout << e.what() << std::endl;
  }
  return 0;
}

実行結果は

Zero divisor cannot create inverse matrix. //零因子の逆行列は存在しない。
[[0,0,0],[0,0,0],[0,0,1]] //要素型がintの場合、intとintの除算が反映される。商が整数になり余りが切り捨てられる実装のため注意が必要。

[[1/5,4/5,4/5],[3/5,7/5,12/5],[1/5,-1/5,-1/5]] //要素型を分数にすることもできる。
[[1/1,0/1,0/1],[0/1,1/1,0/1],[0/1,0/1,1/1]]

のようになる。


気が向いたら、constexprや、pImplで実装をしてみたい。