From 7c907f6f3e157e3dee655433dbe5b04b5af38f5b Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Tue, 20 Jan 2015 08:51:13 +0000
Subject: [PATCH] matrix: determinant computing using LU-decomposition

---
 src/lol/math/matrix.h | 36 +++++++++++++++++-------------------
 src/t/math/matrix.cpp | 30 +++++++++++++++---------------
 2 files changed, 32 insertions(+), 34 deletions(-)

diff --git a/src/lol/math/matrix.h b/src/lol/math/matrix.h
index 48afa468..78ae50ff 100644
--- a/src/lol/math/matrix.h
+++ b/src/lol/math/matrix.h
@@ -509,25 +509,6 @@ T cofactor(mat_t<T, 2, 2> const &m, int i, int j)
     return (i ^ j) ? -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 LU-decomposition
  */
@@ -557,6 +538,23 @@ void lu_decomposition(mat_t<T, N, N> const &m, mat_t<T, N, N> & L, mat_t<T, N, N
     }
 }
 
+/*
+ * 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)
+{
+    mat_t<T, N, N> L, U;
+    lu_decomposition(m, L, U);
+
+    T det = 1;
+    for (int i = 0 ; i < N ; ++i)
+        det *= U[i][i];
+
+    return det;
+}
+
 /*
  * Compute square matrix inverse
  */
diff --git a/src/t/math/matrix.cpp b/src/t/math/matrix.cpp
index 49a471b9..3306e7b0 100644
--- a/src/t/math/matrix.cpp
+++ b/src/t/math/matrix.cpp
@@ -51,19 +51,19 @@ lolunit_declare_fixture(MatrixTest)
         float d1, d2;
 
         d1 = determinant(tri2);
-        lolunit_assert_equal(d1, 2.0f);
+        lolunit_assert_doubles_equal(d1, 2.0f, 1e-5);
         d2 = determinant(inv2);
-        lolunit_assert_equal(d2, -1.0f);
+        lolunit_assert_doubles_equal(d2, -1.0f, 1e-5);
 
         d1 = determinant(tri3);
-        lolunit_assert_equal(d1, 6.0f);
+        lolunit_assert_doubles_equal(d1, 6.0f, 1e-5);
         d2 = determinant(inv3);
-        lolunit_assert_equal(d2, 1.0f);
+        lolunit_assert_doubles_equal(d2, 1.0f, 1e-5);
 
         d1 = determinant(tri4);
-        lolunit_assert_equal(d1, 24.0f);
+        lolunit_assert_doubles_equal(d1, 24.0f, 1e-5);
         d2 = determinant(inv4);
-        lolunit_assert_equal(d2, -1.0f);
+        lolunit_assert_doubles_equal(d2, -1.0f, 1e-5);
     }
 
     lolunit_declare_test(Multiplication)
@@ -166,17 +166,17 @@ lolunit_declare_fixture(MatrixTest)
 
         mat3 m2 = m0 * m1;
 
-        lolunit_assert_equal(m2[0][0], 1.0f);
-        lolunit_assert_equal(m2[1][0], 0.0f);
-        lolunit_assert_equal(m2[2][0], 0.0f);
+        lolunit_assert_doubles_equal(m2[0][0], 1.0f, 1e-5);
+        lolunit_assert_doubles_equal(m2[1][0], 0.0f, 1e-5);
+        lolunit_assert_doubles_equal(m2[2][0], 0.0f, 1e-5);
 
-        lolunit_assert_equal(m2[0][1], 0.0f);
-        lolunit_assert_equal(m2[1][1], 1.0f);
-        lolunit_assert_equal(m2[2][1], 0.0f);
+        lolunit_assert_doubles_equal(m2[0][1], 0.0f, 1e-5);
+        lolunit_assert_doubles_equal(m2[1][1], 1.0f, 1e-5);
+        lolunit_assert_doubles_equal(m2[2][1], 0.0f, 1e-5);
 
-        lolunit_assert_equal(m2[0][2], 0.0f);
-        lolunit_assert_equal(m2[1][2], 0.0f);
-        lolunit_assert_equal(m2[2][2], 1.0f);
+        lolunit_assert_doubles_equal(m2[0][2], 0.0f, 1e-5);
+        lolunit_assert_doubles_equal(m2[1][2], 0.0f, 1e-5);
+        lolunit_assert_doubles_equal(m2[2][2], 1.0f, 1e-5);
     }
 
     lolunit_declare_test(Inverse4x4)