Matrix determinant for 4x4

changelog

@@ -1,5 +1,6 @@
 Version 9.0.5
 - Added exists to compound_factory
+- Added sub_matrix, fix and extend determinant calculation
 
 Version 9.0.4
 - Fix various stopping and reconstruction errors

include/cif++/matrix.hpp

@@ -124,6 +124,23 @@ class matrix_expression
 
 		return os;
 	}
+
+	template <typename M2>
+	constexpr bool operator==(const matrix_expression<M2> &m) const
+	{
+		bool same = false;
+		if (dim_m() == m.dim_m() and dim_n() == m.dim_n())
+		{
+			same = true;
+			for (std::size_t i = 0; same and i < m.dim_m(); ++i)
+			{
+				for (std::size_t j = 0; same and j < m.dim_n(); ++j)
+					same = operator()(i, j) == m(i, j);
+			}
+		}
+
+		return same;
+	}
 };
 
 // --------------------------------------------------------------------
@@ -594,6 +611,35 @@ auto operator*(const matrix_expression<M1> &m1, const matrix_expression<M2> &m2)
 
 // --------------------------------------------------------------------
 
+template <typename M2>
+class sub_matrix : public matrix_expression<sub_matrix<M2>>
+{
+  public:
+	sub_matrix(const M2 &m, int i, int j)
+		: m_m(m)
+		, m_i(i)
+		, m_j(j)
+	{
+	}
+
+	constexpr std::size_t dim_m() const { return m_m.dim_m() - 1; } ///< Return dimension m
+	constexpr std::size_t dim_n() const { return m_m.dim_n() - 1; } ///< Return dimension n
+
+	/** Access to the value of element [ @a i, @a j ] */
+	constexpr auto operator()(std::size_t i, std::size_t j) const
+	{
+		return m_m(
+			i >= m_i ? i + 1 : i,
+			j >= m_j ? j + 1 : j);
+	}
+
+  private:
+	const M2 &m_m;
+	std::size_t m_i, m_j;
+};
+
+// --------------------------------------------------------------------
+
 /** Generic routine to calculate the determinant of a matrix
  * 
  * @note This is currently only implemented for fixed matrices of size 3x3
@@ -605,11 +651,23 @@ auto determinant(const M &m);
 template <typename F = float>
 auto determinant(const matrix3x3<F> &m)
 {
-	return (m(0, 0) * (m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1)) +
-			m(0, 1) * (m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2)) +
-			m(0, 2) * (m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0)));
+	return (m(0, 0) * ((m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1))) +
+			m(0, 1) * ((m(1, 2) * m(2, 0) - m(1, 0) * m(2, 2))) +
+			m(0, 2) * ((m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0))));
 }
 
+/** Implementation of the determinant function for fixed size matrices of size 4x4 */
+template <typename F = float>
+F determinant(const matrix4x4<F> &m)
+{
+	return m(0, 0) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 0))) -
+	       m(0, 1) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 1))) +
+	       m(0, 2) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 2))) -
+	       m(0, 3) * determinant(matrix3x3<F>(sub_matrix<decltype(m)>(m, 0, 3)));
+}
+
+// --------------------------------------------------------------------
+
 /** Generic routine to calculate the inverse of a matrix
  * 
  * @note This is currently only implemented for fixed matrices of size 3x3

include/cif++/point.hpp

@@ -30,6 +30,7 @@
 #include <cmath>
 #include <complex>
 #include <cstdint>
+#include <format>
 #include <functional>
 #include <valarray>
 
@@ -365,11 +366,18 @@ class quaternion_type
 	}
 
 	/// \brief test for all zero values
-	constexpr operator bool() const
+	constexpr explicit operator bool() const
 	{
 		return a != 0 or b != 0 or c != 0 or d != 0;
 	}
 
+	/// \brief for debugging e.g.
+	friend std::ostream &operator<<(std::ostream &os, const quaternion_type &rhs)
+	{
+		os << std::format("{{ a: {}, b: {}, c: {}, d: {} }}", rhs.a, rhs.b, rhs.c, rhs.d);
+		return os;
+	}
+
   private:
 	value_type a, b, c, d;
 };

test/CMakeLists.txt

@@ -49,6 +49,7 @@ list(
 	spinner
 	reconstruction
 	validate-pdbx
+	matrix
 )
 
 add_library(test-main OBJECT "${CMAKE_CURRENT_SOURCE_DIR}/test-main.cpp")

test/matrix-test.cpp

@@ -0,0 +1,101 @@
+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2025 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "cif++/matrix.hpp"
+#include "test-main.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include <cif++.hpp>
+
+TEST_CASE("m1")
+{
+	cif::matrix3x3<int> m = cif::identity_matrix<int>(3);
+
+	CHECK(cif::determinant(m) == 1);
+}
+
+TEST_CASE("m2")
+{
+	cif::matrix4x4<int> m = cif::identity_matrix<int>(4);
+
+	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
+	CHECK(ms == cif::identity_matrix<int>(3));
+}
+
+TEST_CASE("m3")
+{
+	cif::matrix4x4<int> m{
+		{ 1, 2, 3, 4,      //
+			5, 6, 7, 8,    //
+			9, 10, 11, 12, //
+			13, 14, 15, 16 }
+	};
+	cif::sub_matrix<cif::matrix4x4<int>> ms(m, 1, 1);
+
+	cif::matrix3x3<int> t{
+		{ 1, 3, 4, 9, 11, 12, 13, 15, 16 }
+	};
+
+	CHECK(ms == t);
+}
+
+TEST_CASE("m4")
+{
+	cif::matrix4x4<int> m{
+		{
+			-2,
+			3,
+			1,
+			0,
+			4,
+			1,
+			-3,
+			2,
+			0,
+			-1,
+			2,
+			5,
+			3,
+			2,
+			0,
+			-4,
+		}
+	};
+
+	std::cout << m << "\n\n";
+
+	std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0)) << "\n\n";
+	std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1)) << "\n\n";
+	std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2)) << "\n\n";
+	std::cout << cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3)) << "\n\n";
+
+	std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 0))) << "\n\n";
+	std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 1))) << "\n\n";
+	std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 2))) << "\n\n";
+	std::cout << cif::determinant(cif::matrix3x3<int>(cif::sub_matrix<decltype(m)>(m, 0, 3))) << "\n\n";
+
+	CHECK(cif::determinant(m) == 332);
+}