fork download

copy

namespace h3d {
    /** 使用する浮動小数点数型。 */
	using FLOAT = double;
 
	/**
	 * 二つの数が近いかどうかを判定する。
	 * @param lhs 左側(Left Hand Side)の数。
	 * @param rhs 右側(Right Hand Side)の数。
	 * @return 近いならtrue, 遠いならfalse。
	 */
	template<typename X> 
	inline bool nearEqual(const X& lhs, const X& rhs);
};
 
#include <iostream>
 
namespace h3d {
	using std::ostream;
 
	/**
	 * N次元ベクトルクラス。
	 * @param N ベクトルの次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 * この型で実体化したnearEqualを使う。
	 */
	template<unsigned int N, typename C = FLOAT> 
	class Vector {
	public:
		/** 成分(Component)の配列。N個。 */
		C c[N];
 
		/**
		 * 各成分をゼロクリアする。
		 */
		inline void clear();
 
		/**
		 * 次元を一つ拡張(EXTend)する。
		 * @param ext_comp 追加する成分(COMPonent)。
		 * @return 拡張したベクトル。
		 */
		inline Vector<N + 1, C> ext(const C& ext_comp = 1.0) const;
 
		/**
		 * 内積(INner PROduct)を計算する。
		 * @param rhs 右側(Right Hand Side)のベクトル。
		 * @return 計算した内積。
		 */
		inline C inpro(const Vector& rhs) const;
 
		/**
		 * 長さを取得する。
		 * @return 取得した長さ。
		 */
		inline double norm() const;
 
		/**
		 * 正規化し、長さを1にする。
		 * @return 正規化したベクトル。
		 */
		inline Vector normalize() const;
 
		/**
		 * 二つのベクトルが等しくないかどうかを判定する。
		 * @param rhs 右側(Right Hand Side)のベクトル。
		 * @return 等しくないならtrue, 等しいならfalse。
		 */
		inline bool operator !=(const Vector& rhs) const;
 
		/**
		 * スカラー値で乗算する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算したベクトル。
		 */
		inline Vector operator *(const C& rhs) const;
 
		/**
		 * スカラー値で乗算して、代入する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算して、代入したベクトル。
		 */
		inline Vector& operator *=(const C& rhs);
 
		/**
		 * ベクトルを加算する。
		 * @param rhs 加数となる右側(Right Hand Side)のベクトル。
		 * @return 加算したベクトル。
		 */
		inline Vector operator +(const Vector& rhs) const;
 
		/**
		 * ベクトルを加算して、代入する。
		 * @param rhs 加数となる右側(Right Hand Side)のベクトル。
		 * @return 加算して、代入したベクトル。
		 */
		inline Vector& operator +=(const Vector& rhs);
 
		/**
		 * 各成分の符号を反転する。
		 * @return 符号を反転したベクトル。
		 */
		inline Vector operator -() const;
 
		/**
		 * ベクトルを減算する。
		 * @param rhs 減数となる右側(Right Hand Side)のベクトル。
		 * @return 減算したベクトル。
		 */
		inline Vector operator -(const Vector& rhs) const;
 
		/**
		 * ベクトルを減算して、代入する。
		 * @param rhs 減数となる右側(Right Hand Side)のベクトル。
		 * @return 減算して、代入したベクトル。
		 */
		inline Vector& operator -=(const Vector& rhs);
 
		/**
		 * スカラー値で除算する。
		 * @param rhs 除数となる右側(Right Hand Side)のスカラー値。
		 * @return 除算したベクトル。
		 */
		inline Vector operator /(const C& rhs) const;
 
		/**
		 * スカラー値で除算して、代入する。
		 * @param rhs 除数となる右側(Right Hand Side)のスカラー値。
		 * @return 除算して、代入したベクトル。
		 */
		inline Vector& operator /=(const C& rhs);
 
		/**
		 * 二つのベクトルが等しいかどうかを判定する。
		 * @param rhs 右側(Right Hand Side)のベクトル
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const Vector& rhs) const;
 
		/**
		 * 3次元ベクトルの外積(OUTer PROduct)を計算する。
		 * @param rhs 右側(Right Hand Side)のベクトル。
		 * @return 計算した外積。
		 */
		inline Vector outpro(const Vector& rhs) const;
 
		/**
		 * 次元を一つ切り縮める(TRUNCate)。
		 * @return 切り縮めたベクトル。
		 */
		inline Vector<N - 1, C> trunc() const;
 
		/**
		 * ゼロベクトルを取得する。
		 * @return 取得したゼロベクトル。
		 */
		static constexpr Vector ZERO();
	};
 
	/**
	 * N次元ベクトルを出力する。
	 * @param N ベクトルの次元数。
	 * @param C 成分(Component)の型。
	 * @param out 出力先となるストリーム。
	 * @param vec 出力するベクトル。
	 * @return 出力したストリーム。
	 */
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Vector<N, C>& vec);
 
	/** 浮動小数点数の2次元ゼロベクトル。 */
	template<> 
	constexpr Vector<2> Vector<2>::ZERO() {
		return Vector<2>{0.0, 0.0};
	}
 
	/** 浮動小数点数の3次元ゼロベクトル。 */
	template<> 
	constexpr Vector<3> Vector<3>::ZERO() {
		return Vector<3>{0.0, 0.0, 0.0};
	}
 
	/** 浮動小数点数の4次元ゼロベクトル。 */
	template<> 
	constexpr Vector<4> Vector<4>::ZERO() {
		return Vector<4>{0.0, 0.0, 0.0, 0.0};
	}
};
 
#include <iostream>
#include <string>
#include <tuple>
 
namespace h3d {
	using std::ostream;
	using std::string;
	using std::tuple;
 
	/**
	 * N次多項式クラス。
	 * 0次(定数項)からN次までのN + 1個の係数を格納する。
	 * 例) 3x^2 - 2x + 4 → c[0] = 4, c[1] = -2, c[2] = 3
	 * @param N 多項式の次数。項の数より一つ小さい。
	 * @param C 係数(Coefficient)の型。省略ならFLOAT。
	 * この型で実体化したnearEqualを使う。
	 */
	template<unsigned int N, typename C = FLOAT> 
	class PolyExp {
	public:
		/** 各項の係数(Coefficient)を格納する N + 1 次元ベクトル。成分のインデックスが次数に対応する。 */
		Vector<N + 1, C> c;
 
