C++で作る行列クラス その１

programming C++

3×3の行列クラスを作ってみます。
こちらも、できるだけ数式に近い形で扱いたいということで演算子のオーバーロードを使って実装してみます。行列クラスをこんな感じで使うのが目標。

	Matrix3D A(
	1,	0,	0,
	0,	1,	0,
	0,	0,	1
	);
	Matrix3D B(
	2,	3,	4,
	3,	5,	6,
	1,	4,	5
	);
	A.a=2;	A.b=2;	A.c=3;
	A.d=1;	A.e=4;	A.f=0;
	A.g=0;	A.h=9;	A.i=7;

	A[1][1]=2;	A[2][0]=-2;
	A(1, 1)=3;	A(2, 0)=-5;

	cout <<A <<endl;
	cout <<"A+B=\n"<< A+B <<endl;
	A=2*A;
	cout << inverse(A)*A <<endl;
メンバへのアクセスの方法

A.a,A.b,...cdeefghiとA[i][j]の両方でアクセスできるようにするには、同じ領域に違う名前を付けなければならないのでunionを使います。一方は9個の変数、もう一方は二次元配列で定義します。
配列の個数を予め定めないといけないので、static constで行と列の値を決めておきます。今回は正方行列ですが、行と列でサイズがことなる場合も考えるため別々に定義しておくことにします。

class Matrix3D{
public:
	static const int ROW = 3;
	static const int COL = 3;
	union{
		struct{
			float	a,b,c,
				d,e,f,
				g,h,i;
		};
		float val[ROW][COL];
	};

コンストラクタ

行列の要素をコンストラクタで直接設定できると便利なので、9個の引数から地道にメンバに代入するようにします。他のクラスのメンバにすることを想定して、引数がないコンストラクタも同時に作成しておきます。

	Matrix3D(float x00, float x01, float x02,
		  float x10, float x11, float x12,
		  float x20, float x21, float x22){
		  val[0][0]=x00;val[0][1]=x01;val[0][2]=x02;
		  val[1][0]=x10;val[1][1]=x11;val[1][2]=x12;
		  val[2][0]=x20;val[2][1]=x21;val[2][2]=x22;
	}
	Matrix3D(){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				val[i][j]=0;
			}
		}
	}

添え字演算子

要素にアクセスするために2通りの演算子を考えてみます。一つ目はよく利用する[i][j]タイプで、この場合はoperator[]がfloat*を返すようにすればいつも通り使えることになります。もうひとつは(i,j)タイプの演算子ですが、こちらはval[i][j]への参照を返すことで、要素の値を読み書きすることができるようになります。

	float* operator[](int i){
		return val[i];
	}
	float& operator()(int i,int j){
		return val[i][j];
	}

代入演算子

代入のときはそのまま要素をコピーして、*thisを返すとOKです。

	Matrix3D& operator=(const Matrix3D& A){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
			 this->val[i][j]=A.val[i][j];
			}
		}
		return *this;
	}

単項演算子の定義

まず"+A"や"-A"タイプの演算を定義します。+の場合は自分自身を、-は全ての要素の符号を反転させた行列を返します。

	Matrix3D operator+(){return *this;}	//+Matrix3D
	Matrix3D operator-(){		//-Matrix3D A;
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				A.val[i][j]=-val[i][j];
			}
		}
		return A;
	}

次に、"+=","-=","*=","/="を実装します。

	Matrix3D& operator+=(const Matrix3D& A){
		for(int i=0;i<ROW;i++){
		   for(int j=0;j<COL;j++){
				val[i][j]+=A.val[i][j];
		   }
		}
		return *this;// IN:Matrix  OUT:this
	}
	Matrix3D& operator*=(float k){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				val[i][j]*=k;
			}
		}
		return *this;// IN:scalar  OUT:this
	}

比較演算子

これもベクトルの場合と同様に全ての要素が同じかどうかを判定します。!=のときは==の判定結果の逆を返します。

	bool operator==(const Matrix3D& A ) const{
		bool result=true;
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				result &= (val[i][j]==A.val[i][j]);
			}
		}
		return result;
	}
	bool operator!=(const Matrix3D& A ) const{
		return !((*this) == A);
	}

二項演算子の定義

