math: make cofactor computation and matrix inversion simpler and more generic.

src/lol/math/matrix.h

@@ -126,8 +126,6 @@ private:
 };
 
 static_assert(sizeof(imat2) == 16, "sizeof(imat2) == 16");
-static_assert(sizeof(umat2) == 16, "sizeof(umat2) == 16");
-
 #if LOL_FEATURE_CXX11_UNRESTRICTED_UNIONS
 static_assert(sizeof(f16mat2) == 8, "sizeof(f16mat2) == 8");
 #endif
@@ -249,8 +247,6 @@ private:
 };
 
 static_assert(sizeof(imat3) == 36, "sizeof(imat3) == 36");
-static_assert(sizeof(umat3) == 36, "sizeof(umat3) == 36");
-
 #if LOL_FEATURE_CXX11_UNRESTRICTED_UNIONS
 static_assert(sizeof(f16mat3) == 18, "sizeof(f16mat3) == 18");
 #endif
@@ -419,21 +415,15 @@ private:
 };
 
 static_assert(sizeof(imat4) == 64, "sizeof(imat4) == 64");
-static_assert(sizeof(umat4) == 64, "sizeof(umat4) == 64");
-
 #if LOL_FEATURE_CXX11_UNRESTRICTED_UNIONS
 static_assert(sizeof(f16mat4) == 32, "sizeof(f16mat4) == 32");
 #endif
 static_assert(sizeof(mat4) == 64, "sizeof(mat4) == 64");
 static_assert(sizeof(dmat4) == 128, "sizeof(dmat4) == 128");
 
-template<typename T> T determinant(mat_t<T,2,2> const &);
-template<typename T> T determinant(mat_t<T,3,3> const &);
-template<typename T> T determinant(mat_t<T,4,4> const &);
-
-template<typename T> mat_t<T,2,2> inverse(mat_t<T,2,2> const &);
-template<typename T> mat_t<T,3,3> inverse(mat_t<T,3,3> const &);
-template<typename T> mat_t<T,4,4> inverse(mat_t<T,4,4> const &);
+/*
+ * Transpose any matrix
+ */
 
 template<typename T, int COLS, int ROWS>
 static inline mat_t<T, ROWS, COLS> transpose(mat_t<T, COLS, ROWS> const &m)
@@ -445,6 +435,74 @@ static inline mat_t<T, ROWS, COLS> transpose(mat_t<T, COLS, ROWS> const &m)
     return ret;
 }
 
+/*
+ * Compute a square submatrix, useful for minors and cofactor matrices
+ */
+
+template<typename T, int N>
+mat_t<T, N - 1, N - 1> submatrix(mat_t<T, N, N> const &m, int i, int j)
+{
+    ASSERT(i >= 0); ASSERT(j >= 0); ASSERT(i < N); ASSERT(j < N);
+
+    mat_t<T, N - 1, N - 1> ret;
+    for (int i2 = 0; i2 < N - 1; ++i2)
+        for (int j2 = 0; j2 < N - 1; ++j2)
+            ret[i2][j2] = m[i2 + (i2 >= i)][j2 + (j2 >= j)];
+
+    return ret;
+}
+
+/*
+ * Compute square matrix cofactor
+ */
+
+template<typename T, int N>
+T cofactor(mat_t<T, N, N> const &m, int i, int j)
+{
+    ASSERT(i >= 0); ASSERT(j >= 0); ASSERT(i < N); ASSERT(j < N);
+    T tmp = determinant(submatrix(m, i, j));
+    return ((i + j) & 1) ? -tmp : tmp;
+}
+
+/*
+ * Compute square matrix determinant, with a specialisation for 1×1 matrices
+ */
+
+template<typename T, int N>
+T determinant(mat_t<T, N, N> const &m)
+{
+    T ret = (T)0;
+    for (int i = 0; i < N; ++i)
+        ret += m[i][0] * cofactor(m, i, 0);
+    return ret;
+}
+
+template<typename T>
+T const & determinant(mat_t<T, 1, 1> const &m)
+{
+    return m[0][0];
+}
+
+/*
+ * Compute square matrix inverse
+ */
+
+template<typename T, int N>
+mat_t<T, N, N> inverse(mat_t<T, N, N> const &m)
+{
+    mat_t<T, N, N> ret;
+
+    T d = determinant(m);
+    if (d)
+    {
+        /* Transpose result matrix on the fly */
+        for (int i = 0; i < N; ++i)
+            for (int j = 0; j < N; ++j)
+                ret[i][j] = cofactor(m, j, i) / d;
+    }
+    return ret;
+}
+
 /*
  * Matrix-vector and vector-matrix multiplication
  */