		/**
		 * 多項式を作る。
		 */
		inline PolyExp() = default;
 
		/**
		 * 係数ベクトルから多項式を作る。
		 * @param coef_vec 係数ベクトル(COEFficient VECtor)。
		 * 定数項以外の係数はゼロクリアする。
		 */
		inline explicit PolyExp(const Vector<N + 1, C>& coef_vec);
 
		/**
		 * 定数項の係数から多項式を作ることで、CからPolyExpへの型変換を行う。
		 * @param const_coef 定数項の係数(COEFficient)。
		 * 定数項以外の係数はゼロクリアする。
		 */
		inline PolyExp(const C& const_coef);
 
		/**
		 * 二つの多項式が等しくないかどうかを判定する。
		 * @param rhs 右側(Right Hand Side)の多項式。
		 * @return 等しくないならtrue, 等しいならfalse。
		 */
		inline bool operator !=(const PolyExp& rhs) const;
 
		/**
		 * スカラー値で乗算する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算した多項式。
		 */
		inline PolyExp operator *(const C& rhs) const;
 
		/**
		 * 多項式で乗算する。
		 * @param rhs 乗数となる右側(Right Hand Side)の多項式。
		 * @return 乗算した多項式。
		 */
		inline PolyExp operator *(const PolyExp& rhs) const throw(string);
 
		/**
		 * スカラー値で乗算して、代入する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算して、代入した多項式。
		 */
		inline PolyExp& operator *=(const C& rhs);
 
		/**
		 * 多項式で乗算して、代入する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算して、代入した多項式。
		 */
		inline PolyExp& operator *=(const PolyExp& rhs) throw(string);
 
		/**
		 * 多項式を加算する。
		 * @param rhs 加数となる右側(Right Hand Side)の多項式。
		 * @return 加算した多項式。
		 */
		inline PolyExp operator +(const PolyExp& rhs) const;
 
		/**
		 * 多項式を加算して、代入する。
		 * @param rhs 加数となる右側(Right Hand Side)の多項式。
		 * @return 加算して、代入した多項式。
		 */
		inline PolyExp& operator +=(const PolyExp& rhs);
 
		/**
		 * 各係数の符号を反転する。
		 * @return 符号を反転した多項式。
		 */
		inline PolyExp operator -() const;
 
		/**
		 * 多項式を減算する。
		 * @param rhs 減数となる右側(Right Hand Side)の多項式。
		 * @return 減算した多項式。
		 */
		inline PolyExp operator -(const PolyExp& rhs) const;
 
		/**
		 * 多項式を減算して、代入する。
		 * @param rhs 減数となる右側(Right Hand Side)の多項式。
		 * @return 減算して、代入した多項式。
		 */
		inline PolyExp& operator -=(const PolyExp& rhs);
 
		/**
		 * スカラー値で除算する。
		 * @param rhs 除数となる右側(Right Hand Side)のスカラー値。
		 * @return 除算した多項式。
		 */
		inline PolyExp operator /(const C& rhs) const;
 
		/**
		 * スカラー値で除算して、代入する。
		 * @param rhs 除数となる右側(Right Hand Side)のスカラー値。
		 * @return 除算して、代入した多項式。
		 */
		inline PolyExp& operator /=(const C& rhs);
 
		/**
		 * 二つの多項式が等しいかどうかを判定する。
		 * @param rhs 右側(Right Hand Side)の多項式
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const PolyExp& rhs) const;
 
		/**
		 * 実数根(RooTS)を計算する。
		 * @return 計算した実数根。昇順にソートされている。
		 * 一つ目が実数根の個数。二つ目が実数根を格納したベクトル。
		 */
		inline tuple<unsigned int, Vector<N, C>> rts() const;
	};
 
	/**
	 * 多項式の実数根を計算するクラス。
	 * @param N 多項式の次数。
	 * @param C 係数(Coefficient)の型。
	 */
	template<unsigned int N, typename C> 
	class PolyExpRoots {
	public:
		/** インスタンスを作ることはできない。 */
		PolyExpRoots() = delete;
 
		/**
		 * 多項式の実数根を計算する。
		 * @param poly_exp 計算する多項式(POLYnominal EXPression)。
		 * @return 計算した実数根。昇順にソートされている。
		 * 一つ目が実数根の個数。二つ目が実数根を格納したベクトル。
		 */
		inline static tuple<unsigned int, Vector<N, C>> calculate(const PolyExp<N, C>& poly_exp);
	};
 
	/**
	 * 2次多項式の実数根を計算するクラス。
	 * @param C 係数(Coefficient)の型。
	 */
	template<typename C> 
	class PolyExpRoots<2, C> {
	public:
		/** インスタンスを作ることはできない。 */
		PolyExpRoots() = delete;
 
		/**
		 * 2次多項式の実数根を計算する。
		 * @param poly_exp 計算する2次多項式(POLYnominal EXPression)。
		 * @return 計算した実数根。昇順にソートされている。
		 * 一つ目が実数根の個数。二つ目が実数根を格納したベクトル。
		 */
		inline static tuple<unsigned int, Vector<2, C>> calculate(const PolyExp<2, C>& poly_exp);
	};
 
	/**
	 * 3次多項式の実数根を計算するクラス。
	 * @param C 係数(Coefficient)の型。
	 */
	template<typename C> 
	class PolyExpRoots<3, C> {
	public:
		/** インスタンスを作ることはできない。 */
		PolyExpRoots() = delete;
 
		/**
		 * 3次多項式の実数根を計算する。
		 * @param poly_exp 計算する3次多項式(POLYnominal EXPression)。
		 * @return 計算した実数根。昇順にソートされている。
		 * 一つ目が実数根の個数。二つ目が実数根を格納したベクトル。
		 */
		inline static tuple<unsigned int, Vector<3, C>> calculate(const PolyExp<3, C>& poly_exp);
	};
 
	/**
	 * N次多項式を出力する。
	 * @param N 多項式の次元数。
	 * @param out 出力先となるストリーム。
	 * @param poly_exp 出力する多項式。
	 * @return 出力したストリーム。
	 */
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const PolyExp<N, C>& poly_exp);
};
 
