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;
 
	template<unsigned int N, typename C> class Vector;
 
	/**
	 * N次元ベクトル外積計算クラス。
	 * @param N ベクトルの次元数。
	 * @param C 成分(Component)の型。
	 */
	template<unsigned int N, typename C> 
	class VectorOuterProduct {
	public:
		using TYPE = Vector<N, C>;
 
		/** インスタンスを作ることはできない。 */
		VectorOuterProduct() = delete;
 
		/**
		 * N次元ベクトルの外積(OUTer PROduct)を計算する。
		 * @param lhs 左側(Left Hand Side)のベクトル。
		 * @param rhs 右側(Right Hand Side)のベクトル。
		 * @return 計算した外積。
		 */
		inline static TYPE calculate(const Vector<N, C>& lhs, const Vector<N, C>& rhs);
	};
 
	/**
	 * 2次元ベクトル外積計算クラス。
	 * @param C 成分(Component)の型。
	 */
	template<typename C> 
	class VectorOuterProduct<2, C> {
	public:
		using TYPE = C;
 
		/** インスタンスを作ることはできない。 */
		VectorOuterProduct() = delete;
 
		/**
		 * 2次元ベクトルの外積(OUTer PROduct)を計算する。
		 * @param lhs 左側(Left Hand Side)の2次元ベクトル。
		 * @param rhs 右側(Right Hand Side)の2次元ベクトル。
		 * @return 計算した外積。
		 */
		inline static TYPE calculate(const Vector<2, C>& lhs, const Vector<2, C>& rhs);
	};
 
	/**
	 * 3次元ベクトル外積計算クラス。
	 * @param C 成分(Component)の型。
	 */
	template<typename C> 
	class VectorOuterProduct<3, C> {
	public:
		using TYPE = Vector<3, C>;
 
		/** インスタンスを作ることはできない。 */
		VectorOuterProduct() = delete;
 
		/**
		 * 3次元ベクトルの外積(OUTer PROduct)を計算する。
		 * @param lhs 左側(Left Hand Side)の3次元ベクトル。
		 * @param rhs 右側(Right Hand Side)の3次元ベクトル。
		 * @return 計算した外積。
		 */
		inline static TYPE calculate(const Vector<3, C>& lhs, const Vector<3, C>& rhs);
	};
 
	/**
	 * 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;
 
		/**
		 * 長さの自乗を計算する。
		 * @return 計算した長さの自乗。
		 */
		inline double normsqr() 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 typename VectorOuterProduct<N, C>::TYPE 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次多項式(POLYnominal EXPression)クラス。
	 * 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 乗算した多項式。
		 * @throw string 次数が足りないならメッセージをスローする。
		 */
		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);
};
 
