行列式をテンプレートメタプログラミングで頑張る

「行列式って何？」という人は線形代数の教科書でも読んでください（ぉ　まあ最近は教科書を買わなくても教科書レベルは大抵webにあるので調べれば簡単に出てきます。

行列式を解く場合は大抵はLU分解して対角要素を掛け算して出すと思います。しかし今回は以下の漸化式(?)を用います。
$det(M_{ij}) = \sum_j (-1)^j \Delta_0j$
ここでΔは余因子行列である。

こちらの方が計算コストが多いのだが以下の理由により採用した

今回使うのが行列要素が整数であり、LU分解だと割り算が出てくるので誤差が出る可能性がある。
LU分解だと対角要素に０があると場合分けが面倒だから
元の行列を汚したくなかったから
そこまででかい行列は扱わない
なんとなく関数型言語っぽく作れそうだから（ぉ

まあそんなわけでとりあえず実装してみた。

template< class T, int D >
struct Determinant{
	T operator()( T** matrix )const{
		//余因子行列の作成
		T* cofactor[D-1];
		for( int d = 0; d < D-1; ++d ){
			cofactor[d] = matrix[d+1] + 1;
		}
		T sum(0);
		for( int d = 0; d < D; ++d ){
			const T matD0 = matrix[d][0];
			if( matD0 != 0 ){
				//余因子行列のdeterminantを求める
				const T det = matD0 * Determinant< T, D-1 >()( cofactor );
				sum = ((d&1)== 0) ? sum + det : sum - det;
			}
			//次の余因子行列の作成
			if( d < D-1 ){
				cofactor[d] = matrix[d] + 1;
			}
		}
		return sum;
	}
};

//２次元の場合の特殊化
template< class T >
struct Determinant< T, 2 >{
	T operator()( T** matrix )const{
		return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
	}
};

漸化式なので微妙にテンプレートメタプログラミングになってます（＾＾；　余因子行列はポインタを微妙にずらしていくことにより実現。これもＣ言語ならではの妙技ですね。
しかしこれ、ループが２個もあり、ｄの値による場合分けがあったりとイマイチ関数型言語っぽくありません。以前どこかのサイトで「ループ使ったら負けかなと思っている」という名言(?)を見たことがありますがまさにそんな感じ。というわけで今回は徹底的にメタメタしてみました。

template< class T, int D >
struct Determinant{
	T operator()( T** matrix )const{
		T* cofactor[D-1];
		return sum< D-1 >( matrix, cofactor );
	}
private:
	typedef Determinant< T, D-1 > DM;

	template< int C > T sum( T** m, T** c )const{
		c[C-1] = m[C] + 1;
		const T s = sum< C-1 >( m, c );
		c[C-1] = m[C-1] + 1;
		//return ( m[C][0] != 0 ) ? add< C&1 >( s, m[C][0] * DM()( c ) ) : s;
		return add< C&1 >( s, m[C][0] * DM()( c ) );
	}
	template<> T sum< 0 >( T** m, T** c )const{
		//return ( m[0][0] != 0 ) ? m[0][0] * DM()( c ) : 0;
		return m[0][0] * DM()( c );
	}

	template< int PLUS > T add( const T& sum, const T& add )const{ return sum + add; }
	template<> T add< 1 >( const T& sum, const T& add )const{ return sum - add; }
};

//２次元の場合の特殊化・・・省略

結果、見事にifとループを消し去りました。コメントアウトの部分は、先頭が０の場合は余因子行列を計算する必要がない(計算しても０を掛けるので消える)のでその分の高速化を入れた場合です。三項演算子が事実上の分岐なので嫌だったのでそこも消してみました。

ちょっとは変態的プログラマに近づけたかな？（ぇ