まずは行列と行列の演算として"+,-,*"の3つを実装します。

//Matrix3D+Matrix3D
inline Matrix3D operator+(const Matrix3D& A,const Matrix3D& B){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D C;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			C.val[i][j]=A.val[i][j]+B.val[i][j];
		}
	}
	return C;
}
//product C=AB//Matrix3D*Matrix3D
inline Matrix3D operator*(const Matrix3D& A,const Matrix3D& B){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D C;//C=O
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			for(int k=0;k<Matrix3D::ROW;k++)C.val[i][j]+=A.val[i][k]*B.val[k][j];
		}
	}
	return C;
}//operator*

積の演算子"*"では行列とベクトルの演算を定義することもできます。

//product w=Av//Matrix3D*Vector3D
inline Vector3D operator*(const Matrix3D& A,const Vector3D& v){//IN:Vector3D Matrix3D OUT:Vector3D 
	Vector3D w;//w=O

	for(int i=0;i<Matrix3D::ROW;i++){
		w.val[i] = 0;
		for(int j=0;j<Matrix3D::COL;j++){
			w.val[i] += v.val[j]*A.val[i][j];
		}
	}
	return w;
}//operator*

最後に行列と数値の演算として、"float*Matrix","Matrix*float","Matrix/float"を実装すれば演算子は完成です。

画面への出力

コンソールに行列の値を表示するために演算子"<<"を定義してみます。

#include <iostream>
inline std::ostream& operator<<(std::ostream& s, const Matrix3D& A){
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++)s <<A.val[i][j]<<"\t";
			s <<"\n";
		}
	return s;
}

行列に関する関数

ここまでで一通りの行列演算はできるようになったはずです。あとは行列に関する関数をどれだけ追加するかという問題になります。3×3の行列なのでトレースや行列式、逆行列を計算するには直に書いてしまった方が分かりやすいかもしれません。
トレース、行列式、逆行列、転置を追加してひとまず完成↓
ベクトルクラスはC++で作るベクトルクラス その２ - white wheelsのメモで作ったものを流用

#ifndef Matrix3D_H
#define Matrix3D_H
#include <stdexcept>
#include "Vector3D.h"

