求逆矩阵的c++实现

int inc(int i, int N)
{
if(i >= N)
return i % N;
else if(i < 0)
return i + N;
return i;
}

template <typename T, int N>
T CalcDet(T (&m)

)
{
if(2 == N) return m[0][0] * m[1][1] - m[0][1] * m[1][0];

T det = 0;
T pos, neg;

for(int i = 0; i < N; ++i) {
pos = 1;
for(int j = 0; j < N; ++j) {
pos *= m[j][inc(i+j, N)];
}
det += pos;
}

for(int i = N-1; i >= 0; --i) {
neg = 1;
for(int j = 0; j < N; ++j) {
neg *= m[j][inc(i-j, N)];
}
det -= neg;
}

return det;
}

template <typename T, int N>
void Transpose(T (&m)

)
{
T temp;
for(int i = 0; i < N; ++i) {
for(int j = i + i; j < N; ++j) {
temp = m[i][j];
m[i][j] = m[j][i];
m[j][i] = temp;
}
}
}

template <typename T, int N>
bool Inverse(const T (&src)

, double (&dst)

)
{
if(N <= 2) return false;

T det = CalcDet(src);
T submat[N-1][N-1];

for(int i = 0; i < N; ++i) {
for(int j = 0; j < N; ++j) {
int k, l, s, t;
for(k = 0, s = 0; s < N-1; ++s, ++k) {
if(s == i) ++k;
for(l = 0, t = 0; t < N-1; ++t, ++l) {
if(t == j) ++l;
submat[s][t] = src[k][l];
}
}
double sign = (i + j & 1) == 0 ? 1.0 : -1.0;
dst[i][j] = sign * CalcDet(submat) / det;
}
}
Transpose(dst);
return true;
}