#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:
		/** 行ベクトル(Row)の配列。N個。 */
		Vector<N, C> r[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 Matrix 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 Vector<N> operator *=(const Vector<N>& 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_vec 定数ベクトル(CONstant VECtor)。省略ならゼロベクトル。
		 * @return 被約系にした行列と定数ベクトル。一つ目は行列, 二つ目は定数ベクトル。
		 */
		inline tuple<Matrix, Vector<N, C>> reduce(const Vector<N, C>& con_vec = 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;
 
		/**
		 * 転置(TRansPosed)行列を作る。
		 * @return 作った転置行列。
		 */
		inline Matrix trp() 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次元列ベクトルを右側のN次元正方行列で乗算する。
	 * @param N 行列の次元数。
	 * @param C 成分(Component)の型。
	 * @param lhs 左側のN次元列ベクトル。
	 * @param rhs 右側のN次元正方行列。
	 * @return 乗算したN次元正方行列。
	 */
	template<unsigned int N, typename C> 
	inline Matrix<N, C> operator *(const Vector<N, C>& lhs, const Matrix<N, C>& rhs);
 
	/**
	 * 左側のN次元列ベクトルを右側のN次元行ベクトルで乗算する。
	 * @param N 行列の次元数。
	 * @param C 成分(Component)の型。
	 * @param lhs 左側のN次元列ベクトル。
	 * @param rhs 右側のN次元行ベクトル。
	 * @return 乗算したN次元正方行列。
	 */
	template<unsigned int N, typename C> 
	inline Matrix<N, C> operator *(const Vector<N, C>& lhs, const Vector<N, C>& rhs);
 
	/**
	 * 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>{
			Vector<2>{1.0, 0.0}, 
			Vector<2>{0.0, 1.0}, 
		};
	}
 
	/** 浮動小数点数の3次元単位行列。 */
	template<> 
	constexpr Matrix<3> Matrix<3>::IDENTITY() {
		return Matrix<3>{
			Vector<3>{1.0, 0.0, 0.0}, 
			Vector<3>{0.0, 1.0, 0.0}, 
			Vector<3>{0.0, 0.0, 1.0}, 
		};
	}
 
	/** 浮動小数点数の4次元単位行列。 */
	template<> 
	constexpr Matrix<4> Matrix<4>::IDENTITY() {
		return Matrix<4>{
			Vector<4>{1.0, 0.0, 0.0, 0.0}, 
			Vector<4>{0.0, 1.0, 0.0, 0.0}, 
			Vector<4>{0.0, 0.0, 1.0, 0.0}, 
			Vector<4>{0.0, 0.0, 0.0, 1.0}, 
		};
	}
};
 
#include <iostream>
 
namespace h3d {
	using std::ostream;
 
	/**
	 * ボックス。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 */
	template<unsigned int N, typename C = FLOAT> 
	class Box {
	public:
		/** 平面(Plane)の配列。2×N個。 */
		Vector<N + 1, C> p[2 * N];
 
		/**
		 * 二つのボックスが等しくないかどうかを比較する。
		 * @param rhs 右側(Right Hand Side)のボックス
		 * @return 等しくないならtrue, 等しいならfalse。
		 */
		inline bool operator !=(const Box& rhs) const;
 
		/**
		 * 二つのボックスが等しいかどうかを比較する。
		 * @param rhs 右側(Right Hand Side)のボックス
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const Box& rhs) const;
	};
 
	/**
	 * 球。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 */
	template<unsigned int N, typename C = FLOAT> 
	class Sphere {
	public:
		/** 中心点(Center)。 */
		Vector<N, C> c;
		/** 半径(Radius)。 */
		C r;
 
		/**
		 * 二つの球が等しくないかどうかを比較する。
		 * @param rhs 右側(Right Hand Side)の球
		 * @return 等しくないならtrue, 等しいならfalse。
		 */
		inline bool operator !=(const Sphere& rhs) const;
 
		/**
		 * 二つの球が等しいかどうかを比較する。
		 * @param rhs 右側(Right Hand Side)の球
		 * @return 等しいならtrue, 等しくないならfalse。
		 */
		inline bool operator ==(const Sphere& rhs) const;
	};
 
	/**
	 * 角度を弧度に変換する。
	 * @param C 成分(Component)の型。
	 * @param ang 変換する角度(ANGle)。
	 * @return 変換した弧度。
	 */
	template<typename C> 
	inline C angleToRadian(const C& ang);
 
	/**
	 * 多角形を三角形に分割する。
	 * @param I 多角形の頂点の反復子の型。
	 * @param P 分割した三角形ごとに、頂点を出力するラムダの型。
	 * @param poly_begin 多角形の始端の頂点を指す反復子。
	 * @param poly_end 多角形の終端の頂点を指す反復子。
	 * @param put_tri 分割した三角形ごとに、頂点の組を出力するラムダ。
	 * @param put_tri_vert1 一つ目の頂点。
	 * @param put_tri_vert2 二つ目の頂点。
	 * @param put_tri_vert3 三つ目の頂点。
	 */
	template<typename I, typename P> 
	inline void dividePolygonIntoTriangles(I poly_begin, I poly_end, P put_tri);
 
	/**
	 * 球がボックスと交差しているかどうかを判定する。
	 * @param C 成分(Component)の型。
	 * @param sphr 判定する球(SPHeRe)。
	 * @param box ボックス。
	 * @return 交差しているならtrue, 交差していないならfalse。
	 */
	template<unsigned int N, typename C> 
	inline bool intersecting(const Sphere<N, C>& sphr, const Box<N, C>& box);
 
	/**
	 * 境界ボックスを作る。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 * @param I 頂点の反復子の型。
	 * @param G 頂点を取得するラムダの型。
	 * @param vert_begin 始端の頂点を指す反復子。
	 * @param vert_end 終端の頂点を指す反復子。
	 * @param get_vert 頂点を取得するラムダ。
	 * @param get_vert_iter 頂点を指す反復子。
	 * @param get_vert_return 頂点ベクトル。
	 * @return 作った境界ボックス。
	 */
	template<unsigned int N, typename C = FLOAT, typename I, typename G> 
	inline Box<N, C> makeBoundingBox(I vert_begin, I vert_end, G get_vert);
 
	/**
	 * 境界球を作る。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 * @param I 頂点の反復子の型。
	 * @param G 頂点を取得するラムダの型。
	 * @param vert_begin 始端の頂点を指す反復子。
	 * @param vert_end 終端の頂点を指す反復子。
	 * @param get_vert 頂点を取得するラムダ。
	 * @param get_vert_iter 頂点を指す反復子。
	 * @param get_vert_return 頂点ベクトル。
	 * @return 作った境界球。
	 */
	template<unsigned int N, typename C = FLOAT, typename I, typename G> 
	inline Sphere<N, C> makeBoundingSphere(I vert_begin, I vert_end, G get_vert);
 
	/** 
	 * 原点(Origin POINT)。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 */
	template<unsigned int N, typename C = FLOAT> 
	inline constexpr Vector<N, C> O_POINT();
 
	/**
	 * N次元ボックスを出力する。
	 * @param N ベクトルの次元数。
	 * @param C 成分(Component)の型。
	 * @param out 出力先となるストリーム。
	 * @param box 出力するボックス。
	 * @return 出力したストリーム。
	 */
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Box<N, C>& box);
 
	/**
	 * N次元球を出力する。
	 * @param N ベクトルの次元数。
	 * @param C 成分(Component)の型。
	 * @param out 出力先となるストリーム。
	* @param sphr 出力する球(SPHeRe)。
	 * @return 出力したストリーム。
	 */
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Sphere<N, C>& sphr);
 
	/**
	 * 主成分分析を行い、頂点の集合に対して自然に沿う座標系を計算する。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 * @param I 頂点の反復子の型。
	 * @param G 頂点を取得するラムダの型。
	 * @param vert_begin 始端の頂点を指す反復子。
	 * @param vert_end 終端の頂点を指す反復子。
	 * @param get_vert 頂点を取得するラムダ。
	 * @param get_vert_iter 頂点を指す反復子。
	 * @param get_vert_return 頂点ベクトル。
	 * @return 計算した座標系を格納した行列。
	 * 各座標軸は正規化されて、主軸から順に格納されている。
	 */
	template<unsigned int N, typename C = FLOAT, typename I, typename G> 
	inline Matrix<N, C> principalComponentAnalysis(I vert_begin, I vert_end, G get_vert);
 
	/** 
	 * X軸ベクトル。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 */
	template<unsigned int N, typename C = FLOAT> 
	inline constexpr Vector<N, C> X_AXIS();
 
	/** 
	 * Y軸ベクトル。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 */
	template<unsigned int N, typename C = FLOAT> 
	inline constexpr Vector<N, C> Y_AXIS();
 
	/** 
	 * Z軸ベクトル。
	 * @param N 空間の次元数。
	 * @param C 成分(Component)の型。省略ならFLOAT。
	 */
	template<unsigned int N, typename C = FLOAT> 
	inline constexpr Vector<N, C> Z_AXIS();
 
	template<> 
	constexpr Vector<2> O_POINT<2>() { return Vector<2>{0.0, 0.0}; };
 
	template<> 
	constexpr Vector<3> O_POINT<3>() { return Vector<3>{0.0, 0.0, 0.0}; };
 
	template<> 
	constexpr Vector<2> X_AXIS<2>() { return Vector<2>{1.0, 0.0}; };
 
	template<> 
	constexpr Vector<3> X_AXIS<3>() { return Vector<3>{1.0, 0.0, 0.0}; };
 
	template<> 
	constexpr Vector<2> Y_AXIS<2>() { return Vector<2>{0.0, 1.0}; };
 
	template<> 
	constexpr Vector<3> Y_AXIS<3>() { return Vector<3>{0.0, 1.0, 0.0}; };
 
	template<> 
	constexpr Vector<3> Z_AXIS<3>() { return Vector<3>{0.0, 0.0, 1.0}; };
};
 
#include <iomanip>
#include <iostream>
#include <string>
 
namespace h3d {
	using std::cout;
	using std::dec;
	using std::endl;
	using std::hex;
	using std::setfill;
	using std::setw;
	using std::setprecision;
	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 Test4 : public Test {
	public:
		virtual void run() const throw(string) override;
	};
};
 
////////////////////////////////////////////////////////////////////////////////
 
namespace h3d {
	template<typename X> 
	inline bool nearEqual(const X& lhs, const X& rhs) {
		return lhs == rhs;
	}
 
	template<> 
	inline bool nearEqual<double>(const double& lhs, const double& rhs) {
		static const double THRESHOLD = 0.00000000000001;
		double dif = lhs - rhs;
		return dif > -THRESHOLD && dif < THRESHOLD;
	}
 
	template<> 
	inline bool nearEqual<float>(const float& lhs, const float& rhs) {
		static const float THRESHOLD = 0.000001;
		float dif = lhs - rhs;
		return dif > -THRESHOLD && dif < THRESHOLD;
	}
};
 
#include <algorithm>
#include <cmath>
#include <iostream>
 
namespace h3d {
	using std::copy;
	using std::fill;
	using std::ostream;
	using std::sqrt;
 
	template<typename C> 
	inline typename VectorOuterProduct<2, C>::TYPE VectorOuterProduct<2, C>::calculate(const Vector<2, C>& lhs, const Vector<2, C>& rhs) {
		return lhs.c[0] * rhs.c[1] - lhs.c[1] * rhs.c[0];
	}
 
	template<typename C> 
	inline typename VectorOuterProduct<3, C>::TYPE VectorOuterProduct<3, C>::calculate(const Vector<3, C>& lhs, const Vector<3, C>& rhs) {
		Vector<3, C> res;
		for (auto i = 0; i < 3; ++i) {
			auto j = (i + 1) % 3, k = (j + 1) % 3;
			res.c[i] = (Vector<2, C>{lhs.c[j], lhs.c[k]}).outpro((Vector<2, C>{rhs.c[j], rhs.c[k]}));
		}
		return res;
	}
 
	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(normsqr());
	}
 
	template<unsigned int N, typename C> 
	inline double Vector<N, C>::normsqr() const {
		return 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 typename VectorOuterProduct<N, C>::TYPE Vector<N, C>::outpro(const Vector& rhs) const {
		return VectorOuterProduct<N, C>::calculate(*this, rhs);
	}
 
	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 << "}";
		return out;
	}
};
 