namespace h3d {
	/**
	 * polygonToTrianglesのラムダの説明用の関数型。
	 * @param tri_vert1_ind 分割した三角形の一つ目の頂点インデックス。
	 * @param tri_vert2_ind 分割した三角形の二つ目の頂点インデックス。
	 * @param tri_vert3_ind 分割した三角形の三つ目の頂点インデックス。
	 */
	using POLY2TRI_LAMBDA = void (*)(const unsigned int& tri_vert1_ind, const unsigned int& tri_vert2_ind, const unsigned int& tri_vert3_ind);
 
	/** 原点。 */
	constexpr Vector<3> O_POINT{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};
 
	/**
	 * 角度を弧度に変換する。
	 * @param ang 変換する角度(ANGle)。
	 * @return 変換した弧度。
	 */
	inline FLOAT angleToRadian(const FLOAT& ang);
 
	/**
	 * 多角形を三角形に分割する。
	 * @param poly_begin 多角形の始端の頂点インデックスを指す反復子。
	 * @param poly_end 多角形の終端の頂点インデックスを指す反復子。
	 * @param lambda 分割した三角形ごとに呼び出す{@link h3d::POLY2TRI_LAMBDA ラムダ}。
	 */
	template<typename PI, typename L> 
	inline void polygonToTriangles(PI poly_begin, PI poly_end, L lambda);
};
 
#include <iostream>
#include <string>
#include <tuple>
 
namespace h3d {
	using std::ostream;
	using std::string;
	using std::tuple;
 
	/**
	 * N次元正方行列クラス。
	 * @param N 行列の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 * この型で実体化したnearEqualを使う。
	 */
	template<unsigned int N, typename C = FLOAT> 
	class Matrix {
	public:
		/** 成分(Component)の配列。N×N個。 */
		C c[N][N];
 
		/**
		 * 各成分をゼロクリアする。
		 */
		inline void clear();
 
		/**
		 * 行列式(DETerminant)を計算する。
		 * @return 計算した行列式。
		 */
		inline C det() const;
 
		/**
		 * 固有値(EIGenVALueS)を計算する。
		 * @return 計算した固有値を格納したベクトル。昇順にソートされている。
		 * @throw string 行列が非正則ならメッセージをスローする。
		 */
		inline Vector<N, C> eigvals() const throw(string);
 
		/**
		 * 固有ベクトル(EIGenVECtorS)を計算する。
		 * @return 計算した固有ベクトルを格納したベクトル。
		 * 各固有ベクトルは正規化されている。
		 * @throw string 行列が非正則ならメッセージをスローする。
		 */
		inline Vector<N, Vector<N, C>> eigvecs() const throw(string);
 
		/**
		 * 次元を一つ拡張(EXTend)する。
		 * @param ext_mat この行列(MATrix)から最後の行と列の成分をコピーする。
		 * @return 拡張した行列。
		 */
		inline Matrix<N + 1, C> ext(const Matrix<N + 1, C>& ext_mat = Matrix<N + 1, C>::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 乗数となる右側(Right Hand Side)の行列。
		 * @return 乗算した行列。
		 */
		inline Matrix operator *(const Matrix& rhs) const;
 
		/**
		 * ベクトルで乗算する。
		 * @param rhs 乗数となる右側(Right Hand Side)のベクトル。
		 * @return 計算したベクトル。
		 */
		inline Vector<N> operator *(const Vector<N>& rhs) const;
 
		/**
		 * スカラー値で乗算する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算した行列。
		 */
		inline Matrix operator *(const C& rhs) const;
 
		/**
		 * 行列で乗算して、代入する。
		 * @param rhs 乗数となる右側(Right Hand Side)の行列。
		 * @return 乗算して、代入した行列。
		 */
		inline Matrix operator *=(const Matrix& rhs);
 
		/**
		 * スカラー値で乗算して、代入する。
		 * @param rhs 乗数となる右側(Right Hand Side)のスカラー値。
		 * @return 乗算して、代入した行列。
		 */
		inline Matrix operator *=(const C& rhs);
 
		/**
		 * 行列を加算する。
		 * @param rhs 加数となる右側(Right Hand Side)の行列。
		 * @return 加算した行列。
		 */
		inline Matrix operator +(const Matrix& rhs) const;
 
		/**
		 * 行列を加算して、代入する。
		 * @param rhs 加数となる右側(Right Hand Side)の行列。
		 * @return 加算して、代入した行列。
		 */
		inline Matrix operator +=(const Matrix& rhs);
 
		/**
		 * すべての成分の符号を反転する。
		 * @return 符号を反転した行列。
		 */
		inline Matrix operator -() const;
 
		/**
		 * 行列を減算する。
		 * @param rhs 減数となる右側(Right Hand Side)の行列。
		 * @return 減算した行列。
		 */
		inline Matrix operator -(const Matrix& rhs) const;
 
		/**
		 * 行列を減算して、代入する。
		 * @param rhs 減数となる右側(Right Hand Side)の行列。
		 * @return 減算して、代入した行列。
		 */
		inline Matrix operator -=(const Matrix& rhs);
 
		/**
		 * スカラー値で除算する。
		 * @param rhs 除数となる右側(Right Hand Side)のスカラー値。
		 * @return 除算した行列。
		 */
		inline Matrix operator /(const C& rhs) const;
 
		/**
		 * スカラー値で除算して、代入する。
		 * @param rhs 除数となる右側(Right Hand Side)のスカラー値。
		 * @return 除算して、代入した行列。
		 */
		inline Matrix operator /=(const C& rhs);
 
		/**
		 * 二つの行列が等しいかどうかを判定する。
		 * @param rhs 右側(Right Hand Side)の行列。
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const Matrix& rhs) const;
 
		/**
		 * 被約形にする。連立1次方程式の解を計算するときに使う。
	     * @param con_col 定数ベクトル。省略ならゼロベクトル。
		 * @return 被約系にした行列と定数ベクトル。一つ目は行列, 二つ目は定数ベクトル。
		 */
		inline tuple<Matrix, Vector<N, C>> reduce(const Vector<N, C>& con_col = Vector<N, C>::ZERO()) const;
 
		/**
		 * 行と列を取り除いて、小行列(SUBmatrix)を作る。
		 * @param row 取り除く行。
		 * @param col 取り除く列。
		 * @return 作った行列。
		 */
		inline Matrix<N - 1, C> sub(const unsigned int& row, const unsigned int& col) const;
 
		/**
		 * 次元を一つ切り縮める(TRUNCate)。
		 * @return 切り縮めた行列。
		 */
		inline Matrix<N - 1, C> trunc() const;
	};
 
