前回は行列式と余因子行列を定義した。
最後の今回は、転置行列、逆行列、要素型との除算を実装する。また、いくつかの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で実装をしてみたい。