#include <algorithm>
#include <cmath>
#include <functional>
#include <iostream>
#include <tuple>
 
namespace h3d {
	using std::acos;
	using std::cbrt;
	using std::cos;
	using std::get;
	using std::greater;
	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)) continue;
			for (auto j = 0; j < N + 1; ++j) {
				if (nearEqual(rhs.c.c[j], 0.0)) continue;
				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 {
		tuple<unsigned int, Vector<N, C>> res(PolyExpRoots<N, C>::calculate(*this));
		if (get<0>(res) >= 2) sort(get<1>(res).c, get<1>(res).c + get<0>(res), greater<C>());
		return res;
	}
 
	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_a2m = 1.0 / (a + a);
		// 判別式がゼロなら重根。正なら二つの実数根を持つ。
		if (nearEqual(d, 0.0)) {
			res_rts_cnt = 1;
			res_rts.c[0] = -b * inv_a2m;
		}
		else if (d > 0.0) {
			res_rts_cnt = 2;
			res_rts.c[0] = (-b + sqrt(d)) * inv_a2m;
			res_rts.c[1] = (-b - sqrt(d)) * inv_a2m;
		}
		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 sqr_a = a * a, cb_a = sqr_a * a;
		C p = -sqr_a / 3.0 + b, q = cb_a * 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 {
				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;
		}
		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 << "}";
		return out;
	}
};
 
