C++TMPで逆行列を計算（する人を作る）

matrix.cpp

// リスト構造
template<int N, typename Cdr> struct cons {};
struct nil {};

// リストから特定の要素を除去
template<int N, typename List> struct exclude { typedef nil val; };
template<int N, typename  Cdr> struct exclude<N,cons<N,Cdr>>
{
	typedef typename exclude<N,Cdr>::val val;
};
template<int N, int M, typename Cdr> struct exclude<N,cons<M,Cdr>>
{
	typedef cons<M,typename exclude<N,Cdr>::val> val;
};

// テンプレートでN*N の行列に対する演算処理を作成
enum: int { N = 3 };
typedef double matrix_t[N][N];

// 数値から数値のリストを作成
template<int M> struct numlist    { typedef cons<M,typename numlist<M+1>::val> val; };
template<>      struct numlist<N> { typedef cons<N,nil>                        val; };

// 行列式を計算するプログラムを作成
struct det
{
	template<typename Rows, typename Cols> struct rows
	{
		double operator()(const matrix_t) { return 1; }
	};
	template<int Num, typename Cdr, typename Cols> struct rows<cons<Num,Cdr>,Cols>
	{
		double operator()(const matrix_t src) {
			return cols<1,Num,Cdr,Cols,Cols>()(src);
		}
	};
	
	template<int Sign, int Row, typename Rows, typename Cols, typename Car> struct cols
	{
		double operator()(const matrix_t) { return 0; }
	};
	template<int Sign, int Row, typename Rows, typename Cols, int Col, typename Cdr>
	struct cols<Sign,Row,Rows,Cols,cons<Col,Cdr>>
	{
		double operator()(const matrix_t src) {
			return Sign
					* src[Row-1][Col-1]
					* rows<Rows,typename exclude<Col,Cols>::val>()(src)
					+ cols<Sign * -1,Row,Rows,Cols,Cdr>         ()(src);
		}
	};
	
	double operator()(const matrix_t src)
	{
		return rows<numlist<1>::val,numlist<1>::val>()(src);
	}
};

// 逆行列を計算するプログラムを構築
struct inverse
{
	template<typename Rows, typename Cols> struct rows
	{
		void operator()(const matrix_t src, matrix_t dest) {}
	};
	template<int Num, typename Cdr, typename Cols> struct rows<cons<Num,Cdr>,Cols>
	{
		
		void operator()(const matrix_t src, matrix_t dest) {
			cols<Num,Cols>()(src, dest);
			rows<Cdr,Cols>()(src, dest);
		}
	};

	template<int Row, typename Cols> struct cols
	{
		void operator()(const matrix_t src, matrix_t dest) {}
	};
	template<int Row, int Col, typename Cdr> struct cols<Row,cons<Col,Cdr>>
	{

		void operator()(const matrix_t src, matrix_t dest) {
			dest[Col-1][Row-1] = det::rows<
					typename exclude<Row,typename numlist<1>::val>::val,
					typename exclude<Col,typename numlist<1>::val>::val>()(src)
					/ det()(src)
					* (1 - (Col + Row) % 2 * 2) // <(^_^;)
					;
			cols<Row,Cdr>()(src, dest);
		}
	};

	void operator()(const matrix_t src, matrix_t dest)
	{
		rows<numlist<1>::val,numlist<1>::val>()(src, dest);
	}
};

#include <algorithm>
#include <random>
#include <time.h>
#include <stdio.h>

// 行列の積
void multiply(double (&out)[N][N],
		const double (&lhs)[N][N],
		const double (&rhs)[N][N])
{
	for(int k = 0; k < N; ++k)
	for(int j = 0; j < N; ++j)
	{
		double sum = 0;
		for(int i = 0; i < N; ++i) {
			sum += lhs[j][i] *  rhs[i][k];
		}
		out[j][k] = sum;
	}
}

// 行列をダンプ
void describe(const double (&src)[N][N])
{
	puts("---[");
	for(int k = 0; k < N; ++k) {
		for(int j = 0; j < N; ++j) printf("%+.3f; ", src[j][k]);
		puts("");
	}
	puts("]---");
}

int main(int argc, char **argv)
{
	std::mt19937 gen(time(NULL));
	for(int i = 0; i < 10; ++i) {
		double a[N][N];
		double b[N][N];
		double c[N][N];
		double *j = reinterpret_cast<double*>(a);
		std::generate(j, j + N*N, gen);
		inverse()(a, b);
		multiply(c, a, b);
		describe(c);
	}
	return 0;
}

operator() の部分、ちゃんと最適化でインライン展開されるのかな、、？