2010-03-17

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():逆行列は存在しません

2010-03-17

C++で作るベクトルクラスその２

programming C++

3次元版を作っていて気づいたこと。
ベクトルクラスを作るときに、x,y,zというメンバを使ってコードを書くか、配列v[i]を利用してコードを書くか迷うところですが、3次元ならx,y,zを使って書いた方が分かりやすいかもしれません。↓

Vector3D.h

class Vector3D{
public:
	//メンバ変数
	float x;
	float y;
	float z;
	//コンストラクタ
	Vector3D();
	Vector3D(float x,float y,float z);
	//代入演算子
	Vector3D& Vector3D::operator=(const Vector3D& v);
	//単項演算子
	Vector3D& operator+=(const Vector3D& v);
	Vector3D& operator-=(const Vector3D& v);
	Vector3D& operator*=(float k);
	Vector3D& operator/=(float k);
	Vector3D operator+()const;
	Vector3D operator-()const;
	//添え字演算子
	float& operator[](int i);
	//比較演算子
	bool operator==(const Vector3D& v ) const;
	bool operator!=(const Vector3D& v ) const;
	//べクトルの長さ
	float norm()const;
	//正規化
	void normalize();
};

//ベクトル演算
//Vector3D+Vector3D
Vector3D operator+(const Vector3D& u,const Vector3D& v);
//Vector3D-Vector3D
Vector3D operator-(const Vector3D& u,const Vector3D& v);
//float*Vector3D
Vector3D operator*(float k,const  Vector3D& v);
//Vector3D*float
Vector3D operator*(const Vector3D& v,float k);
//Vector3D/float
Vector3D operator/(const Vector3D& v,float k);
//内積 Vector3D*Vector3D
float operator*(const Vector3D& u,const Vector3D& v);
//外積 Vector3D%Vector3D
Vector3D operator%(const Vector3D& u,const Vector3D& v);

//2つのベクトルのなす角度
float angle(const Vector3D& u,const Vector3D& v);

//出力
#include <iostream>
inline std::ostream& operator<<(std::ostream& s, const Vector3D& v);

constメンバ関数

定数ベクトルを利用するときにはconstなインスタンスを使うことになります。constなインスタンスが呼び出すことのできるのは次のようなconstメンバ関数だけになります。定数ベクトルに関する演算のパターンを考えるとどのメンバ関数をconst化するべきかが分かります。

	Vector3D operator+()const;
	Vector3D operator-()const;

例えば、定数ベクトルの長さを取得したり、定数ベクトルを比較するときには、メンバの値を書き換えずに済むのでconstを付けても問題はなさそうです。一方、定数ベクトルに*=という演算を施す場合には各要素の値が変わります。したがって、operator*=をconstメンバにしてはいけないことが分かります。

	const Vector3D I1(1,0,0);
	const Vector3D I2(0,1,0);
	cout << +I1 <<endl;		//constメンバ関数の呼び出し
	cout << -I1 <<endl;		//constメンバ関数の呼び出し
	cout << a*I1 <<endl;		//引数がconst
	I1*=3;				//値を書き換えているので、コンパイルエラー、constは付けない
	cout << I2*I1 <<endl;		//constメンバ関数の呼び出し
	if(I1 == a){			//constメンバ関数の呼び出し
 		//処理
	}
	cout<<"angle(I1,I2)="<<angle(I1,I2)<<"\n";	//2つの引数がconst

inline関数

関数を呼び出すとそのたびにオーバーヘッドが生じます。関数の実装部分に"inline"を付けると、関数の呼び出し部分でその関数の処理内容を展開して書き下すことができます。
小さな処理で多く呼ばれるタイプの関数にはinlineを付けると効果的です。
関数の呼び出し部分にその処理を書き下すためには、コンパイル時にヘッダ内の宣言部分で実装内容を知らないといけません。そのため、inline関数の実装部分は、クラス宣言内に書くか、クラス宣言の外であるならヘッダ内で実装を書いておく必要があります。(.cppにinline関数の実装部分を記述するとコンパイルエラー)
inline関数ばかりならVector3D.hに全て定義を書いてしまうのが良いと思います。

比較演算子

ベクトルの比較は全ての要素が同じならtrueを返すようにします。!=は==の逆を返しましょう。

Vector3D.h

inline bool Vector3D::operator==(const Vector3D& v ) const{
    return (x == v.x) && (y == v.y) && (z == v.z);
}
inline bool Vector3D::operator!=(const Vector3D& v ) const{
    return !(*this == v);
}

添え字演算子

要素にアクセスするのに、配列と同じ書き方ができると便利です。範囲外の処理(例えばxを返すなど)も書く必要があります。

