fork download

namespace h3d {
    template <unsigned int N> class Matrix;
	template <unsigned int N> class Vector;
	using FLOAT = double;
	class Test;
	class TestSet;
	class Test1;
	class Test2;
	class Test3;
};
 
#include <iostream>
 
namespace h3d {
	using std::ostream;
 
	/**
	 * N次元ベクトルクラス。
	 * @param N ベクトルの次元数。
	 */
	template <unsigned int N> 
	class Vector {
	public:
		/** 成分(Component)の配列。N個。 */
		FLOAT c[N];
 
		/**
		 * 各成分をゼロクリアする。
		 */
		inline void clear();
 
		/**
		 * 次元を一つ拡張(EXTend)する。
		 * @param c 追加する成分。
		 * @return 拡張したベクトル。
		 */
		inline Vector<N + 1> ext(const FLOAT& c = 1.0) const;
 
		/**
		 * 内積(INner PROduct)を計算する。
		 * @param rhs 右側のベクトル。
		 * @return 計算した内積。
		 */
		inline FLOAT inpro(const Vector& rhs) const;
 
		/**
		 * 長さを取得する。
		 * @return 取得した長さ。
		 */
		inline double norm() const;
 
		/**
		 * 正規化して、長さを1にする。
		 * @return 正規化したベクトル。
		 */
		inline Vector normalize() const;
 
		/**
		 * 二つのベクトルが等しくないかどうかを判定する。
		 * @param rhs 右側のベクトル。
		 * @return 等しくないならtrue, 等しいならfalse。
		 */
		inline bool operator !=(const Vector& rhs) const;
 
		/**
		 * スカラー値で乗算する。
		* @param rhs 乗数となる右側のスカラー値。
		 * @return 乗算したベクトル。
		 */
		inline Vector operator *(const FLOAT& rhs) const;
 
		/**
		 * スカラー値で乗算して、代入する。
		 * @param rhs 乗数となる右側のスカラー値。
		 * @return 乗算して、代入したベクトル。
		 */
		inline Vector& operator *=(const FLOAT& rhs);
 
		/**
		 * ベクトルを加算する。
		 * @param rhs 加数となる右側のベクトル。
		 * @return 加算したベクトル。
		 */
		inline Vector operator +(const Vector& rhs) const;
 
		/**
		 * ベクトルを加算して、代入する。
		 * @param rhs 加数となる右側のベクトル。
		 * @return 加算して、代入したベクトル。
		 */
		inline Vector& operator +=(const Vector& rhs);
 
		/**
		 * 各成分の符号を反転する。
		 * @return 符号を反転したベクトル。
		 */
		inline Vector operator -() const;
 
		/**
		 * ベクトルを減算する。
		 * @param rhs 減数となる右側のベクトル。
		 * @return 減算したベクトル。
		 */
		inline Vector operator -(const Vector& rhs) const;
 
		/**
		 * ベクトルを減算して、代入する。
		 * @param rhs 減数となる右側のベクトル。
		 * @return 減算して、代入したベクトル。
		 */
		inline Vector& operator -=(const Vector& rhs);
 
		/**
		 * スカラー値で除算する。
		 * @param rhs 除数となる右側のスカラー値。
		 * @return 除算したベクトル。
		 */
		inline Vector operator /(const FLOAT& rhs) const;
 
		/**
		 * スカラー値で除算して、代入する。
		 * @param rhs 除数となる右側のスカラー値。
		 * @return 除算して、代入したベクトル。
		 */
		inline Vector& operator /=(const FLOAT& rhs);
 
		/**
		 * 二つのベクトルが等しいかどうかを判定する。
		 * @param rhs 右側のベクトル
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const Vector& rhs) const;
 
		/**
		 * 3次元ベクトルの外積(OUTer PROduct)を計算する。
		 * @param rhs 右側のベクトル。
		 * @return 計算した外積。
		 */
		inline Vector outpro(const Vector& rhs) const;
 
		/**
		 * 次元を一つ切り縮める(TRUNCate)。
		 * @return 切り縮めたベクトル。
		 */
		inline Vector<N - 1> trunc() const;
 
		/**
		 * ゼロベクトルを取得する。
		 * @return 取得したゼロベクトル。
		 */
		static constexpr Vector ZERO();
	};
 
	/**
	 * ベクトルを出力する。
	 * @param N ベクトルの次元数。
	 * @param out 出力先となるストリーム。
	 * @param vec 出力するベクトル。
	 * @return 出力したストリーム。
	 */
	template <unsigned int N> 
	inline ostream& operator <<(ostream& out, const Vector<N>& vec);
 
	template <> 
	constexpr Vector<2> Vector<2>::ZERO() {
		return Vector<2>{0.0, 0.0};
	}
 
	template <> 
	constexpr Vector<3> Vector<3>::ZERO() {
		return Vector<3>{0.0, 0.0, 0.0};
	}
 
	template <> 
	constexpr Vector<4> Vector<4>::ZERO() {
		return Vector<4>{0.0, 0.0, 0.0, 0.0};
	}
};
 
namespace h3d {
	/** 原点(Origin PoinT)。 */
	constexpr Vector<3> O_PT		{0.0,	0.0,	0.0};
	/** X軸。 */
	constexpr Vector<3> X_AXIS		{1.0,	0.0,	0.0};
	/** Y軸。 */
	constexpr Vector<3> Y_AXIS		{0.0,	1.0,	0.0};
	/** Z軸。 */
	constexpr Vector<3> Z_AXIS		{0.0,	0.0,	1.0};
};
 
namespace h3d {
	/**
	 * 拡大・縮小行列を作る。
	 * @param coef 各座標の係数(COEFficient)。
	 * coef > 1.0 なら拡大, coef < 1.0 なら縮小。
	 * @return 作った拡大・縮小行列。
	 */
	Matrix<3> scaling(const Vector<3>& coef);
 
	/**
	 * 回転行列を作る。
	 * @param axis この軸の周りを回転する。
	 * @param rad 回転する弧度(RADian)。
	 * @return 作った回転行列。
	 */
	Matrix<3> rotation(const Vector<3>& axis, const FLOAT& rad);
 
	/**
	 * 平行移動行列を作る。
	 * @param os 変換後の位置から変換前の位置を引いた差分(OffSet)。
	 * @return 作った平行移動行列。
	 */
	Matrix<4> translation(const Vector<3>& os);
 