class Matrix3D{
public:
	static const int ROW = 3;
	static const int COL = 3;
	//メンバ変数
	union{
		struct{
			float	a,b,c,
				d,e,f,
				g,h,i;
		};
		float val[ROW][COL];
	};
	//コンストラクタ
	Matrix3D(){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				val[i][j]=0;
			}
		}
	}
	Matrix3D(float x00, float x01, float x02,
		  float x10, float x11, float x12,
		  float x20, float x21, float x22){
		  val[0][0]=x00;val[0][1]=x01;val[0][2]=x02;
		  val[1][0]=x10;val[1][1]=x11;val[1][2]=x12;
		  val[2][0]=x20;val[2][1]=x21;val[2][2]=x22;
	}
	//代入演算子
	Matrix3D& operator=(const Matrix3D& A){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
			 this->val[i][j]=A.val[i][j];
			}
		}
		return *this;
	}
	Matrix3D operator+(){return *this;}	//+Matrix3D
	Matrix3D operator-(){				//-Matrix3D
		Matrix3D A;
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				A.val[i][j]=-val[i][j];
			}
		}
		return A;
	}
	// +=
	Matrix3D& operator+=(const Matrix3D& A){
		for(int i=0;i<ROW;i++){
		   for(int j=0;j<COL;j++){
				val[i][j]+=A.val[i][j];
		   }
		}
		return *this;// IN:Matrix  OUT:this
	}
	// -=
	Matrix3D& operator-=(const Matrix3D& A){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				val[i][j]-=A.val[i][j];
			}
		}
		return *this;// IN:Matrix  OUT:this
	}
	// *=
	Matrix3D& operator*=(float k){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				val[i][j]*=k;
			}
		}
		return *this;// IN:scalar  OUT:this
	}
	// /=
	Matrix3D& operator/=(float k){
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				val[i][j]/=k;
			}
		}
		return *this;// IN:scalar  OUT:this
	}
	//比較演算子
	bool operator==(const Matrix3D& A ) const{
		bool result=true;
		for(int i=0;i<ROW;i++){
			for(int j=0;j<COL;j++){
				result &= (val[i][j]==A.val[i][j]);
			}
		}
		return result;
	}
	bool operator!=(const Matrix3D& A ) const{
		return !((*this) == A);
	}
	//添え字演算子
	float* operator[](int i){
		return val[i];
	}
	float& operator()(int i,int j){
		return val[i][j];
	}

};
//Matrix3D+Matrix3D
inline Matrix3D operator+(const Matrix3D& A,const Matrix3D& B){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D C;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			C.val[i][j]=A.val[i][j]+B.val[i][j];
		}
	}
	return C;
}
//Matrix3D-Matrix3D
inline Matrix3D operator-(const Matrix3D& A,const Matrix3D& B){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D C;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			C.val[i][j]=A.val[i][j]-B.val[i][j];
		}
	}
	return C;
}
//float*Matrix3D
inline Matrix3D operator*(float k,const  Matrix3D& A){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D B;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			B.val[i][j]=A.val[i][j]*k;
		}
	}
	return B;
}
//Matrix3D*float
inline Matrix3D operator*(const Matrix3D& A,float k){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D B;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			B.val[i][j]=A.val[i][j]*k;
		}
	}
	return B;
}
//Matrix3D/float
inline Matrix3D operator/(const Matrix3D& A,float k){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D B;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			B.val[i][j]=A.val[i][j]/k;
		}
	}
	return B;
}
//product C=AB//Matrix3D*Matrix3D
inline Matrix3D operator*(const Matrix3D& A,const Matrix3D& B){//IN:Matrix3D Matrix3D OUT:Matrix3D
	Matrix3D C;//C=O
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			for(int k=0;k<Matrix3D::ROW;k++)C.val[i][j]+=A.val[i][k]*B.val[k][j];
		}
	}
	return C;
}//operator*
//product w=Av//Matrix3D*Vector3D
inline Vector3D operator*(const Matrix3D& A,const Vector3D& v){//IN:Vector3D Matrix3D OUT:Vector3D 
	Vector3D w;//w=O

	for(int i=0;i<Matrix3D::ROW;i++){
		w.val[i] = 0;
		for(int j=0;j<Matrix3D::COL;j++){
			w.val[i] += v.val[j]*A.val[i][j];
		}
	}
	return w;
}//operator*
//トレースを取得する
inline float trace(const Matrix3D& A){
	float tr=0;
	for(int i=0;i<Matrix3D::ROW;i++){
		tr+=A.val[i][i];
	}
	return tr;
}
//転置行列を取得する
inline Matrix3D transpose(const Matrix3D& A){
	Matrix3D AT;
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++){
			AT.val[i][j]=A.val[j][i];
		}
	}
	return AT;
}
//行列式を取得する
inline float det(const Matrix3D& A){
	//Sarrusの方法
	float d = 0;
	d+=A.val[0][0]*A.val[1][1]*A.val[2][2];
	d+=A.val[0][1]*A.val[1][2]*A.val[2][0];
	d+=A.val[0][2]*A.val[1][0]*A.val[2][1];
	d-=A.val[0][2]*A.val[1][1]*A.val[2][0];
	d-=A.val[0][0]*A.val[1][2]*A.val[2][1];
	d-=A.val[0][1]*A.val[1][0]*A.val[2][2];
	return d;
}
//逆行列を取得する
inline Matrix3D inverse(const Matrix3D& A){
	Matrix3D Ai;
	float d=det(A);

	if(d == 0){
		throw std::overflow_error("inverse():逆行列は存在しません");//error処理
	}
	//余因子展開,cramerの公式
	Ai.val[0][0]=(A.val[1][1]*A.val[2][2]-A.val[1][2]*A.val[2][1])/d;
	Ai.val[0][1]=-(A.val[0][1]*A.val[2][2]-A.val[0][2]*A.val[2][1])/d;
	Ai.val[0][2]=(A.val[0][1]*A.val[1][2]-A.val[1][1]*A.val[0][2])/d;

	Ai.val[1][0]=-(A.val[1][0]*A.val[2][2]-A.val[1][2]*A.val[2][0])/d;
	Ai.val[1][1]=(A.val[0][0]*A.val[2][2]-A.val[0][2]*A.val[2][0])/d;
	Ai.val[1][2]=-(A.val[0][0]*A.val[1][2]-A.val[1][0]*A.val[0][2])/d;

	Ai.val[2][0]=(A.val[1][0]*A.val[2][1]-A.val[1][1]*A.val[2][0])/d;
	Ai.val[2][1]=-(A.val[0][0]*A.val[2][1]-A.val[0][1]*A.val[2][0])/d;
	Ai.val[2][2]=(A.val[0][0]*A.val[1][1]-A.val[0][1]*A.val[1][0])/d;
	return Ai;
}
//出力
#include <iostream>
inline std::ostream& operator<<(std::ostream& s, const Matrix3D& A){
	for(int i=0;i<Matrix3D::ROW;i++){
		for(int j=0;j<Matrix3D::COL;j++)s <<A.val[i][j]<<"\t";
			s <<"\n";
		}
	return s;
}
#endif