//添え字演算子
inline float& Vector3D::operator[](int i){
	if(i == 0){
		return x;
	}
	else if(i == 1){
		return y;
	}
	else if(i == 2){
		return z;
	}
	else{
		return x;
	}
}

外積

3次元ベクトルでは外積の計算がよく使われるので追加してみました。演算子は例えば%を使うといいかもしれません。

//外積 Vector3D%Vector3D
inline Vector3D operator%(const Vector3D& u,const Vector3D& v){//IN:Vector3D Vector3D OUT:Vector3D
	Vector3D w;
	w.x=u.y*v.z-u.z*v.y;
	w.y=u.z*v.x-u.x*v.z;
	w.z=u.x*v.y-u.y*v.x;
	return w;
}

ベクトルに関する関数

一通りベクトルの基本演算を定義することができたので次はベクトルに関する関数を追加してみます。関数を定義するには、ベクトルのメンバにする方法と外部の関数で定義する方法があります。
グローバル関数で定義すると見た目は見やすいですが、グローバル関数が増えすぎると管理が煩雑になってしまいます。クラスのメソッドでは、引数がない関数はともかく2つのベクトルに関する計算は見た目の点でちょっとマイナスになるかもしれません。

メンバ関数として定義

Vector3D u(1,1,1);
Vector3D v(1,2,1);
cout<< u.angle(v) << endl;

グローバル関数として定義

Vector3D u(1,1,1);
Vector3D v(1,2,1);
cout<< angle(u,v) << endl;

結局作成したものはこんな感じになりました↓

Vector3D.h

#ifndef Vector3D_H
#define Vector3D_H

class Vector3D{
public:
	//メンバ変数
	float x;
	float y;
	float z;
	//コンストラクタ
	Vector3D();
	Vector3D(float x,float y,float z);
	//代入演算子
	Vector3D& Vector3D::operator=(const Vector3D& v);
	//単項演算子
	Vector3D& operator+=(const Vector3D& v);
	Vector3D& operator-=(const Vector3D& v);
	Vector3D& operator*=(float k);
	Vector3D& operator/=(float k);
	Vector3D operator+()const;
	Vector3D operator-()const;
	//添え字演算子
	float& operator[](int i);
    //比較演算子
	bool operator==(const Vector3D& v ) const;
	bool operator!=(const Vector3D& v ) const;
	//べクトルの長さ
	float norm()const;
	//正規化
	void normalize();
};
//ベクトル演算
//Vector3D+Vector3D
Vector3D operator+(const Vector3D& u,const Vector3D& v);
//Vector3D-Vector3D
Vector3D operator-(const Vector3D& u,const Vector3D& v);
//float*Vector3D
Vector3D operator*(float k,const  Vector3D& v);
//Vector3D*float
Vector3D operator*(const Vector3D& v,float k);
//Vector3D/float
Vector3D operator/(const Vector3D& v,float k);
//内積 Vector3D*Vector3D
float operator*(const Vector3D& u,const Vector3D& v);
//外積 Vector3D%Vector3D
Vector3D operator%(const Vector3D& u,const Vector3D& v);
//2つのベクトルのなす角度
float angle(const Vector3D& u,const Vector3D& v);
//出力
#include <iostream>
inline std::ostream& operator<<(std::ostream& s, const Vector3D& v);
//*----------------------メンバ関数の実装--------------------------*//
#include <cmath>
//コンストラクタ
inline Vector3D::Vector3D(){ x = y = z = 0; }
inline Vector3D::Vector3D(float x,float y,float z){
	this->x=x;		this->y=y;		this->z=z;
}
//代入演算子
inline Vector3D& Vector3D::operator=(const Vector3D& v){
	this->x=v.x;	this->y=v.y;	this->z=v.z;
	return *this;
}
//単項演算子
inline Vector3D& Vector3D::operator+=(const Vector3D& v){
	 this->x += v.x;	this->y += v.y;		this->z += v.z;
	 return *this;
}
inline Vector3D& Vector3D::operator-=(const Vector3D& v){
	 this->x -= v.x;	this->y -= v.y;		this->z -= v.z;
	 return *this;
}
inline Vector3D& Vector3D::operator*=(float k){
	 this->x *= k;		this->y *= k;		this->z *= k;
	 return *this;
}
inline Vector3D& Vector3D::operator/=(float k){
	this->x /= k;		this->y /= k;		this->z /= k;
	return *this;
}
inline Vector3D Vector3D::operator+()const{		//+Vector3D
	return *this;
}
inline Vector3D Vector3D::operator-()const{		//-Vector3D
	return Vector3D(-x,-y,-z);
}
//添え字演算子
inline float& Vector3D::operator[](int i){
	if(i == 0){
		return x;
	}
	else if(i == 1){
		return y;
	}
	else if(i == 2){
		return z;
	}
	else{
		return x;
	}
}
//比較演算子
inline bool Vector3D::operator==(const Vector3D& v ) const{
    return (x == v.x) && (y == v.y) && (z == v.z);
}
inline bool Vector3D::operator!=(const Vector3D& v ) const{
    return !(*this == v);
}
//べクトルの長さ
inline float Vector3D::norm()const{
	return pow(x*x+y*y+z*z,0.5f);
}
//正規化
inline void Vector3D::normalize(){
	*this /= norm();
}
//*----------------------グローバル関数の実装--------------------------*//
//二項演算子の定義
//Vector3D+Vector3D
inline Vector3D operator+(const Vector3D& u,const Vector3D& v){
	Vector3D w;
	w.x=u.x+v.x;
	w.y=u.y+v.y;
	w.z=u.z+v.z;
	return w;
}
//Vector3D-Vector3D
inline Vector3D operator-(const Vector3D& u,const Vector3D& v){
	Vector3D w;
	w.x=u.x-v.x;
	w.y=u.y-v.y;
	w.z=u.z-v.z;
	return w;
}
//float*Vector3D
inline Vector3D operator*(float k,const  Vector3D& v){
	return Vector3D(k*v.x,k*v.y,k*v.z);
}
//Vector3D*float
inline Vector3D operator*(const Vector3D& v,float k){
	return Vector3D(v.x*k,v.y*k,v.z*k);
}
//Vector3D/float
inline Vector3D operator/(const Vector3D& v,float k){
	return Vector3D(v.x/k,v.y/k,v.z/k);
}
//内積 Vector3D*Vector3D
inline float operator*(const Vector3D& u,const Vector3D& v){
	return u.x*v.x+u.y*v.y+u.z*v.z;
}
//外積 Vector3D%Vector3D
inline Vector3D operator%(const Vector3D& u,const Vector3D& v){
	Vector3D w;
	w.x=u.y*v.z-u.z*v.y;
	w.y=u.z*v.x-u.x*v.z;
	w.z=u.x*v.y-u.y*v.x;
	return w;
}
//画面への表示
#include <iostream>
inline std::ostream& operator<<(std::ostream& s, const Vector3D& v){
	return s <<'('<<v.x<<","<<v.y<<","<<v.z<<')';
}