	/**
	 * N次元正方行列式計算クラス。
	 * @param N 行列の次元数。
	 * @param C 成分(Component)の型。
	 */
	template<unsigned int N, typename C> 
	class MatrixDeterminant {
	public:
		/** インスタンスを作ることはできない。 */
		MatrixDeterminant() = delete;
 
		/**
		 * 行列式を計算する。
		 * @param mat 計算する行列。
		 * @return 計算した行列式。
		 */
		inline static C calculate(const Matrix<N, C>& mat);
	};
 
	/**
	 * 2次元正方行列式計算クラス。
	 * @param C 成分(Component)の型。
	 */
	template<typename C> 
	class MatrixDeterminant<2, C> {
	public:
		/** インスタンスを作ることはできない。 */
		MatrixDeterminant() = delete;
 
		/**
		 * 2次元正方行列式を計算する。
		 * @param mat 計算する2次元行列。
		 * @return 計算した行列式。
		 */
		inline static C calculate(const Matrix<2, C>& mat);
	};
 
	/**
	 * N次元正方行列を出力する。
	 * @param N 行列の次元数。
	 * @param C 成分(Component)の型。
	 * @param out 出力先となるストリーム。
	 * @param mat 出力する行列。
	 * @return 出力したストリーム。
	 */
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Matrix<N, C>& mat);
 
	/** 浮動小数点数の2次元単位行列。 */
	template<> 
	constexpr Matrix<2> Matrix<2>::IDENTITY() {
		return Matrix<2>{
			1.0, 0.0, 
			0.0, 1.0, 
		};
	}
 
	/** 浮動小数点数の3次元単位行列。 */
	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, 
		};
	}
 
	/** 浮動小数点数の4次元単位行列。 */
	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 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);
 
	/**
	 * 透視射影行列を作る。
	 * 透視射影行列はカメラ空間からクリップ空間への変換を行う。
	 * @param hori_fov 水平視野角(HORIzontal Field Of View)。弧度。
	 * @param asp_rat アスペクト比(ASPect RATio)。
	 * @param near_dist 近平面までの距離(DISTance)。
	 * @param far_dist 遠平面までの距離(DISTance)。
	 */
	Matrix<4> perspectiveProjection(const double& hori_fov, const double& asp_rat, const double& near_dist, const double& far_dist);
};
 
#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);
};
 
#include <string>
 
namespace h3d {
	using std::string;
 
	class Test { public: virtual void run() const throw(string) = 0; };
};
 
#include <memory>
#include <string>
#include <vector>
 
namespace h3d {
	using std::shared_ptr;
	using std::string;
	using std::vector;
 
	class TestSet {
	public:
		TestSet();
		void run() const throw(string);
	protected:
		vector<shared_ptr<Test>> tests;
	};
};
 
#include <string>
 
namespace h3d {
	using std::string;
 
	class Test2 : public Test {
	public:
		virtual void run() const throw(string) override;
	};
};
 
#include <string>
 
namespace h3d {
	using std::string;
 
	class Test5 : public Test {
	public:
		virtual void run() const throw(string) override;
	};
};
 
 
////////////////////////////////////////////////////////////////////////////////
 
#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;
	}
};
 
namespace h3d {
	template<typename X> 
	inline bool nearEqual(const X& lhs, const X& rhs) {
		return lhs == rhs;
	}
 
	template<> 
	inline bool nearEqual<FLOAT>(const FLOAT& lhs, const FLOAT& rhs) {
		static const FLOAT THRESHOLD = 0.00000000000001;
		FLOAT dif = lhs - rhs;
		return dif > -THRESHOLD && dif < THRESHOLD;
	}
};
 
#include <cmath>
 
namespace h3d {
	inline FLOAT angleToRadian(const FLOAT& ang) { return M_PI * ang / 180.0; }
 
	template<typename PI, typename L> 
	inline void polygonToTriangles(PI poly_begin, PI poly_end, L lambda) {
		// アルゴリズムを示す。
		// 1. 始端から始める。
		// 2. 現在位置から二つ先が終端以上なら終了する。
		// 3. 始端なら、現在位置のものを一つ目の頂点インデックスとして、5に飛ぶ。
		// 4. 始端ではないなら、一つ前の位置のものを一つ目の頂点インデックスとする。
		// 5. 一つ次の位置に進んで、現在位置のものを二つ目の頂点インデックスとする。
		// 6. 一つ次の位置に進んで、現在位置のものを三つ目の頂点インデックスとする。
		// 7. ラムダを呼ぶ。
		// 8. 一つ前の位置に戻って、2に飛ぶ。
		for (auto iter = poly_begin; iter + 2 < poly_end;) {
			unsigned int tri_vert1_ind, tri_vert2_ind, tri_vert3_ind;
			if (iter == poly_begin) tri_vert1_ind = *iter;
			else tri_vert1_ind = *(iter - 1);
			tri_vert2_ind = *++iter;
			tri_vert3_ind = *(++iter)--;
			lambda(tri_vert1_ind, tri_vert2_ind, tri_vert3_ind);
		}
	}
};
 
#include <algorithm>
#include <iostream>
#include <tuple>
#include <utility>
 
namespace h3d {
	using std::fill;
	using std::get;
	using std::make_tuple;
	using std::swap;
	using std::tuple;
 
	template<unsigned int N, typename C> 
	inline void Matrix<N, C>::clear() {
		// 各成分をゼロクリアする。
		fill((C*)this->c, (C*)this->c + N * N, 0.0);
	}
 