#include <iostream>
#include <tuple>
#include <utility>
 
namespace h3d {
	using std::get;
	using std::make_tuple;
	using std::ostream;
	using std::swap;
	using std::tuple;
 
	template<unsigned int N, typename C> 
	inline void Matrix<N, C>::clear() {
		for (auto i = 0; i < N; ++i) this->r[i].clear();
	}
 
	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->r[row].c[col]);
				if (row == col) comp.c.c[1] = -1.0;
				ch_mat.r[row].c[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 Matrix<N, C> Matrix<N, C>::eigvecs() const throw(string) {
		Matrix 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.r[j].c[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.r[j].c[j], 0.0)) continue;
				// 固有ベクトルの成分ごとに反復する。
				for (auto k = 0; k < N; ++k) {
					// 主対角に対応する成分は-1、それ以外は対応する行の成分。
					if (k == j) res.r[i].c[k] = -1.0;
					else res.r[i].c[k] = rdc_coef_mat.r[k].c[j];
				}
				break;
			}
			// 正規化する。
			res.r[i] = res.r[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.r[row].c[col] = this->r[row].c[col];
				else res.r[row].c[col] = ext_mat.r[row].c[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.r[row].c[col] = inv_det * sub(col, row).det();
				// (-1)^(row + col)
				if (((row + col) & 0x1) == 1) 
					res.r[row].c[col] = -res.r[row].c[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 res;
		// 右側の行列を転置する。
		Matrix trp_rhs(rhs.trp());
		// 各成分ごとに反復する。
		for (auto row = 0; row < N; ++row) {
			for (auto col = 0; col < N; ++col) 
				// 左側は対応する行、右側は対応する列の内積を計算する。
				res.r[row].c[col] = this->r[row].inpro(trp_rhs.r[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] = this->r[i].inpro(rhs);
		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) res.r[row] = this->r[row] * 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 Vector<N> Matrix<N, C>::operator *=(const Vector<N>& 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) res.r[row] = this->r[row] + rhs.r[row];
		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) res.r[row] = -this->r[row];
		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) {
				if (this->r[row] != rhs.r[row]) {
					res = false;
					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_vec) const {
		Matrix<N, C> res_mat(*this);
		Vector<N, C> res_con_vec(con_vec);
		// カレント行インデックスを初期化する。
		auto i = 0;
		// 各列ごとに反復する。
		for (auto j = 0; j < N; ++j) {
			// ゼロではない成分を探索する。
			auto k = i;
			for (; k < N; ++k) {
				if (!nearEqual(res_mat.r[k].c[j], 0.0)) break;
			}
			// ゼロではない成分が見つからなければ次の列に移る。
			if (k >= N) continue;
			// 見つかった行がカレント行でなければ、二つの行を入れ替える。
			// 例)  ↓カレント列
			//     | 0 5|5| ←カレント行
			//     |*2 4|6| ←見つかった行
			//        ↓二つの行を入れ替える。
			//     |*2 4|6| ←カレント行
			//     | 0 5|5|
			if (k != i) {
				swap(res_mat.r[k], res_mat.r[i]);
				swap(res_con_vec.c[k], res_con_vec.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.r[i].c[j];
			res_mat.r[i] *= inv_ij_comp;
			res_con_vec.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.r[r].c[j];
				if (nearEqual(coef, 0.0)) continue;
				// ゼロにしたい行にカレント行×係数を加える。
				res_mat.r[r] += res_mat.r[i] * coef;
				res_con_vec.c[r] += res_con_vec.c[i] * coef;
			}
			++i;
		}
		return make_tuple(res_mat, res_con_vec);
	}
 
	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.r[sub_row].c[sub_col] = this->r[sup_row].c[sup_col];
				++sub_col;
			}
			++sub_row;
		}
		return sub;
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> Matrix<N, C>::trp() const {
		Matrix res;
		// 各成分ごとに反復して、行と列を転置してコピーする。
		for (auto row = 0; row < N; ++row) {
			for (auto col = 0; col < N; ++col) 
				res.r[row].c[col] = this->r[col].c[row];
		}
		return res;
	}
 
	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.r[row].c[col] = this->r[row].c[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.r[0].c[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) {
		// 再帰的な行列式計算の終着点。
		return mat.r[0].outpro(mat.r[1]);
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> operator *(const Vector<N, C>& lhs, const Matrix<N, C>& rhs) {
		Matrix<N, C> res;
		// 右側の行列を転置する。
		Matrix<N, C> trp_rhs(rhs.trp());
		// 各成分ごとに反復して、行に対応する左側の成分を対応する右側の列で乗算する。
		for (auto row = 0; row < N; ++row) {
			for (auto col = 0; col < N; ++col) {
				res.r[row].c[col] = 0.0;
				for (auto i = 0; i < N; i++) res.r[row].c[col] += lhs.c[row] * trp_rhs.r[col].c[i];
			}
		}
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline Matrix<N, C> operator *(const Vector<N, C>& lhs, const Vector<N, C>& rhs) {
		Matrix<N, C> res;
		// 各成分ごとに反復して、行に対応する左側の成分を列に対応する右側の成分で乗算する。
		for (auto row = 0; row < N; ++row) {
			for (auto col = 0; col < N; ++col) res.r[row].c[col] = lhs.c[row] * rhs.c[col];
		}
		return res;
	}
 
	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.r[row].c[col] << ", ";
			out << "}, ";
		}
		out << "}";
		return out;
	}
};
 
#include <algorithm>
#include <cmath>
#include <iostream>
 
namespace h3d {
	using std::max;
	using std::min;
	using std::ostream;
	using std::sqrt;
 
	template<unsigned int N, typename C> 
	inline bool Box<N, C>::operator !=(const Box& rhs) const {
		return !operator ==(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline bool Box<N, C>::operator ==(const Box& rhs) const {
		bool res = true;
		for (auto i = 0; i < 2 * N; i++) {
			if (this->p[i] != rhs.p[i]) {
				res = false;
				break;
			}
		}
		return res;
	}
 
	template<unsigned int N, typename C> 
	inline bool Sphere<N, C>::operator !=(const Sphere& rhs) const {
		return !operator ==(rhs);
	}
 
	template<unsigned int N, typename C> 
	inline bool Sphere<N, C>::operator ==(const Sphere& rhs) const {
		return this->c == rhs.c && nearEqual(this->r, rhs.r);
	}
 
	template<typename C> 
	inline C angleToRadian(const C& ang) { return M_PI * ang / 180.0; }
 
	template<typename I, typename P> 
	inline void dividePolygonIntoTriangles(I poly_begin, I poly_end, P put_tri) {
		// アルゴリズムを示す。
		// 1. 始端から始める。
		// 2. 現在位置から二つ先が終端以上なら終了する。
		// 3. 始端なら、現在位置のものを一つ目の頂点インデックスとして、5に飛ぶ。
		// 4. 始端ではないなら、一つ前の位置のものを一つ目の頂点インデックスとする。
		// 5. 一つ次の位置に進んで、現在位置のものを二つ目の頂点インデックスとする。
		// 6. 一つ次の位置に進んで、現在位置のものを三つ目の頂点インデックスとする。
		// 7. ラムダを呼ぶ。
		// 8. 一つ前の位置に戻って、2に飛ぶ。
		for (auto iter = poly_begin; iter + 2 < poly_end;) {
			I tri_vert1, tri_vert2, tri_vert3;
			if (iter == poly_begin) tri_vert1 = iter;
			else tri_vert1 = iter - 1;
			tri_vert2 = ++iter;
			tri_vert3 = (++iter)--;
			put_tri(tri_vert1, tri_vert2, tri_vert3);
		}
	}
 
	template<unsigned int N, typename C> 
	inline bool intersecting(const Sphere<N, C>& sphr, const Box<N, C>& box) {
		bool res = true;
		for (auto i = 0; i < 2 * N; i++) {
			C d = box.p[i].trunc().inpro(sphr.c) + box.p[i].c[N] + sphr.r;
			if (!nearEqual(d, 0.0) && d < 0.0) {
				res = false;
				break;
			}
		}
		return res;
	}
 
	template<unsigned int N, typename C, typename I, typename G> 
	inline Box<N,C> makeBoundingBox(I vert_begin, I vert_end, G get_vert) {
		Box<N, C> res;
		// 頂点の集合に対して自然に沿う座標系を計算する。
		Matrix<N, C> cs(principalComponentAnalysis<N, C>(vert_begin, vert_end, get_vert));
		// 各座標軸ごとに反復して、平面を作る。
		for (auto i = 0; i < N; ++i) {
			C min_d, max_d;
			for (auto iter = vert_begin; iter != vert_end; ++iter) {
				Vector<N, C> vert(get_vert(iter));
				if (iter == vert_begin) min_d = max_d = vert.inpro(cs.r[i]);
				else {
					min_d = min(min_d, vert.inpro(cs.r[i]));
					max_d = max(max_d, vert.inpro(cs.r[i]));
				}
			}
			res.p[2 * i + 0] = cs.r[i].ext(-min_d);
			res.p[2 * i + 1] = (-cs.r[i]).ext(max_d);
		}
		return res;
	}
 
	template<unsigned int N, typename C, typename I, typename G> 
	inline Sphere<N, C> makeBoundingSphere(I vert_begin, I vert_end, G get_vert) {
		// 頂点の集合に対して自然に沿う座標系を計算する。
		Matrix<N, C> cs(principalComponentAnalysis<N, C>(vert_begin, vert_end, get_vert));
		// 各頂点を反復して、主軸平面との距離が最小と最大の頂点を取得する。
		Vector<N, C> min, max;
		C min_d, max_d;
		for (auto iter = vert_begin; iter != vert_end; ++iter) {
			Vector<N, C> vert(get_vert(iter));
			C d = vert.inpro(cs.r[0]);
			if (iter == vert_begin) {
				min = max = vert;
				min_d = max_d = d;
			}
			else {
				if (d < min_d) {
					min = vert;
					min_d = d;
				}
				else if (d > max_d) {
					max = vert;
					max_d = d;
				}
			}
		}
		// 仮の中心点と半径を計算する。
		Vector<N, C> c((min + max) / 2.0);
		C sqr_r = (min - c).normsqr();
		C r = sqrt(sqr_r);
		// 包含しない頂点があれば、包含するように調整する。
		for (auto iter = vert_begin; iter != vert_end; ++iter) {
			Vector<N, C> vert(get_vert(iter));
			C sqr_rd = (vert - c).normsqr();
			if (sqr_rd > sqr_r) {
				Vector<N, C> g(c - (vert - c) * r / sqrt(sqr_rd));
				c = (g + vert) / 2.0;
				sqr_r = (vert - c).normsqr();
				r = sqrt(sqr_r);
			}
		}
		return Sphere<N, C>{c, r};
	}
 
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Box<N, C>& box) {
		out << "{";
		for (auto i = 0; i < 2 * N; ++i) {
			out << box.p[i];
			out << ", ";
		}
		out << "}";
		return out;
	}
 
	template<unsigned int N, typename C> 
	inline ostream& operator <<(ostream& out, const Sphere<N, C>& sphr) {
		out << "{";
		out << sphr.c << ", " << sphr.r;
		out << "}";
		return out;
	}
 
	template<unsigned int N, typename C, typename I, typename G> 
	inline Matrix<N, C> principalComponentAnalysis(I vert_begin, I vert_end, G get_vert) {
		// 頂点の合計と個数を計算する。
		Vector<N, C> sum;
		sum.clear();
		C n = 0.0;
		for (auto iter = vert_begin; iter != vert_end; ++iter) {
			sum += get_vert(iter);
			n += 1.0;
		}
		// 平均位置を計算する。
		Vector<N, C> avg(sum / n);
		// 共分散行列を作る。
		Matrix<N, C> vcm;
		vcm.clear();
		for (auto iter = vert_begin; iter != vert_end; ++iter) {
			Vector<N, C> psa = get_vert(iter) - avg;
			vcm += psa * psa;
		}
		vcm /= n;
		// 共分散行列の固有ベクトルを座標系として返す。
		return vcm.eigvecs();
	}
};
 
#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>
#include <string>
 
namespace h3d {
	using std::shared_ptr;
	using std::string;
 
	TestSet::TestSet() {
		this->tests.push_back(shared_ptr<Test>(new Test4));
	}
 
	void TestSet::run() const throw(string) { for (auto iter : this->tests) iter->run(); }
};
 
#include <cmath>
#include <vector>
 
namespace h3d {
	using std::sqrt;
	using std::vector;
 
	void Test4::run() const throw(string) {
		Vector<3> pnts3[] = {
			Vector<3>{-1.0, -2.0, 1.0}, 
			Vector<3>{ 1.0,  0.0, 2.0}, 
			Vector<3>{ 2.0, -1.0, 3.0}, 
			Vector<3>{ 2.0, -1.0, 2.0}, 
		};
		auto get_vert = [&](const Vector<3>* iter) { return *iter; };
		ASSERT(principalComponentAnalysis<3>(pnts3, pnts3 + 4, get_vert) == (Matrix<3>{
			Vector<3>{-0.83339666730631745128, -0.33030475324794755787,  -0.44311258716553691972}, 
			Vector<3>{-0.25734203789881582303,  0.94146015975144747845,  -0.21777934505104812324}, 
			Vector<3>{ 0.48910639993033028228, -0.067465084173902054032, -0.86961105786702130871}, 
		}))
		ASSERT(makeBoundingSphere<3>(pnts3, pnts3 + 4, get_vert) == (Sphere<3>{
			Vector<3>{0.5, -1.5, 2.0}, 1.8708286933869706647}))
 
		Box<3> box{{
			Vector<4>{ 1.0,  0.0,  0.0, 1.0}, 
			Vector<4>{-1.0,  0.0,  0.0, 1.0}, 
			Vector<4>{ 0.0,  1.0,  0.0, 2.0}, 
			Vector<4>{ 0.0, -1.0,  0.0, 2.0}, 
			Vector<4>{ 0.0,  0.0,  1.0, 3.0}, 
			Vector<4>{ 0.0,  0.0, -1.0, 3.0}, 
		}};
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 2.1,  0.0,  0.0}, 1.0}), box) == false)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 2.0,  0.0,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 1.9,  0.0,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{-1.9,  0.0,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{-2.0,  0.0,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{-2.1,  0.0,  0.0}, 1.0}), box) == false)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  3.1,  0.0}, 1.0}), box) == false)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  3.0,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  2.9,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0, -2.9,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0, -3.0,  0.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0, -3.1,  0.0}, 1.0}), box) == false)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  0.0,  4.1}, 1.0}), box) == false)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  0.0,  4.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  0.0,  3.9}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  0.0, -3.9}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  0.0, -4.0}, 1.0}), box) == true)
		ASSERT(intersecting((Sphere<3>{Vector<3>{ 0.0,  0.0, -4.1}, 1.0}), box) == false)
	}
};
 
#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;
}