#define PI 3.1415926535
//2つのベクトルのなす角
inline float angle(const Vector3D& u,const Vector3D& v){
	float cos =u*v/(u.norm()*v.norm());
	return float(acos(cos)/PI*180);
}

#endif

作成したクラスでテストしたのがこちら。

int main()
{
	Vector3D u(1,2,3);
	Vector3D v(4,1,0);
	u[2]=3;
	cout <<"u[2]="<<u[2]<<endl;	//3
	
	cout<<"u="<<u<<endl;		//(1,2,3)
	cout<<"v="<<v<<endl;		//(4,1,0)

	cout<<"-v="<<-v<<endl;		//(-4,-1,0)
	cout<<"u+v="<<u+v<<endl;	//(5,3,3)
	cout<<"u-v="<<u-v<<endl;	//(-3,1,3)
	cout<<"u*v="<<u*v<<endl;	//6
	cout<<"u%v="<<u%v<<endl;		//(-3,12,-7)

	u+=Vector3D(1,1,1);
	cout<<"u+(1,1,1)="<<u<<endl;	//u=(2,3,4)
	u-=Vector3D(1,1,1);
	cout<<"u-(1,1,1)="<<u<<endl;	//u=(1,2,3)
	u*=3;
	cout<<"u*3="<<u<<endl;			//u=(3,6,9)
	u/=3;
	cout<<"u/3="<<u<<endl;			//u=(1,2,3)

	cout <<"norm(u)^2="<<u.norm()*u.norm()<<endl;	//14	
	Vector3D a(1,2,3);
	Vector3D b;
	b = a;
	cout << (a == b) <<endl;		//1 true
	cout << (a != b) <<endl;		//0 false

	const Vector3D I1(1,0,0);
	const Vector3D I2(0,1,0);
	cout<<"angle(I1,I2)="<<angle(I1,I2)<<endl;		//(1,2,3)
	return 0;
}

u[2]=3
u=(1,2,3)
v=(4,1,0)
-v=(-4,-1,-0)
u+v=(5,3,3)
u-v=(-3,1,3)
u*v=6
u%v=(-3,12,-7)
u+(1,1,1)=(2,3,4)
u-(1,1,1)=(1,2,3)
u*3=(3,6,9)
u/3=(1,2,3)
norm(u)^2=14
1
0
angle(I1,I2)=90