From 7f83e821ca17ef333e3d10d26c24cc2300987203 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vladim=C3=ADr=20Vondru=C5=A1?= <mosra@centrum.cz>
Date: Thu, 30 Dec 2010 02:07:52 +0100
Subject: [PATCH] Recursive computation of matrix determinant.

---
 src/Math/Matrix.h            | 44 ++++++++++++++++++++++++++++++++++++
 src/Math/Test/MatrixTest.cpp | 14 ++++++++++++
 src/Math/Test/MatrixTest.h   |  1 +
 3 files changed, 59 insertions(+)

diff --git a/src/Math/Matrix.h b/src/Math/Matrix.h
index 2a1a121e5..7b0c596a6 100644
--- a/src/Math/Matrix.h
+++ b/src/Math/Matrix.h
@@ -156,10 +156,54 @@ template<class T, size_t size> class Matrix {
             return out;
         }
 
+        /**
+         * @brief Determinant
+         *
+         * Computed recursively using Laplace's formula expanded down to 2x2
+         * matrix, where the determinant is computed directly. Complexity is
+         * O(n!), the same as when computing the determinant directly.
+         */
+        T determinant() const {
+            T out(0);
+
+            for(size_t col = 0; col != size; ++col)
+                out += ((col & 1) ? -1 : 1)*at(0, col)*ij(0, col).determinant();
+
+            return out;
+        }
+
     private:
         T _data[size*size];
 };
 
+#ifndef DOXYGEN_GENERATING_OUTPUT
+/* Barebone template specialization for 2x2 matrix (only for determinant computation) */
+template<class T> class Matrix<T, 2> {
+    public:
+        inline Matrix(bool identity = true) {
+            memset(_data, 0, 4*sizeof(T));
+            if(identity) {
+                set(0, 0, T(1));
+                set(1, 1, T(1));
+            }
+        }
+        inline Matrix<T, 2>& operator=(const Matrix<T, 2>& other) {
+            if(&other != this) setData(other.data());
+            return *this;
+        }
+        inline void setData(const T* data) { memcpy(_data, data, 4*sizeof(T)); }
+        inline T at(size_t row, size_t col) const { return _data[col*2+row]; }
+        inline void set(size_t row, size_t col, T value) { _data[col*2+row] = value; }
+
+        T determinant() const {
+            return at(0, 0)*at(1, 1) - at(0, 1)*at(1, 0);
+        }
+
+    private:
+        T _data[4];
+};
+#endif
+
 }}
 
 #endif
diff --git a/src/Math/Test/MatrixTest.cpp b/src/Math/Test/MatrixTest.cpp
index bf6463ae2..5522df25d 100644
--- a/src/Math/Test/MatrixTest.cpp
+++ b/src/Math/Test/MatrixTest.cpp
@@ -189,4 +189,18 @@ void MatrixTest::ij() {
     QVERIFY(Matrix4(original).ij(2, 1) == Matrix3(skipped));
 }
 
+void MatrixTest::determinant() {
+    int m[] = {
+        1, 2, 2, 1, 0,
+        2, 3, 2, 1, -2,
+        1, 1, 1, 1, 0,
+        2, 0, 0, 1, 2,
+        3, 1, 0, 1, -2
+    };
+
+    int d = Matrix<int, 5>(m).determinant();
+
+    QVERIFY(d == -2);
+}
+
 }}}
diff --git a/src/Math/Test/MatrixTest.h b/src/Math/Test/MatrixTest.h
index 28bab0ab4..3b190881d 100644
--- a/src/Math/Test/MatrixTest.h
+++ b/src/Math/Test/MatrixTest.h
@@ -32,6 +32,7 @@ class MatrixTest: public QObject {
         void multiplyVector();
         void transposed();
         void ij();
+        void determinant();
 };
 
 }}}