	/**
	 * ビュー行列を作る。
	 * ビュー行列は、ワールド空間からカメラ空間への変換を行う。
	 * @param eye_pos 目の位置(POSition)。
	 * @param cent_pt 視界の中央(CENTer)にある点(PoinT)。
	 * @param up_dir 上方向(UP DIRection)。省略ならY軸。
	 * @return 作ったビュー行列。
	 */
	Matrix<4> view(const Vector<3>& eye_pos, const Vector<3>& cent_pt, const Vector<3>& up_dir = Y_AXIS);
 
	/**
	 * 透視射影(PERSpective PROJection)行列を作る。
	 * 透視射影行列は、カメラ空間からクリップ空間への変換を行う。
	 * @param hori_fov 水平視野角(HORIzontal Field Of View)。弧度。
	 * @param asp_rat アスペクト比(ASPect RATio)。
	 * @param near_dist 近平面までの距離(DISTance)。
	 * @param far_dist 遠平面までの距離(DISTance)。
	 */
	Matrix<4> persProj(const double& hori_fov, const double& asp_rat, const double& near_dist, const double& far_dist);
};
 
namespace h3d {
	/**
	 * 角度(ANGle)を弧度(RADian)に変換する。
	 * @param ang 変換する角度(ANGle)。
	 * @return 変換した弧度。
	 */
	inline FLOAT ang2rad(const FLOAT& ang);
};
 
#include <iostream>
 
namespace h3d {
	using std::ostream;
	using std::string;
 
	/**
	 * N次元正方行列クラス。
	 * @param N 行列の次元数。
	 */
	template <unsigned int N> 
	class Matrix {
	public:
		/** 成分(Component)の配列。N×N個。 */
		FLOAT c[N][N];
 
		/**
		 * 各成分をゼロクリアする。
		 */
		inline void clear();
 
		/**
		 * 行列式(DETerminant)を計算する。
		 * @return 計算した行列式。
		 */
		inline FLOAT det() const;
 
		/**
		 * 次元を一つ拡張(EXTend)する。
		 * @param mat この行列から最後の行と列の成分をコピーする。
		 * @return 拡張した行列。
		 */
		inline Matrix<N + 1> ext(const Matrix<N + 1>& mat = Matrix<N + 1>::IDENTITY()) const;
 
		/**
		 * 単位行列を取得する。
		 * @return 単位行列。
		 */
		static constexpr Matrix IDENTITY();
 
		/**
		 * 逆行列(INVerse matrix)を計算する。
		 * @return 計算した逆行列。
		 * @throw string 行列が正則ではないならメッセージをスローする。
		 */
		inline Matrix inv() const throw(string);
 
		/**
		 * 二つの行列が等しくないかどうかを判定する。
		 * @param rhs 右側の行列。
		 * @return 等しくないならtrue, 等しいならfalse。
		 */
		inline bool operator !=(const Matrix& rhs) const;
 
		/**
		 * 行列で乗算する。
		 * @param rhs 乗数となる右側の行列。
		 * @return 乗算した行列。
		 */
		inline Matrix operator *(const Matrix& rhs) const;
 
		/**
		 * ベクトルとの積を計算する。
		 * @param rhs 乗数となる右側のベクトル。
		 * @return 計算したベクトル。
		 */
		inline Vector<N> operator *(const Vector<N>& rhs) const;
 
		/**
		 * スカラー値で乗算する。
		 * @param rhs 乗数となる右側のスカラー値。
		 * @return 乗算した行列。
		 */
		inline Matrix operator *(const FLOAT& rhs) const;
 
		/**
		 * 行列で乗算して、代入する。
		 * @param rhs 乗数となる右側の行列。
		 * @return 乗算して、代入した行列。
		 */
		inline Matrix operator *=(const Matrix& rhs);
 
		/**
		 * スカラー値で乗算して、代入する。
		 * @param rhs 乗数となる右側のスカラー値。
		 * @return 乗算して、代入した行列。
		 */
		inline Matrix operator *=(const FLOAT& rhs);
 
		/**
		 * 行列を加算する。
		 * @param rhs 加数となる右側の行列。
		 * @return 加算した行列。
		 */
		inline Matrix operator +(const Matrix& rhs) const;
 
		/**
		 * 行列を加算して、代入する。
		 * @param rhs 加数となる右側の行列。
		 * @return 加算して、代入した行列。
		 */
		inline Matrix operator +=(const Matrix& rhs);
 
		/**
		 * すべての成分の符号を反転する。
		 * @return 符号を反転した行列。
		 */
		inline Matrix operator -() const;
 
		/**
		 * 行列を減算する。
		 * @param rhs 減数となる右側の行列。
		 * @return 減算した行列。
		 */
		inline Matrix operator -(const Matrix& rhs) const;
 
		/**
		 * 行列を減算して、代入する。
		 * @param rhs 減数となる右側の行列。
		 * @return 減算して、代入した行列。
		 */
		inline Matrix operator -=(const Matrix& rhs);
 
		/**
		 * スカラー値で除算する。
		 * @param rhs 除数となる右側のスカラー値。
		 * @return 除算した行列。
		 */
		inline Matrix operator /(const FLOAT& rhs) const;
 
		/**
		 * スカラー値で除算して、代入する。
		 * @param rhs 除数となる右側のスカラー値。
		 * @return 除算して、代入した行列。
		 */
		inline Matrix operator /=(const FLOAT& rhs);
 
		/**
		 * 二つの行列が等しいかどうかを判定する。
		 * @param rhs 右側の行列。
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const Matrix& rhs) const;
 
		/**
		 * 行と列を取り除いて、小行列(SUBmatrix)を作る。
		 * @param row 取り除く行。
		 * @param col 取り除く列。
		 * @return 作った行列。
		 */
		inline Matrix<N - 1> sub(const unsigned int& row, const unsigned int& col) const;
 
		/**
		 * 次元を一つ切り縮める(TRUNCate)。
		 * @return 切り縮めた行列。
		 */
		inline Matrix<N - 1> trunc() const;
	};
 
	/**
	 * 行列を出力する。
	 * @param N 行列の次元数。
	 * @param out 出力先となるストリーム。
	 * @param mat 出力する行列。
	 * @return 出力したストリーム。
	 */
	template <unsigned int N> 
	inline ostream& operator <<(ostream& out, const Matrix<N>& mat);
 
	template <> 
	constexpr Matrix<2> Matrix<2>::IDENTITY() {
		return Matrix<2>{
			1.0, 0.0, 
			0.0, 1.0, 
		};
	}
 
	template <> 
	constexpr Matrix<3> Matrix<3>::IDENTITY() {
		return Matrix<3>{
			1.0, 0.0, 0.0, 
			0.0, 1.0, 0.0, 
			0.0, 0.0, 1.0, 
		};
	}
 