作成した行列クラスの利用例です。

#include "Matrix3D.h"
#include <iostream>
using namespace std;
int main(){
	Matrix3D A(
	1,	0,	0,
	0,	1,	0,
	0,	0,	1
	);

	A.a=2;	A.b=2;	A.c=3;
	A.d=1;	A.e=4;	A.f=0;
	A.g=0;	A.h=-1;	A.i=-2;

	A[1][1]=2;	A[2][0]=-2;
	A(1, 1)=3;	A(2, 0)=-5;
	cout <<"A=\n"<< A <<endl;

	Matrix3D B(
	2,	3,	1,
	3,	1,	0,
	1,	-4,	1
	);
	cout <<"B=\n"<< B <<endl;

	Matrix3D C(
	-1,	2,	3,
	-1,	-1,	1,
	2,	1,	2
	);

	cout <<"C=\n"<< C <<endl;
	cout <<"A+B=\n"<< A+B <<endl;
	cout <<"A-B=\n"<< A-B <<endl;
	C=2*C;
	cout <<"C=2*C:C=\n"<< C <<endl;
	cout <<"C/2=\n"<< C/2 <<endl;

	Matrix3D D(
	-2,	1,	1,
	0,	2,	1,
	1,	4,	3
	);
	cout <<"CD=\n"<<C<<"*\n"<<D<<"=\n"<<C*D<<endl;
	cout <<"trC="<< trace(C) <<endl;
	cout <<"transposeC=\n"<< transpose(C) <<endl;
	cout <<"detC="<< det(C) <<endl;

	Vector3D p=Vector3D(1,2,3);
	cout <<"p="<<p<<endl;
	cout <<"Cp="<<C*p<<endl;

	Matrix3D E(
	3,	2,	1,
	1,	3,	1,
	2,	2,	1
	);
	cout <<"E^-1=\n"<<inverse(E)*E<<endl;
	Matrix3D P(
	-1,	2,	1,
	0,	3,	1,
	0,	2,	5
	);
	Matrix3D Q(
	-1,	2,	1,
	0,	3,	1,
	0,	2,	5
	);

	cout <<"P==Q:" <<(P==Q)<<endl;
	cout <<"P!=Q:" <<(P!=Q)<<endl;

	Matrix3D F;
	try{
	//F=F/0;
	cout <<"F^-1=\n"<<inverse(F)<<endl;
	}
	catch( exception& e )
	{
	cout <<"error:"<<e.what()<<endl;
	}
  return 0;
}


A=
2 2 3
1 3 0
-5 -1 -2

B=
2 3 1
3 1 0
1 -4 1

C=
-1 2 3
-1 -1 1
2 1 2

A+B=
4 5 4
4 4 0
-4 -5 -1

A-B=
0 -1 2
-2 2 0
-6 3 -3

C=2*C:C=
-2 4 6
-2 -2 2
4 2 4

C/2=
-1 2 3
-1 -1 1
2 1 2

CD=
-2 4 6
-2 -2 2
4 2 4
*
-2 1 1
0 2 1
1 4 3
=
10 30 20
6 2 2
-4 24 18

trC=0
transposeC=
-2 -2 4
4 -2 2
6 2 4

detC=112
p=(1,2,3)
Cp=(24,0,20)
E^-1=
1 0 0
0 1 0
0 0 1

P==Q:1
P!=Q:0
error:inverse():逆行列は存在しません