src/math/matrix.cpp

@@ -17,126 +17,6 @@
 namespace lol
 {
 
-/*
- * Return the determinant of a 2×2 matrix.
- */
-static inline float det2(float a, float b,
-                         float c, float d)
-{
-    return a * d - b * c;
-}
-
-/*
- * Return the determinant of a 3×3 matrix.
- */
-static inline float det3(float a, float b, float c,
-                         float d, float e, float f,
-                         float g, float h, float i)
-{
-    return a * (e * i - h * f)
-         + b * (f * g - i * d)
-         + c * (d * h - g * e);
-}
-
-/*
- * Return the cofactor of the (i,j) entry in a 2×2 matrix.
- */
-static inline float cofact(mat2 const &m, int i, int j)
-{
-    float tmp = m[(i + 1) & 1][(j + 1) & 1];
-    return ((i + j) & 1) ? -tmp : tmp;
-}
-
-/*
- * Return the cofactor of the (i,j) entry in a 3×3 matrix.
- */
-static inline float cofact(mat3 const &m, int i, int j)
-{
-    return det2(m[(i + 1) % 3][(j + 1) % 3],
-                m[(i + 2) % 3][(j + 1) % 3],
-                m[(i + 1) % 3][(j + 2) % 3],
-                m[(i + 2) % 3][(j + 2) % 3]);
-}
-
-/*
- * Return the cofactor of the (i,j) entry in a 4×4 matrix.
- */
-static inline float cofact(mat4 const &m, int i, int j)
-{
-    return det3(m[(i + 1) & 3][(j + 1) & 3],
-                m[(i + 2) & 3][(j + 1) & 3],
-                m[(i + 3) & 3][(j + 1) & 3],
-                m[(i + 1) & 3][(j + 2) & 3],
-                m[(i + 2) & 3][(j + 2) & 3],
-                m[(i + 3) & 3][(j + 2) & 3],
-                m[(i + 1) & 3][(j + 3) & 3],
-                m[(i + 2) & 3][(j + 3) & 3],
-                m[(i + 3) & 3][(j + 3) & 3]) * (((i + j) & 1) ? -1.0f : 1.0f);
-}
-
-template<> float determinant(mat2 const &m)
-{
-    return det2(m[0][0], m[0][1],
-                m[1][0], m[1][1]);
-}
-
-template<> mat2 inverse(mat2 const &m)
-{
-    mat2 ret;
-    float d = determinant(m);
-    if (d)
-    {
-        d = 1.0f / d;
-        for (int j = 0; j < 2; j++)
-            for (int i = 0; i < 2; i++)
-                ret[j][i] = cofact(m, i, j) * d;
-    }
-    return ret;
-}
-
-template<> float determinant(mat3 const &m)
-{
-    return det3(m[0][0], m[0][1], m[0][2],
-                m[1][0], m[1][1], m[1][2],
-                m[2][0], m[2][1], m[2][2]);
-}
-
-template<> mat3 inverse(mat3 const &m)
-{
-    mat3 ret;
-    float d = determinant(m);
-    if (d)
-    {
-        d = 1.0f / d;
-        for (int j = 0; j < 3; j++)
-            for (int i = 0; i < 3; i++)
-                ret[j][i] = cofact(m, i, j) * d;
-    }
-    return ret;
-}
-
-template<> float determinant(mat4 const &m)
-{
-    float ret = 0;
-    for (int n = 0; n < 4; n++)
-        ret += m[n][0] * cofact(m, n, 0);
-    return ret;
-}
-
-template<> mat4 inverse(mat4 const &m)
-{
-    mat4 ret;
-    float d = determinant(m);
-    if (d)
-    {
-        d = 1.0f / d;
-        for (int j = 0; j < 4; j++)
-            for (int i = 0; i < 4; i++)
-                ret[j][i] = cofact(m, i, j) * d;
-    }
-    return ret;
-}
-
 template<> mat3 mat3::scale(float x)
 {
     mat3 ret(1.0f);