	template <> 
	constexpr Matrix<4> Matrix<4>::IDENTITY() {
		return Matrix<4>{
			1.0, 0.0, 0.0, 0.0, 
			0.0, 1.0, 0.0, 0.0, 
			0.0, 0.0, 1.0, 0.0, 
			0.0, 0.0, 0.0, 1.0, 
		};
	}
};
 
namespace h3d {
	/**
	 * 二つの数が近いかどうかを判定する。
	 * @param x 一つ目の数。
	 * @param y 二つ目の数。
	 * @return 近いならtrue, 遠いならfalse。
	 */
	inline bool nearEqual(const FLOAT& x, const FLOAT& y);
};
 
#include <iostream>
#include <string>
 
namespace h3d {
	using std::cout;
	using std::endl;
	using std::string;
 
	#define ASSERT(pred) assert(#pred, (pred));
 
	#define PRINT(val) cout << #val << "=" << (val) << endl;
 
	/**
	 * アサーションを実行する。
	 * 成功なら標準出力に結果を出力する。
	 * @param pred_str 判定する述語(PREDicate)を記述した文字列(STRing)。
	 * @param pred_res 判定する述語(PREDicate)の結果(RESult)。
	 * trueなら成功，falseなら失敗と判定する。
	 * @throw string 失敗ならメッセージをスローする。
	 */
	inline void assert(const string& pred_str, const bool& pred_res) throw(string);
};
 
namespace h3d {
	class Test { public: virtual void run() const = 0; };
};
 
#include <memory>
#include <vector>
 
namespace h3d {
	using std::shared_ptr;
	using std::vector;
 
	class TestSet {
	public:
		TestSet();
		void run() const;
	protected:
		vector<shared_ptr<Test>> tests;
	};
};
 
namespace h3d {
	class Test1 : public Test {
	public:
		virtual void run() const override;
	};
};
 
namespace h3d {
	class Test2 : public Test {
	public:
		virtual void run() const override;
	};
};
 
namespace h3d {
	class Test3 : public Test {
	public:
		virtual void run() const override;
	};
};
 
////////////////////////////////////////////////////////////////////////////////
 
#include <cmath>
 
namespace h3d {
	using std::tan;
 
	Matrix<3> scaling(const Vector<3>& coef) {
		// 主対角成分に係数を配置する。
		// |a 0 0|
		// |0 b 0|
		// |0 0 c|
		Matrix<3> res(Matrix<3>::IDENTITY());
		for (auto i = 0; i < 3; i++) res.c[i][i] = coef.c[i];
		return res;
	}
 
	Matrix<3> rotation(const Vector<3>& axis, const FLOAT& rad) {
		// 公式に従って行列を作る。
		Vector<3> nmz_axis = axis.normalize();
		FLOAT xx = nmz_axis.c[0] * nmz_axis.c[0];
		FLOAT xy = nmz_axis.c[0] * nmz_axis.c[1];
		FLOAT xz = nmz_axis.c[0] * nmz_axis.c[2];
		FLOAT yy = nmz_axis.c[1] * nmz_axis.c[1];
		FLOAT yz = nmz_axis.c[1] * nmz_axis.c[2];
		FLOAT zz = nmz_axis.c[2] * nmz_axis.c[2];
		FLOAT s = sin(rad);
		FLOAT c = cos(rad);
		FLOAT d = 1.0 - c;
		FLOAT xxd = xx * d;
		FLOAT xyd = xy * d;
		FLOAT zs = nmz_axis.c[2] * s;
		FLOAT xzd = xz * d;
		FLOAT ys = nmz_axis.c[1] * s;
		FLOAT yyd = yy * d;
		FLOAT yzd = yz * d;
		FLOAT xs = nmz_axis.c[0] * s;
		FLOAT zzd = zz * d;
		return Matrix<3>{
			c + xxd,  xyd - zs, xzd + ys, 
			xyd + zs, c + yyd,  yzd - xs, 
			xzd - ys, yzd + xs, c + zzd, 
		};
	}
 
	Matrix<4> translation(const Vector<3>& os) {
		// 単位行列の第4列に平行移動量を配置する。。
		// |1 0 0 x|
		// |0 1 0 y|
		// |0 0 1 z|
		// |0 0 0 1|
		Matrix<4> res(Matrix<4>::IDENTITY());
		for (auto i = 0; i < 3; i++) res.c[i][3] = os.c[i];
		return res;
	}
 
	Matrix<4> view(const Vector<3>& eye_pos, const Vector<3>& cent_pt, const Vector<3>& up_dir) {
		// ビュー行列は、目の位置までの平行移動→視線回転→上方向回転の順序で変換を行う。
 
		// 目の位置が原点にくるように平行移動行列を作る。
		Matrix<4> eye_pos_trans(translation(-eye_pos));
 
		// 視線についての回転行列を作る。視線が逆Z軸になるように回転する。
 
		Matrix<3> ec_iz_rot(Matrix<3>::IDENTITY());
		Vector<3> nmz_ec = (cent_pt - eye_pos).normalize();
		// 視線と逆Z軸の法線を計算する。
		Vector<3> ec_iz_nml = nmz_ec.outpro(-Z_AXIS);
		// 視線と逆Z軸が平行ではないなら、法線まわりに視線と逆Z軸がなす角度で回転する。
		if (ec_iz_nml != Vector<3>::ZERO()) 
			ec_iz_rot = rotation(ec_iz_nml, acos(nmz_ec.inpro(-Z_AXIS)));
		// 視線がZ軸なら、Y軸まわりに反転する。
		else if (nmz_ec == Z_AXIS) ec_iz_rot = rotation(Y_AXIS, M_PI);
 
		// 上方向についての回転行列を作る。上方向がY軸になるように回転する。
 
		Matrix<3> u_y_rot(Matrix<3>::IDENTITY());
		Vector<3> nmz_up_dir = up_dir.normalize();
		// 上方向とY軸の法線を計算する。
		Vector<3> u_y_nml = nmz_up_dir.outpro(Y_AXIS);
		// 上方向とY軸が平行ではないなら、法線まわりに上方向とY軸がなす角度で回転する。
		if (u_y_nml != Vector<3>::ZERO()) 
			u_y_rot = rotation(u_y_nml, acos(nmz_up_dir.inpro(Y_AXIS)));
		// 上方向が逆Y軸なら、Z軸まわりに反転する。
		else if (nmz_up_dir == -Y_AXIS) u_y_rot = rotation(Z_AXIS, M_PI);
 
		// 作った3つの行列を連結する。右から順に配置する。
		return (u_y_rot * ec_iz_rot).ext() * eye_pos_trans;
	}
 