	template<unsigned int N, typename C> 
	inline C Matrix<N, C>::det() const {
		return MatrixDeterminant<N, C>::calculate(*this);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Matrix<N, C>::eigvals() const throw(string) {
		// 特性多項式を作る。
		Matrix<N, PolyExp<N, C>> ch_mat;
		for (auto row = 0; row < N; row++) {
			for (auto col = 0; col < N; col++) {
				PolyExp<N, C> comp(this->c[row][col]);
				if (row == col) comp.c.c[1] = -1.0;
				ch_mat.c[row][col] = comp;
			}
		}
		PolyExp<N, C> ch_poly_exp(ch_mat.det());
		// 特性多項式の実数根を計算する。
		tuple<unsigned int, Vector<N, C>> rts = ch_poly_exp.rts();
		// 行列が非正則なら、N個の実数根を持たない(?)。
		if (get<0>(rts) < N) throw string("行列が非正則。");
		return get<1>(rts);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, Vector<N, C>> Matrix<N, C>::eigvecs() const throw(string) {
		Vector<N, Vector<N, C>> res;
		// 固有値を計算する。
		Vector<N, C> eigvals_res(eigvals());
		// 固有値ごとに反復する。
		for (auto i = 0; i < N; i++) {
			// 被約係数行列を作る。
			Matrix rdc_coef_mat(*this);
			for (auto j = 0; j < N; j++) rdc_coef_mat.c[j][j] -= eigvals_res.c[i];
			rdc_coef_mat = get<0>(rdc_coef_mat.reduce());
 
			// (M - λI)V = 0 を解く。
 
			// 主対角成分ごとに反復する。
			for (auto j = 0; j < N; j++) {
				// 先導成分を含むならスキップする。
				if (!nearEqual(rdc_coef_mat.c[j][j], 0.0)) continue;
				// 固有ベクトルの成分ごとに反復する。
				for (auto k = 0; k < N; k++) {
					// 主対角に対応する成分は1、それ以外は対応する行の成分の符号を反転したもの。
					if (k == j) res.c[i].c[k] = 1.0;
					else res.c[i].c[k] = -rdc_coef_mat.c[k][j];
				}
				break;
			}
			// 正規化する。
			res.c[i] = res.c[i].normalize();
		}
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N + 1, C> Matrix<N, C>::ext(const Matrix<N + 1, C>& ext_mat) const {
		Matrix<N + 1, C> 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] = ext_mat.c[row][col];
			}
		}
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::inv() const throw(string) {
		// まず行列式を計算し、正則であることを確かめる。
		C det_res = det();
		if (nearEqual(det_res, 0.0)) throw string("行列が非正則。");
		// 行列式の逆数を計算しておく。
		C inv_det = 1.0 / det_res;
		Matrix<N, C> 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, typename C> 
	inline bool Matrix<N, C>::operator !=(const Matrix& rhs) const {
		return !operator ==(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator *(const Matrix& rhs) const {
		Matrix<N, C> 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, typename C> 
	inline Vector<N> Matrix<N, C>::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, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator *(const C& rhs) const {
		Matrix<N, C> 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, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator *=(const Matrix& rhs) {
		return *this = operator *(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator *=(const C& rhs) {
		return *this = operator *(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator +(const Matrix& rhs) const {
		Matrix<N, C> 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, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator +=(const Matrix& rhs) {
		return *this = operator +(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator -() const {
		Matrix<N, C> 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, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator -(const Matrix& rhs) const {
		return operator +(-rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator -=(const Matrix& rhs) {
		return *this = operator -(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator /(const C& rhs) const {
		return operator *(1.0 / rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::operator /=(const C& rhs) {
		return *this = operator /(rhs);
	}
 
	template<unsigned int N, typename C> 
	bool Matrix<N, C>::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, typename C> 
	inline tuple<Matrix<N, C>, Vector<N, C>> Matrix<N, C>::reduce(const Vector<N, C>& con_col) const {
		Matrix<N, C> res_mat(*this);
		Vector<N, C> res_con_col(con_col);
		// カレント行インデックスを初期化する。
		auto i = 0;
		// 各列ごとに反復する。
		for (auto j = 0; j < N; j++) {
			// ゼロではない成分を探索する。
			auto k = i;
			for (; k < N; k++) {
				if (!nearEqual(res_mat.c[k][j], 0.0)) break;
			}
			// ゼロではない成分が見つからなければ次の列に移る。
			if (k >= N) continue;
			// 見つかった行がカレント行でなければ、二つの行を入れ替える。
			// 例)  ↓カレント列
			//     | 0 5|5| ←カレント行
			//     |*2 4|6| ←見つかった行
			//        ↓二つの行を入れ替える。
			//     |*2 4|6| ←カレント行
			//     | 0 5|5|
			if (k != i) {
				for (auto s = 0; s < N; s++) 
					swap(res_mat.c[k][s], res_mat.c[i][s]);
				swap(res_con_col.c[k], res_con_col.c[i]);
			}
			// カレント行の先導成分を1にするために、先導成分の逆数を各成分に掛ける。
			// 例)  ↓カレント列
			//     |*2 4|6| ←カレント行
			//     | 2 6|9|
			//        ↓先導成分2の逆数である1/2を各成分に掛ける。
			//     |*1 2|3| ←カレント行
			//     | 2 6|9|
			C inv_ij_comp = 1.0 / res_mat.c[i][j];
			for (auto s = j; s < N; s++) 
				res_mat.c[i][s] *= inv_ij_comp;
			res_con_col.c[i] *= inv_ij_comp;
 
			// カレント列の各成分が、カレント行以外はゼロになるようにする。
			// 例)  ↓カレント列
			//     | 1 2|3| ←カレント行
			//     |*2 6|9| ←ゼロにしたい行
			//        ↓ゼロにしたい行にカレント行×-2を加える。
			//     | 1 2|3| ←カレント行
			//     |*0 2|3| ←ゼロにしたい行
 
			// 各行ごとに反復する。
			for (auto r = 0; r < N; r++) {
				// カレント行ならスキップする。
				if (r == i) continue;
				// 係数を計算しておく。ゼロならスキップ。
				C coef = -res_mat.c[r][j];
				if (nearEqual(coef, 0.0)) continue;
				// 各列ごとに反復して、ゼロにしたい行の成分にカレント行の成分×係数を加える。
				for (auto s = j; s < N; s++) 
					res_mat.c[r][s] += res_mat.c[i][s] * coef;
				res_con_col.c[r] += res_con_col.c[i] * coef;
			}
			i++;
		}
		return make_tuple(res_mat, res_con_col);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N - 1, C> Matrix<N, C>::sub(const unsigned int& row, const unsigned int& col) const {
		Matrix<N - 1, C> 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, typename C> 
	inline Matrix<N - 1, C> Matrix<N, C>::trunc() const {
		Matrix<N - 1, C> 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, typename C> 
	inline C MatrixDeterminant<N, C>::calculate(const Matrix<N, C>& mat) {
		C res = 0.0;
		// 1行目の各成分ごとに反復する。
		for (auto col = 0; col < N; col++) {
			// 成分に、それと対応する小行列式を乗算する。
			C cofac = mat.c[0][col] * mat.sub(0, col).det();
			// (-1)^col
			if ((col & 0x1) == 1) cofac = -cofac;
			// 結果に余因子を加算する。
			res += cofac;
		}
		return res;
	}
 
	template<typename C> 
	inline C MatrixDeterminant<2, C>::calculate(const Matrix<2, C>& mat) {
		// 再帰的な行列式計算の終着点。
		// 2×2の行列式を計算する。
		return mat.c[0][0] * mat.c[1][1] - mat.c[0][1] * mat.c[1][0];
	}
 
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Matrix<N, C>& 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>
#include <tuple>
 
namespace h3d {
	using std::acos;
	using std::cbrt;
	using std::cos;
	using std::make_tuple;
	using std::ostream;
	using std::pow;
	using std::sort;
	using std::tuple;
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>::PolyExp(const Vector<N + 1, C>& coef_vec) : c(coef_vec) {}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>::PolyExp(const C& const_coef) {
		this->c.clear();
		this->c.c[0] = const_coef;
	}
 
	template<unsigned int N, typename C> 
	inline bool PolyExp<N, C>::operator !=(const PolyExp& rhs) const {
		return !operator ==(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C> PolyExp<N, C>::operator *(const C& rhs) const {
		return PolyExp(this->c * rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C> PolyExp<N, C>::operator *(const PolyExp& rhs) const throw(string) {
		PolyExp<N, C> res;
		res.c.clear();
		// 左側の各係数と右側の各係数の組み合わせごとに反復して、乗算する。
		for (auto i = 0; i < N + 1; i++) {
			if (!nearEqual(this->c.c[i], 0.0)) {
				for (auto j = 0; j < N + 1; j++) {
					if (!nearEqual(rhs.c.c[j], 0.0)) {
						if (i + j >= N + 1) throw string("多項式の次数が足りない。");
						res.c.c[i + j] += this->c.c[i] * rhs.c.c[j];
					}
				}
			}
		}
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>& PolyExp<N, C>::operator *=(const C& rhs) {
		return *this = operator *(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>& PolyExp<N, C>::operator *=(const PolyExp& rhs) throw(string) {
		return *this = operator *(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C> PolyExp<N, C>::operator +(const PolyExp& rhs) const {
		return PolyExp(this->c + rhs.c);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>& PolyExp<N, C>::operator +=(const PolyExp& rhs) {
		return *this = operator +(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C> PolyExp<N, C>::operator -() const {
		return PolyExp(-this->c);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C> PolyExp<N, C>::operator -(const PolyExp& rhs) const {
		return operator +(-rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>& PolyExp<N, C>::operator -=(const PolyExp& rhs) {
		return *this = operator -(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C> PolyExp<N, C>::operator /(const C& rhs) const {
		return operator *(1.0 / rhs);
	}
 
	template<unsigned int N, typename C> 
	inline PolyExp<N, C>& PolyExp<N, C>::operator /=(const C& rhs) {
		return *this = operator /(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline bool PolyExp<N, C>::operator ==(const PolyExp& rhs) const {
		return this->c == rhs.c;
	}
 
	template<unsigned int N, typename C> 
	inline tuple<unsigned int, Vector<N, C>> PolyExp<N, C>::rts() const {
		return PolyExpRoots<N, C>::calculate(*this);
	}
 
	template<typename C> 
	inline tuple<unsigned int, Vector<2, C>> PolyExpRoots<2, C>::calculate(const PolyExp<2, C>& poly_exp) {
		unsigned int res_rts_cnt = 0;
		Vector<2, C> res_rts;
		// 各係数をわかりやすい名前の変数に代入する。
		// ax^2 + bx + c
		C a = poly_exp.c.c[2], b = poly_exp.c.c[1], c = poly_exp.c.c[0];
		// 判別式を計算する。
		// d = b^2 - 4ac
		C d = b * b - 4.0 * a * c;
		// 1 / 2a
		C inv_a2 = 1.0 / (a + a);
		// d = 0 なら重根。正なら二つの実数根を持つ。
		if (nearEqual(d, 0.0)) {
			res_rts_cnt = 1;
			res_rts.c[0] = -b * inv_a2;
		}
		else if (d > 0.0) {
			res_rts_cnt = 2;
			res_rts.c[0] = (-b + sqrt(d)) * inv_a2;
			res_rts.c[1] = (-b - sqrt(d)) * inv_a2;
			sort(res_rts.c, res_rts.c + 2);
		}
		return make_tuple(res_rts_cnt, res_rts);
	}
 
	template<typename C> 
	inline tuple<unsigned int, Vector<3, C>> PolyExpRoots<3, C>::calculate(const PolyExp<3, C>& poly_exp) {
		unsigned int res_rts_cnt = 0;
		Vector<3, C> res_rts;
		// 3次の係数が1になるように除算を行う。
		PolyExp<3, C> poly_exp2(poly_exp / poly_exp.c.c[3]);
		// 各係数をわかりやすい名前の変数に代入する。
		// t^3 + at^2 + bt + c
		C a = poly_exp2.c.c[2], b = poly_exp2.c.c[1], c = poly_exp2.c.c[0];
		//          a
		// t = x - ---
		//          3
		// とすると、
		// x^3 + px + q
		// が得られる。ただし、
		//        1
		// p = - ---a^2 + b
		//        3
		//      2         1
		// q = ----a^3 - ---ab + c
		//      27        3
		// である。
		// xの根を求められたら、a/3を引けば、tの根を得られる。
		C a2p = a * a, a3p = a2p * a;
		C p = -a2p / 3.0 + b, q = a3p * 2.0 / 27.0 - a * b / 3.0 + c;
		C pd = p / 3.0, qd = q / 2.0;
		C pd3p = pd * pd * pd, qd2p = qd * qd;
		C dd = pd3p + qd2p;
		if (nearEqual(dd, 0.0) || dd > 0) {
			C r = cbrt(-qd + sqrt(dd));
			if (nearEqual(dd, 0.0)) {
				res_rts_cnt = 2;
				res_rts.c[0] = r + r;
				res_rts.c[1] = -r;
				if (nearEqual(res_rts.c[0], res_rts.c[1])) res_rts_cnt = 1;
				else sort(res_rts.c, res_rts.c + 2);
			}
			else {
				C s = cbrt(-qd - sqrt(dd));
				res_rts_cnt = 1;
				res_rts.c[0] = r + s;
			}
		}
		else {
			res_rts_cnt = 3;
			C theta = acos(-qd / sqrt(-pd3p)) / 3.0;
			C u = sqrt(-pd), u2m = u + u;
			res_rts.c[0] = cos(theta) * u2m;
			res_rts.c[1] = cos(theta + M_PI * 2.0 / 3.0) * u2m;
			res_rts.c[2] = cos(theta - M_PI * 2.0 / 3.0) * u2m;
			sort(res_rts.c, res_rts.c + 3);
		}
		C a3d = a / 3.0;
		for (auto i = 0; i < res_rts_cnt; i++) res_rts.c[i] -= a3d;
		return make_tuple(res_rts_cnt, res_rts);
	}
 
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const PolyExp<N, C>& poly_exp) {
		out << "{";
		for (auto c : poly_exp.c) out << c << ", ";
		out << "}";
	}
};
 
#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_trl(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_trl;
	}
 
	Matrix<4> perspectiveProjection(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;
		// 上平面と近平面が交差するY座標。
		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 <algorithm>
#include <cmath>
#include <iostream>
 
namespace h3d {
	using std::copy;
	using std::fill;
	using std::ostream;
	using std::sqrt;
 
	template<unsigned int N, typename C> 
	inline void Vector<N, C>::clear() {
		// 各成分をゼロクリアする。
		fill(this->c, this->c + N, 0.0);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N + 1, C> Vector<N, C>::ext(const C& ext_comp) const {
		Vector<N + 1, C> res;
		// このベクトルからN個の成分をコピーして、最後は引数の成分をコピーする。
		copy(this->c, this->c + N, res.c);
		res.c[N] = ext_comp;
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline C Vector<N, C>::inpro(const Vector& rhs) const {
		C res = 0.0;
		// 対応する成分同士の積を求めて、合計する。
		for (auto i = 0; i < N; i++) res += this->c[i] * rhs.c[i];
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline double Vector<N, C>::norm() const {
		return sqrt(this->inpro(*this));
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Vector<N, C>::normalize() const {
		return *this / norm();
	}
 
	template<unsigned int N, typename C> 
	inline bool Vector<N, C>::operator !=(const Vector& rhs) const {
		return !operator ==(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Vector<N, C>::operator *(const C& rhs) const {
		Vector<N, C> res;
		// 各成分にスカラー値を乗算する。
		for (auto i = 0; i < N; i++) res.c[i] = this->c[i] * rhs;
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C>& Vector<N, C>::operator *=(const C& rhs) {
		return *this = operator *(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Vector<N, C>::operator +(const Vector& rhs) const {
		Vector<N, C> res;
		// 対応する成分同士で加算する。
		for (auto i = 0; i < N; i++) res.c[i] = this->c[i] + rhs.c[i];
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C>& Vector<N, C>::operator +=(const Vector& rhs) {
		return *this = operator +(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Vector<N, C>::operator -() const {
		Vector<N, C> res;
		// 各成分の符号を反転する。
		for (auto i = 0; i < N; i++) res.c[i] = -this->c[i];
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Vector<N, C>::operator -(const Vector& rhs) const {
		return operator +(-rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C>& Vector<N, C>::operator -=(const Vector& rhs) {
		return *this = operator -(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C> Vector<N, C>::operator /(const C& rhs) const {
		return operator *(1.0 / rhs);
	}
 
	template<unsigned int N, typename C> 
	inline Vector<N, C>& Vector<N, C>::operator /=(const C& rhs) {
		return *this = operator /(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline bool Vector<N, C>::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, typename C> 
	inline Vector<N, C> Vector<N, C>::outpro(const Vector& rhs) const {
		// 各成分について、対応する列を除いた2×2の行列式を計算する。
		// P×Q = < Py・Qz - Pz・Qy, Pz・Qx - Px・Qz, Px・Qy - Py・Qx >
		return Vector<N, C>{
			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, typename C> 
	inline Vector<N - 1, C> Vector<N, C>::trunc() const {
		Vector<N - 1, C> res;
		copy(this->c, this->c + N - 1, res.c);
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Vector<N, C>& vec) {
		out << "{";
		for (auto c : vec.c) out << c << ", ";
		out << "}";
	}
};
 
#include <memory>
#include <string>
 
namespace h3d {
	using std::shared_ptr;
	using std::string;
 
	TestSet::TestSet() {
		this->tests.push_back(shared_ptr<Test>(new Test2));
		this->tests.push_back(shared_ptr<Test>(new Test5));
	}
 
	void TestSet::run() const throw(string) { for (auto iter : this->tests) iter->run(); }
};
 
#include <cmath>
#include <tuple>
 
namespace h3d {
	using std::cos;
	using std::get;
	using std::make_tuple;
	using std::sin;
	using std::sqrt;
	using std::tuple;
 
	void Test2::run() const throw(string) {
		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 = angleToRadian(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}))
 
		ASSERT((Matrix<3>{
			3.0, 2.0,  -3.0, 
			4.0, -3.0, 6.0, 
			1.0, 0.0,  -1.0, 
		}).reduce((Vector<3>{
			5.0, 
			1.0, 
			3.0, 
		})) == make_tuple((Matrix<3>{
			1.0, 0.0, 0.0, 
			0.0, 1.0, 0.0, 
			0.0, 0.0, 1.0, 
		}), (Vector<3>{
			13.0 / 10.0, 
			-2.0, 
			-17.0 / 10.0, 
		})))
		ASSERT((Matrix<3>{
			2.0, 1.0, 3.0, 
			0.0, 1.0, -1.0, 
			1.0, 3.0, -1.0, 
		}).reduce() == make_tuple((Matrix<3>{
			1.0, 0.0, 2.0, 
			0.0, 1.0, -1.0, 
			0.0, 0.0, 0.0, 
		}), (Vector<3>{
			0.0, 
			0.0, 
			0.0, 
		})))
		ASSERT((Matrix<3>{
			4.0, 3.0,  2.0, 
			1.0, -1.0, -3.0, 
			2.0, 3.0,  4.0, 
		}).reduce() == make_tuple((Matrix<3>{
			1.0, 0.0, -1.0, 
			0.0, 1.0, 2.0, 
			0.0, 0.0, 0.0, 
		}), (Vector<3>{
			0.0, 
			0.0, 
			0.0, 
		})))
 
 
		ASSERT((Matrix<2>{
			1.0, 1.0, 
			3.0, -1.0, 
		}).eigvals() == (Vector<2>{-2.0, 2.0}))
		ASSERT((Matrix<3>{
			2.0, 1.0, 0.0, 
			1.0, 1.0, 0.0, 
			0.0, 0.0, -1.0, 
		}).eigvals() == (Vector<3>{-1.0, (3.0 - sqrt(5.0)) / 2.0, (3.0 + sqrt(5.0)) / 2.0}))
		ASSERT((Matrix<3>{
			2.0,  0.0, 0.0, 
			5.0,  2.0, 3.0, 
			-4.0, 3.0, 2.0, 
		}).eigvals() == (Vector<3>{-1.0, 2.0, 5.0}))
 
		ASSERT((Matrix<2>{
			1.0, 1.0, 
			3.0, -1.0, 
		}).eigvecs() == (Vector<2, Vector<2>>{
			(Vector<2>{-1.0 / 3.0 / (sqrt(10.0) / 3.0), 1.0 / (sqrt(10.0) / 3.0)}), 
			(Vector<2>{1.0 / sqrt(2.0), 1.0 / sqrt(2.0)}), 
		}))
		ASSERT((Matrix<3>{
			2.0,  0.0, 0.0, 
			5.0,  2.0, 3.0, 
			-4.0, 3.0, 2.0, 
		}).eigvecs() == (Vector<3, Vector<3>>{
			(Vector<3>{0.0 / sqrt(2.0), -1.0 / sqrt(2.0), 1.0 / sqrt(2.0)}), 
			(Vector<3>{-3.0 / 5.0 / sqrt(2.0), -4.0 / 5.0 / sqrt(2.0), 1.0 / sqrt(2.0)}), 
			(Vector<3>{0.0 / sqrt(2.0), 1.0 / sqrt(2.0), 1.0 / sqrt(2.0)}), 
		}))
	}
};
 
#include <cmath>
#include <tuple>
 
namespace h3d {
	using std::get;
	using std::tuple;
 
	void Test5::run() const throw(string) {
		PolyExp<2> a{Vector<3>{1.0, -2.0, 3.0}};
		ASSERT((PolyExp<2>{Vector<3>{3.0, -2.0, 1.0}}).c == 
			(Vector<3>{3.0, -2.0, 1.0}))
		ASSERT(-(PolyExp<2>{Vector<3>{3.0, -2.0, 1.0}}) == 
			(PolyExp<2>{Vector<3>{-3.0, 2.0, -1.0}}))
		ASSERT((PolyExp<2>{Vector<3>{3.0, -2.0, 1.0}}) * 2.0 == 
			(PolyExp<2>{Vector<3>{6.0, -4.0, 2.0}}))
		ASSERT((PolyExp<2>{Vector<3>{3.0, -2.0, 1.0}}) / 2.0 == 
			(PolyExp<2>{Vector<3>{1.5, -1.0, 0.5}}))
		ASSERT((PolyExp<2>{Vector<3>{3.0, -2.0, 1.0}}) + 
			(PolyExp<2>{Vector<3>{-5.0, 1.0, 0.0}}) == 
			(PolyExp<2>{Vector<3>{-2.0, -1.0, 1.0}}))
		ASSERT((PolyExp<2>{Vector<3>{3.0, -2.0, 1.0}}) - 
			(PolyExp<2>{Vector<3>{-5.0, 1.0, 0.0}}) == 
			(PolyExp<2>{Vector<3>{8.0, -3.0, 1.0}}))
		ASSERT((PolyExp<2>{Vector<3>{4.0, 2.0, 0.0}}) * 
			(PolyExp<2>{Vector<3>{-3.0, 1.0, 0.0}}) == 
			(PolyExp<2>{Vector<3>{-12.0, -2.0, 2.0}}))
 
		tuple<unsigned int, Vector<2>> rts_res2;
		rts_res2 = (PolyExp<2>{Vector<3>{2.0, -4.0, 2.0}}).rts();
		ASSERT(get<0>(rts_res2) == 1)
		ASSERT(nearEqual(get<1>(rts_res2).c[0], 1.0))
		rts_res2 = (PolyExp<2>{Vector<3>{-3.0, -6.0, -3.0}}).rts();
		ASSERT(get<0>(rts_res2) == 1)
		ASSERT(nearEqual(get<1>(rts_res2).c[0], -1.0))
		rts_res2 = (PolyExp<2>{Vector<3>{-2.0, 0.0, 2.0}}).rts();
		ASSERT(get<0>(rts_res2) == 2)
		ASSERT(nearEqual(get<1>(rts_res2).c[0], -1.0))
		ASSERT(nearEqual(get<1>(rts_res2).c[1], 1.0))
		rts_res2 = (PolyExp<2>{Vector<3>{18.0, -3.0, -3.0}}).rts();
		ASSERT(get<0>(rts_res2) == 2)
		ASSERT(nearEqual(get<1>(rts_res2).c[0], -3.0))
		ASSERT(nearEqual(get<1>(rts_res2).c[1], 2.0))
 
		tuple<unsigned int, Vector<3>> rts_res3;
		rts_res3 = (PolyExp<3>{Vector<4>{-2.0, 6.0, -6.0, 2.0}}).rts();
		ASSERT(get<0>(rts_res3) == 1)
		ASSERT(nearEqual(get<1>(rts_res3).c[0], 1.0))
		rts_res3 = (PolyExp<3>{Vector<4>{2.0, 3.0, 0.0, -1.0}}).rts();
		ASSERT(get<0>(rts_res3) == 2)
		ASSERT(nearEqual(get<1>(rts_res3).c[0], -1.0))
		ASSERT(nearEqual(get<1>(rts_res3).c[1], 2.0))
		rts_res3 = (PolyExp<3>{Vector<4>{180.0, 51.0, -12.0, -3.0}}).rts();
		ASSERT(get<0>(rts_res3) == 3)
		ASSERT(nearEqual(get<1>(rts_res3).c[0], -5.0))
		ASSERT(nearEqual(get<1>(rts_res3).c[1], -3.0))
		ASSERT(nearEqual(get<1>(rts_res3).c[2], 4.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;
}