	Matrix<4> persProj(const double& hori_fov, const double& asp_rat, const double& near_dist, const double& far_dist) {
		// 焦点距離(FOCal DISTance)。
		double foc_dist = 1.0 / tan(hori_fov / 2.0);
		// 右平面と近平面が交差するX座標。
		double right_x = near_dist / foc_dist;
		// 上平面と近平面が交差するX座標。
		double top_y = asp_rat * near_dist / foc_dist;
 
		// 公式に従って行列を作る。
		double n2 = near_dist + near_dist;
		double n2f = n2 * far_dist;
		double fsn = far_dist - near_dist;
		double fan = far_dist + near_dist;
		double wid = right_x + right_x;
		double hei = top_y + top_y;
		return Matrix<4>{
			n2 / wid, 0.0,      0.0,        0.0, 
			0.0,      n2 / hei, 0.0,        0.0, 
			0.0,      0.0,      -fan / fsn, -n2f / fsn, 
			0.0,      0.0,      -1.0,       0.0, 
		};
	}
};
 
#include <cmath>
 
namespace h3d {
	inline FLOAT ang2rad(const FLOAT& ang) {
		return M_PI * ang / 180.0;
	}
};
 
#include <iostream>
 
namespace h3d {
	template <unsigned int N> 
	inline void Matrix<N>::clear() {
		// 各成分をゼロクリアする。
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) 
				this->c[row][col] = 0.0;
		}
	}
 
	template <unsigned int N> 
	inline FLOAT Matrix<N>::det() const {
		FLOAT res = 0.0;
		// 1行目の各成分ごとに反復する。
		for (auto col = 0; col < N; col++) {
			// 成分に、それと対応する小行列式を乗算する。
			FLOAT cofac = c[0][col] * sub(0, col).det();
			// (-1)^col
			if ((col & 0x1) == 1) cofac = -cofac;
			// 結果に余因子を加算する。
			res += cofac;
		}
		return res;
	}
 
	template <> 
	inline FLOAT Matrix<2>::det() const {
		// 再帰的な行列式計算の終着点。
		// 2×2の行列式を計算する。
		return this->c[0][0] * this->c[1][1] - this->c[0][1] * this->c[1][0];
	}
 
	template <unsigned int N> 
	inline Matrix<N + 1> Matrix<N>::ext(const Matrix<N + 1>& mat) const {
		Matrix<N + 1> res;
		// 結果の各成分ごとに反復する。
		for (auto row = 0; row < N + 1; row++) {
			for (auto col = 0; col < N + 1; col++) {
				// 行と列がN以内ならこの行列, それ以外なら引数の行列から成分をコピーする。
				if (row < N && col < N) res.c[row][col] = this->c[row][col];
				else res.c[row][col] = mat.c[row][col];
			}
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::inv() const throw(string) {
		// まず行列式を計算し、正則であることを確かめる。
		FLOAT det_res = det();
		if (nearEqual(det_res, 0.0)) throw "行列が正則ではない。";
		// 行列式の逆数を計算しておく。
		FLOAT inv_det = 1.0 / det_res;
		Matrix<N> res;
		// 各成分ごとに反復する。
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) {
				// 行列式の逆数に、対角の小行列式を乗算する。
				res.c[row][col] = inv_det * sub(col, row).det();
				// (-1)^(row + col)
				if (((row + col) & 0x1) == 1) 
					res.c[row][col] = -res.c[row][col];
			}
		}
		return res;
	}
 
	template <unsigned int N> 
	inline bool Matrix<N>::operator !=(const Matrix& rhs) const {
		return !operator ==(rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator *(const Matrix& rhs) const {
		Matrix<N> res;
		// 結果の各成分ごとに反復する。
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) {
				res.c[row][col] = 0.0;
				// 左側は対応する行の各成分, 右側は対応する列の各成分ごとに反復して、乗算する。
				for (auto i = 0; i < N; i++) 
					res.c[row][col] += this->c[row][i] * rhs.c[i][col];
			}
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Vector<N> Matrix<N>::operator *(const Vector<N>& rhs) const {
		Vector<N> res;
		// ベクトルの各成分ごとに反復する。
		for (auto i = 0; i < N; i++) {
			res.c[i] = 0.0;
			// 左側の行列は、対応する行の各成分ごとに反復して、乗算する。
			for (auto j = 0; j < N; j++) 
				res.c[i] += this->c[i][j] * rhs.c[j];
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator *(const FLOAT& rhs) const {
		Matrix<N> res;
		// 各成分にスカラー値を乗算する。
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) 
				res.c[row][col] = this->c[row][col] * rhs;
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator *=(const Matrix& rhs) {
		return *this = operator *(rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator *=(const FLOAT& rhs) {
		return *this = operator *(rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator +(const Matrix& rhs) const {
		Matrix<N> res;
		// 対応する成分同士で加算する。
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) 
				res.c[row][col] = this->c[row][col] + rhs.c[row][col];
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator +=(const Matrix& rhs) {
		return *this = operator +(rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator -() const {
		Matrix<N> res;
		// 各成分の符号を反転する。
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) 
				res.c[row][col] = -this->c[row][col];
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator -(const Matrix& rhs) const {
		return operator +(-rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator -=(const Matrix& rhs) {
		return *this = operator -(rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator /(const FLOAT& rhs) const {
		return operator *(1.0 / rhs);
	}
 
	template <unsigned int N> 
	inline Matrix<N> Matrix<N>::operator /=(const FLOAT& rhs) {
		return *this = operator /(rhs);
	}
 
	template <unsigned int N> 
	bool Matrix<N>::operator ==(const Matrix& rhs) const {
		bool res = true;
		if (&rhs != this) {
			// 各成分ごとに反復する。
			for (auto row = 0; row < N; row++) {
				for (auto col = 0; col < N; col++) {
					// 対応する成分同士が近ければ等しいと判定する。
					if (!nearEqual(this->c[row][col], rhs.c[row][col])) {
						res = false;
						break;
					}
				}
				if (!res) break;
			}
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Matrix<N - 1> Matrix<N>::sub(const unsigned int& row, const unsigned int& col) const {
		Matrix<N - 1> sub;
		auto sub_row = 0;
		// この行列の各成分ごとに反復する。
		// 取り除く行と列については処理をスキップする。
		for (auto sup_row = 0; sup_row < N; sup_row++) {
			if (sup_row == row) continue;
			auto sub_col = 0;
			for (auto sup_col = 0; sup_col < N; sup_col++) {
				if (sup_col == col) continue;
				// 対応する成分をコピーする。
				sub.c[sub_row][sub_col] = this->c[sup_row][sup_col];
				sub_col++;
			}
			sub_row++;
		}
		return sub;
	}
 
	template <unsigned int N> 
	inline Matrix<N - 1> Matrix<N>::trunc() const {
		Matrix<N - 1> res;
		// 結果の各成分ごとに反復する。
		for (auto row = 0; row < N - 1; row++) {
			// 対応する成分をコピーする。
			for (auto col = 0; col < N - 1; col++) 
				res.c[row][col] = this->c[row][col];
		}
		return res;
	}
 
	template <unsigned int N> 
	inline ostream& operator <<(ostream& out, const Matrix<N>& mat) {
		out << "{";
		for (auto row = 0; row < N; row++) {
			out << "{";
			for (auto col = 0; col < N; col++) out << mat.c[row][col] << ", ";
			out << "}, ";
		}
		out << "}";
	}
};
 
#include <algorithm>
#include <cmath>
#include <iostream>
 
namespace h3d {
	using std::copy;
	using std::ostream;
	using std::sqrt;
 
	template <unsigned int N> 
	inline void Vector<N>::clear() {
		// 各成分をゼロクリアする。
		for (auto i = 0; i < N; i++) this->c[i] = 0.0;
	}
 
	template <unsigned int N> 
	inline Vector<N + 1> Vector<N>::ext(const FLOAT& c) const {
		Vector<N + 1> res;
		// このベクトルからN個の成分をコピーして、最後は引数の成分をコピーする。
		copy(this->c, this->c + N, res.c);
		res.c[N] = c;
		return res;
	}
 
	template <unsigned int N> 
	inline FLOAT Vector<N>::inpro(const Vector& rhs) const {
		FLOAT res = 0;
		// 対応する成分同士の積を求めて、合計する。
		for (auto i = 0; i < N; i++) res += this->c[i] * rhs.c[i];
		return res;
	}
 
	template <unsigned int N> 
	inline double Vector<N>::norm() const {
		return sqrt(this->inpro(*this));
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::normalize() const {
		return *this / norm();
	}
 
	template <unsigned int N> 
	inline bool Vector<N>::operator !=(const Vector& rhs) const {
		return !operator ==(rhs);
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::operator *(const FLOAT& rhs) const {
		Vector<N> res;
		// 各成分にスカラー値を乗算する。
		for (auto i = 0; i < N; i++) res.c[i] = this->c[i] * rhs;
		return res;
	}
 
	template <unsigned int N> 
	inline Vector<N>& Vector<N>::operator *=(const FLOAT& rhs) {
		return *this = operator *(rhs);
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::operator +(const Vector& rhs) const {
		Vector<N> res;
		// 対応する成分同士で加算する。
		for (auto i = 0; i < N; i++) res.c[i] = this->c[i] + rhs.c[i];
		return res;
	}
 
	template <unsigned int N> 
	inline Vector<N>& Vector<N>::operator +=(const Vector& rhs) {
		return *this = operator +(rhs);
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::operator -() const {
		Vector<N> res;
		// 各成分の符号を反転する。
		for (auto i = 0; i < N; i++) res.c[i] = -this->c[i];
		return res;
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::operator -(const Vector& rhs) const {
		return operator +(-rhs);
	}
 
	template <unsigned int N> 
	inline Vector<N>& Vector<N>::operator -=(const Vector& rhs) {
		return *this = operator -(rhs);
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::operator /(const FLOAT& rhs) const {
		return operator *(1.0 / rhs);
	}
 
	template <unsigned int N> 
	inline Vector<N>& Vector<N>::operator /=(const FLOAT& rhs) {
		return *this = operator /(rhs);
	}
 
	template <unsigned int N> 
	inline bool Vector<N>::operator ==(const Vector& rhs) const {
		bool res = true;
		if (&rhs != this) {
			// 各成分ごとに反復して、対応する成分同士が近ければ等しいと判定する。
			for (auto i = 0; i < N; i++) {
				if (!nearEqual(this->c[i], rhs.c[i])) {
					res = false;
					break;
				}
			}
		}
		return res;
	}
 
	template <unsigned int N> 
	inline Vector<N> Vector<N>::outpro(const Vector& rhs) const {
		// 各成分について、対応する列を除いた2×2の行列式を計算する。
		// P×Q = < Py・Qz - Pz・Qy, Pz・Qx - Px・Qz, Px・Qy - Py・Qx >
		return Vector<N>{
			this->c[1] * rhs.c[2] - this->c[2] * rhs.c[1], 
			this->c[2] * rhs.c[0] - this->c[0] * rhs.c[2], 
			this->c[0] * rhs.c[1] - this->c[1] * rhs.c[0], 
		};
	}
 
	template <unsigned int N> 
	inline Vector<N - 1> Vector<N>::trunc() const {
		Vector<N - 1> res;
		copy(this->c, this->c + N - 1, res.c);
		return res;
	}
 
	template <unsigned int N> 
	inline ostream& operator <<(ostream& out, const Vector<N>& vec) {
		out << "{";
		for (auto c : vec.c) out << c << ", ";
		out << "}";
	}
};
 
namespace h3d {
	inline bool nearEqual(const FLOAT& x, const FLOAT& y) {
		static const FLOAT THRESHOLD = 0.0000001;
		FLOAT dif = x - y;
		return dif > -THRESHOLD && dif < THRESHOLD;
	}
};
 
#include <iostream>
#include <string>
 
namespace h3d {
	using std::cout;
	using std::endl;
	using std::string;
 
	inline void assert(const string& pred_str, const bool& pred_res) throw(string) {
		if (pred_res) cout << "アサート成功: " << pred_str << endl;
		else throw "アサート失敗: " + pred_str;
	}
};
 
#include <memory>
 
namespace h3d {
	using std::shared_ptr;
 
	TestSet::TestSet() {
		this->tests.push_back(shared_ptr<Test>(new Test1));
		this->tests.push_back(shared_ptr<Test>(new Test2));
		this->tests.push_back(shared_ptr<Test>(new Test3));
	}
 
	void TestSet::run() const { for (auto iter : this->tests) iter->run(); }
};
 
#include <cmath>
 
namespace h3d {
	using std::sqrt;
 
	void Test1::run() const {
		Vector<3> p{2.0, 2.0, 1.0}, q{1.0, -2.0, 0.0};
		ASSERT(nearEqual(p.c[0], 2.0))
		ASSERT(nearEqual(p.c[1], 2.0))
		ASSERT(nearEqual(p.c[2], 1.0))
		ASSERT(nearEqual(q.c[0], 1.0))
		ASSERT(nearEqual(q.c[1], -2.0))
		ASSERT(nearEqual(q.c[2], 0.0))
 
		ASSERT(p.ext() == (Vector<4>{2.0, 2.0, 1.0, 1.0}))
		ASSERT(p.ext(0.0) == (Vector<4>{2.0, 2.0, 1.0, 0.0}))
		ASSERT(p.trunc() == (Vector<2>{2.0, 2.0}))
 
		ASSERT(p == (Vector<3>{2.0, 2.0, 1.0}))
		ASSERT(q == (Vector<3>{1.0, -2.0, 0.0}))
		ASSERT(p != (Vector<3>{1.0, -2.0, 0.0}))
		ASSERT(q != (Vector<3>{2.0, 2.0, 1.0}))
		ASSERT(-p == (Vector<3>{-2.0, -2.0, -1.0}))
		ASSERT(-q == (Vector<3>{-1.0, 2.0, 0.0}))
		ASSERT(p * 2.0 == (Vector<3>{4.0, 4.0, 2.0}))
		ASSERT(q * -3.0 == (Vector<3>{-3.0, 6.0, 0.0}))
		ASSERT(p / 2.0 == (Vector<3>{1.0, 1.0, 0.5}))
		ASSERT(q / -3.0 == (Vector<3>{1.0 / -3.0, -2.0 / -3.0, 0.0}))
 
		ASSERT(p + q == (Vector<3>{3.0, 0.0, 1.0}))
		ASSERT(q + p == (Vector<3>{3.0, 0.0, 1.0}))
		ASSERT(p - q == (Vector<3>{1.0, 4.0, 1.0}))
		ASSERT(q - p == (Vector<3>{-1.0, -4.0, -1.0}))
 
		Vector<3> r;
		r = p;
		ASSERT(r == (Vector<3>{2.0, 2.0, 1.0}))
		ASSERT((r *= 3.0) == (Vector<3>{6.0, 6.0, 3.0}))
		ASSERT((r /= -2.0) == (Vector<3>{-3.0, -3.0, -1.5}))
		ASSERT((r += (Vector<3>{5.0, -1.0, 0.0})) == (Vector<3>{2.0, -4.0, -1.5}))
		ASSERT((r -= (Vector<3>{0.0, -2.0, 4.0})) == (Vector<3>{2.0, -2.0, -5.5}))
		r.clear();
		ASSERT(nearEqual(r.c[0], 0.0))
		ASSERT(nearEqual(r.c[1], 0.0))
		ASSERT(nearEqual(r.c[2], 0.0))
 
		ASSERT(nearEqual(p.inpro(q), -2.0))
		ASSERT(p.outpro(q) == (Vector<3>{2.0, 1.0, -6.0}))
 
		ASSERT(nearEqual(p.norm(), 3.0))
		ASSERT(nearEqual(q.norm(), sqrt(5.0)))
 
		ASSERT(p.normalize() == (Vector<3>{2.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0}))
		ASSERT(q.normalize() == (Vector<3>{1.0 / sqrt(5.0), -2.0 / sqrt(5.0), 0.0}))
	}
};
 
#include <cmath>
 
namespace h3d {
	using std::cos;
	using std::sin;
	using std::sqrt;
 
	void Test2::run() const {
		ASSERT(nearEqual(Matrix<2>::IDENTITY().c[0][0], 1.0))
		ASSERT(nearEqual(Matrix<2>::IDENTITY().c[0][1], 0.0))
		ASSERT(nearEqual(Matrix<2>::IDENTITY().c[1][0], 0.0))
		ASSERT(nearEqual(Matrix<2>::IDENTITY().c[1][1], 1.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[0][0], 1.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[0][1], 0.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[0][2], 0.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[1][0], 0.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[1][1], 1.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[1][2], 0.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[2][0], 0.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[2][1], 0.0))
		ASSERT(nearEqual(Matrix<3>::IDENTITY().c[2][2], 1.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[0][0], 1.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[0][1], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[0][2], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[0][3], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[1][0], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[1][1], 1.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[1][2], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[1][3], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[2][0], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[2][1], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[2][2], 1.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[2][3], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[3][0], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[3][1], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[3][2], 0.0))
		ASSERT(nearEqual(Matrix<4>::IDENTITY().c[3][3], 1.0))
 
		Matrix<2> f{
			1.0,  -2.0, 
			-3.0, 4.0, 
		}, g{
			-9.0, 8.0, 
			7.0,  -6.0, 
		};
 
		ASSERT(nearEqual(f.c[0][0], 1.0))
		ASSERT(nearEqual(f.c[0][1], -2.0))
		ASSERT(nearEqual(f.c[1][0], -3.0))
		ASSERT(nearEqual(f.c[1][1], 4.0))
		ASSERT(nearEqual(g.c[0][0], -9.0))
		ASSERT(nearEqual(g.c[0][1], 8.0))
		ASSERT(nearEqual(g.c[1][0], 7.0))
		ASSERT(nearEqual(g.c[1][1], -6.0))
 
		ASSERT(f.ext() == (Matrix<3>{
			1.0,  -2.0, 0.0, 
			-3.0, 4.0,  0.0, 
			0.0,  0.0,  1.0, 
		}))
		ASSERT(f.ext((Matrix<3>{
			3.0,  -8.0, 5.0, 
			-2.0, 7.0,  4.0, 
			1.0,  -5.0, 0.0, 
		})) == (Matrix<3>{
			1.0,  -2.0, 5.0, 
			-3.0, 4.0,  4.0, 
			1.0,  -5.0, 0.0, 
		}))
		ASSERT(f.trunc() == (Matrix<1>{
			1.0, 
		}))
 
		ASSERT(f == (Matrix<2>{
			1.0,  -2.0, 
			-3.0, 4.0, 
		}))
		ASSERT(g == (Matrix<2>{
			-9.0, 8.0, 
			7.0,  -6.0, 
		}))
		ASSERT(f != (Matrix<2>{
			-9.0, 8.0, 
			7.0,  -6.0, 
		}))
		ASSERT(g != (Matrix<2>{
			1.0,  -2.0, 
			-3.0, 4.0, 
		}))
 
		ASSERT(-f == (Matrix<2>{
			-1.0, 2.0, 
			3.0,  -4.0, 
		}))
		ASSERT(-g == (Matrix<2>{
			9.0,  -8.0, 
			-7.0, 6.0, 
		}))
		ASSERT(f + g == (Matrix<2>{
			-8.0, 6.0, 
			4.0,  -2.0, 
		}))
		ASSERT(f - g == (Matrix<2>{
			10.0,  -10.0, 
			-10.0, 10.0, 
		}))
		ASSERT(f * g == (Matrix<2>{
			-23.0, 20.0, 
			55.0,  -48.0, 
		}))
		ASSERT(g * f == (Matrix<2>{
			-33.0, 50.0, 
			25.0,  -38.0, 
		}))
		ASSERT(f * Matrix<2>::IDENTITY() == (Matrix<2>{
			1.0,  -2.0, 
			-3.0, 4.0, 
		}))
		ASSERT(Matrix<2>::IDENTITY() * g == (Matrix<2>{
			-9.0, 8.0, 
			7.0,  -6.0, 
		}))
		ASSERT(f * 3.0 == (Matrix<2>{
			3.0,  -6.0, 
			-9.0, 12.0, 
		}))
		ASSERT(g / -2.0 == (Matrix<2>{
			4.5,  -4.0, 
			-3.5, 3.0, 
		}))
 
		Matrix<2> h;
		h = f;
		ASSERT(h == (Matrix<2>{
			1.0,  -2.0, 
			-3.0, 4.0, 
		}))
		ASSERT((h += (Matrix<2>{
			0.1,  -0.2, 
			-0.3, 0.4, 
		})) == (Matrix<2>{
			1.1,  -2.2, 
			-3.3, 4.4, 
		}))
		ASSERT((h -= (Matrix<2>{
			-0.9, 0.8, 
			0.7,  -0.6, 
		})) == (Matrix<2>{
			2.0,  -3.0, 
			-4.0, 5.0, 
		}))
		ASSERT((h *= (Matrix<2>{
			-4.0, 3.0, 
			-2.0, 1.0, 
		})) == (Matrix<2>{
			-2.0, 3.0, 
			6.0,  -7.0, 
		}))
		ASSERT((h *= 3.0) == (Matrix<2>{
			-6.0, 9.0, 
			18.0, -21.0, 
		}))
		ASSERT((h /= -2.0) == (Matrix<2>{
			3.0,  -4.5, 
			-9.0, 10.5, 
		}))
		h.clear();
		ASSERT(nearEqual(h.c[0][0], 0.0))
		ASSERT(nearEqual(h.c[0][1], 0.0))
		ASSERT(nearEqual(h.c[1][0], 0.0))
		ASSERT(nearEqual(h.c[1][1], 0.0))
 
		ASSERT(nearEqual((Matrix<2>{
			2.0,  7.0, 
			-3.0, 0.5, 
		}).det(), 22.0))
		ASSERT(nearEqual((Matrix<3>{
			0.0, 0.0, 1.0, 
			0.0, 1.0, 0.0, 
			1.0, 0.0, 0.0, 
		}).det(), -1.0))
		ASSERT(nearEqual((Matrix<3>{
			0.5,              sqrt(3.0) / 2.0, 0.0, 
			-sqrt(3.0) / 2.0, 0.5,             0.0, 
			0.0,              0.0,             1.0, 
		}).det(), 1.0))
		ASSERT(nearEqual((Matrix<3>{
			5,  7,  1, 
			17, 2,  64, 
			10, 14, 2, 
		}).det(), 0.0))
 
		ASSERT((Matrix<3>{
			2.0, 0.0, 0.0, 
			0.0, 3.0, 0.0, 
			0.0, 0.0, 4.0, 
		}).inv() == (Matrix<3>{
			12.0, 0.0, 0.0, 
			0.0,  8.0, 0.0, 
			0.0,  0.0, 6.0, 
		}) / 24.0)
		ASSERT((Matrix<3>{
			2.0, 0.0, 0.0, 
			0.0, 3.0, 0.0, 
			0.0, 0.0, 4.0, 
		}) * ((Matrix<3>{
			12.0, 0.0, 0.0, 
			0.0,  8.0, 0.0, 
			0.0,  0.0, 6.0, 
		}) / 24.0) == Matrix<3>::IDENTITY())
		ASSERT((Matrix<3>{
			1.0, 0.0, 0.0, 
			0.0, 2.0, 2.0, 
			3.0, 0.0, 8.0, 
		}).inv() == (Matrix<3>{
			16.0, 0.0, 0.0, 
			6.0,  8.0, -2.0, 
			-6.0, 0.0, 2.0, 
		}) / 16.0)
		ASSERT((Matrix<3>{
			1.0, 0.0, 0.0, 
			0.0, 2.0, 2.0, 
			3.0, 0.0, 8.0, 
		}) * ((Matrix<3>{
			16.0, 0.0, 0.0, 
			6.0,  8.0, -2.0, 
			-6.0, 0.0, 2.0, 
		}) / 16.0) == Matrix<3>::IDENTITY())
		FLOAT theta_rad = ang2rad(60.0);
		FLOAT cos_res = cos(theta_rad);
		FLOAT sin_res = sin(theta_rad);
		ASSERT((Matrix<3>{
			cos_res, 0.0, -sin_res, 
			0.0,     1.0, 0.0, 
			sin_res, 0.0, cos_res, 
		}).inv() == (Matrix<3>{
			cos_res,  0.0, sin_res, 
			0.0,      1.0, 0.0, 
			-sin_res, 0.0, cos_res, 
		}))
		ASSERT((Matrix<3>{
			cos_res, 0.0, -sin_res, 
			0.0,     1.0, 0.0, 
			sin_res, 0.0, cos_res, 
		}) * (Matrix<3>{
			cos_res,  0.0, sin_res, 
			0.0,      1.0, 0.0, 
			-sin_res, 0.0, cos_res, 
		}) == Matrix<3>::IDENTITY())
		ASSERT((Matrix<4>{
			1.0, 0.0, 0.0, 4.0, 
			0.0, 1.0, 0.0, 3.0, 
			0.0, 0.0, 1.0, 7.0, 
			0.0, 0.0, 0.0, 1.0, 
		}).inv() == (Matrix<4>{
			1.0, 0.0, 0.0, -4.0, 
			0.0, 1.0, 0.0, -3.0, 
			0.0, 0.0, 1.0, -7.0, 
			0.0, 0.0, 0.0, 1.0, 
		}))
		ASSERT((Matrix<4>{
			1.0, 0.0, 0.0, 4.0, 
			0.0, 1.0, 0.0, 3.0, 
			0.0, 0.0, 1.0, 7.0, 
			0.0, 0.0, 0.0, 1.0, 
		}) * (Matrix<4>{
			1.0, 0.0, 0.0, -4.0, 
			0.0, 1.0, 0.0, -3.0, 
			0.0, 0.0, 1.0, -7.0, 
			0.0, 0.0, 0.0, 1.0, 
		}) == Matrix<4>::IDENTITY())
 
		ASSERT((Matrix<2>{
			1.0,  -3.0, 
			-4.0, 6.0, 
		}) * (Vector<2>{2.0, 5.0}) == (Vector<2>{-13.0, 22.0}))
		ASSERT((Matrix<2>{
			0.0,  8.0, 
			-9.0, 3.0, 
		}) * (Vector<2>{1.0, -4.0}) == (Vector<2>{-32.0, -21.0}))
	}
};
 
#include <cmath>
 
namespace h3d {
	using std::sqrt;
 
	void Test3::run() const {
		ASSERT(scaling((Vector<3>{1.0, -2.0, 0.0})) == (Matrix<3>{
			1.0, 0.0,  0.0, 
			0.0, -2.0, 0.0, 
			0.0, 0.0,  0.0, 
		}))
		ASSERT(scaling((Vector<3>{0.0, 9.0, -8.0})) == (Matrix<3>{
			0.0, 0.0, 0.0, 
			0.0, 9.0, 0.0, 
			0.0, 0.0, -8.0, 
		}))
 
		ASSERT(rotation((Vector<3>{1.0, 0.0, 0.0}), ang2rad(30.0)) == (Matrix<3>{
			1.0, 0.0,             0.0, 
			0.0, sqrt(3.0) / 2.0, -0.5, 
			0.0, 0.5,             sqrt(3.0) / 2.0, 
		}))
		ASSERT(rotation((Vector<3>{0.0, 1.0, 0.0}), ang2rad(-45.0)) == (Matrix<3>{
			1.0 / sqrt(2.0), 0.0, -1.0 / sqrt(2.0), 
			0.0,             1.0, 0.0, 
			1.0 / sqrt(2.0), 0.0, 1.0 / sqrt(2.0), 
		}))
		ASSERT(rotation((Vector<3>{0.0, 0.0, 1.0}), ang2rad(60.0)) == (Matrix<3>{
			0.5,             -sqrt(3.0) / 2.0, 0.0, 
			sqrt(3.0) / 2.0, 0.5,              0.0, 
			0.0,             0.0,              1.0, 
		}))
 
		ASSERT(translation((Vector<3>{1.0, -2.0, 0.0})) == (Matrix<4>{
			1.0, 0.0, 0.0, 1.0, 
			0.0, 1.0, 0.0, -2.0, 
			0.0, 0.0, 1.0, 0.0, 
			0.0, 0.0, 0.0, 1.0, 
		}))
		ASSERT(translation((Vector<3>{0.0, 9.0, -8.0})) == (Matrix<4>{
			1.0, 0.0, 0.0, 0.0, 
			0.0, 1.0, 0.0, 9.0, 
			0.0, 0.0, 1.0, -8.0, 
			0.0, 0.0, 0.0, 1.0, 
		}))
 
		ASSERT((view(
			(Vector<3>{0.0, -3.0, 9.0}), (Vector<3>{0.0, -3.0, 0.0})) * 
			(Vector<3>{3.0, 5.0,  7.0}).ext()).trunc() == 
			(Vector<3>{3.0, 8.0,  -2.0}))
		ASSERT((view(
			O_PT, X_AXIS) * 
			(Vector<3>{0.0,  0.0, -5.0}).ext()).trunc() == 
			(Vector<3>{-5.0, 0.0, 0.0}))
		ASSERT((view(
			O_PT, -Y_AXIS) * 
			(Vector<3>{0.0, 0.0, -5.0}).ext()).trunc() == 
			(Vector<3>{0.0, 5.0, 0.0}))
		ASSERT((view(
			O_PT, Z_AXIS) * 
			(Vector<3>{0.0, 0.0, -5.0}).ext()).trunc() == 
			(Vector<3>{0.0, 0.0, 5.0}))
		ASSERT((view(
			O_PT, -Z_AXIS, 
			(Vector<3>{1.0, 0.0, 0.0})) * 
			(Vector<3>{3.0, 0.0, 0.0}).ext()).trunc() == 
			(Vector<3>{0.0, 3.0, 0.0}))
		ASSERT((view(
			O_PT, -Z_AXIS, 
			(Vector<3>{-1.0, 0.0,  0.0})) * 
			(Vector<3>{0.0,  -3.0, 0.0}).ext()).trunc() == 
			(Vector<3>{-3.0, 0.0,  0.0}))
		ASSERT((view(
			O_PT, X_AXIS, 
			(Vector<3>{1.0,  0.0, 0.0})) * 
			(Vector<3>{7.0,  8.0, 0.0}).ext()).trunc() == 
			(Vector<3>{-8.0, 0.0, -7.0}))
 
		ASSERT(persProj(ang2rad(120.0), 0.75, 1.0, 10.0) == (Matrix<4>{
			(1.0 * 2.0) / (sqrt(3.0) * 2.0), 0.0,                                    0.0,                            0.0, 
			0.0,                             (1.0 * 2.0) / (0.75 * sqrt(3.0) * 2.0), 0.0,                            0.0, 
			0.0,                             0.0,                                    -((10.0 + 1.0) / (10.0 - 1.0)), -((1.0 * 10.0 * 2.0) / (10.0 - 1.0)),  
			0.0,                             0.0,                                    -1.0,                           0.0, 
		}))
		ASSERT(persProj(ang2rad(90.0), 0.5, 2.0, 5.0) == (Matrix<4>{
			(2.0 * 2.0) / (2.0 * 2.0), 0.0,                             0.0,                          0.0, 
			0.0,                       (2.0 * 2.0) / (0.5 * 2.0 * 2.0), 0.0,                          0.0, 
			0.0,                       0.0,                             -((5.0 + 2.0) / (5.0 - 2.0)), -((2.0 * 5.0 * 2.0) / (5.0 - 2.0)),  
			0.0,                       0.0,                             -1.0,                         0.0, 
		}))
	}
};
 
#include <iostream>
#include <string>
 
int main() {
	using std::cerr;
	using std::endl;
	using std::string;
 
	try {
		h3d::TestSet().run();
	}
	catch (const string& msg) {
		cerr << msg << endl;
		return 1;
	}
	return 0;
}