Merge pull request #160 from orange-cpp/feaure/gjk-epa-improvement

Feaure/gjk epa improvement

.idea/dictionaries/project.xml

@@ -1,7 +0,0 @@
-<component name="ProjectDictionaryState">
-  <dictionary name="project">
-    <words>
-      <w>vmprotect</w>
-    </words>
-  </dictionary>
-</component>

CMakeLists.txt

@@ -31,10 +31,6 @@ option(OMATH_SUPRESS_SAFETY_CHECKS
 option(OMATH_ENABLE_COVERAGE "Enable coverage" OFF)
 option(OMATH_ENABLE_FORCE_INLINE
        "Will for compiler to make some functions to be force inlined no matter what" ON)
-option(OMATH_VMPROTECT_INTEGRATION
-       "omath will use vmprotect sdk to protect sensitive parts of code from reverse engineering"
-       OFF)
-
 if(VCPKG_MANIFEST_FEATURES)
     foreach(omath_feature IN LISTS VCPKG_MANIFEST_FEATURES)
         if(omath_feature STREQUAL "imgui")
@@ -47,8 +43,6 @@ if(VCPKG_MANIFEST_FEATURES)
             set(OMATH_BUILD_BENCHMARK ON)
         elseif(omath_feature STREQUAL "examples")
             set(OMATH_BUILD_EXAMPLES ON)
-        elseif(omath_feature STREQUAL "vmprotect")
-            set(OMATH_VMPROTECT_INTEGRATION ON)
         endif()
 
     endforeach()
@@ -113,11 +107,6 @@ if(OMATH_IMGUI_INTEGRATION)
 
 endif()
 
-if(OMATH_VMPROTECT_INTEGRATION)
-    find_package(vmprotect_sdk CONFIG REQUIRED)
-    target_link_libraries(${PROJECT_NAME} PUBLIC vmprotect_sdk::vmprotect_sdk)
-endif()
-
 if(OMATH_USE_AVX2)
     target_compile_definitions(${PROJECT_NAME} PUBLIC OMATH_USE_AVX2)
 endif()

CMakePresets.json

@@ -145,7 +145,7 @@
       "hidden": true,
       "inherits": ["linux-base", "vcpkg-base"],
       "cacheVariables": {
-        "VCPKG_MANIFEST_FEATURES": "tests;imgui;avx2;vmprotect;examples"
+        "VCPKG_MANIFEST_FEATURES": "tests;imgui;avx2"
       }
     },
     {

benchmark/benchmark_collision.cpp

@@ -0,0 +1,161 @@
+//
+// Created by Vlad on 3/2/2026.
+//
+#include <benchmark/benchmark.h>
+#include <memory_resource>
+#include <omath/collision/epa_algorithm.hpp>
+#include <omath/collision/gjk_algorithm.hpp>
+#include <omath/engines/source_engine/collider.hpp>
+#include <omath/engines/source_engine/mesh.hpp>
+
+using Mesh = omath::source_engine::Mesh;
+using Collider = omath::source_engine::MeshCollider;
+using Gjk = omath::collision::GjkAlgorithm<Collider>;
+using Epa = omath::collision::Epa<Collider>;
+
+namespace
+{
+    // Unit cube with half-extent 1 — 8 vertices in [-1,1]^3.
+    const std::vector<omath::primitives::Vertex<>> k_cube_vbo = {
+        { { -1.f, -1.f, -1.f }, {}, {} },
+        { { -1.f, -1.f,  1.f }, {}, {} },
+        { { -1.f,  1.f, -1.f }, {}, {} },
+        { { -1.f,  1.f,  1.f }, {}, {} },
+        { {  1.f,  1.f,  1.f }, {}, {} },
+        { {  1.f,  1.f, -1.f }, {}, {} },
+        { {  1.f, -1.f,  1.f }, {}, {} },
+        { {  1.f, -1.f, -1.f }, {}, {} },
+    };
+    const std::vector<omath::Vector3<std::uint32_t>> k_empty_vao{};
+} // namespace
+
+// ---------------------------------------------------------------------------
+// GJK benchmarks
+// ---------------------------------------------------------------------------
+
+// Separated cubes — origin distance 2.1, no overlap.
+// Exercises the early-exit path and the centroid-based initial direction.
+static void BM_Gjk_Separated(benchmark::State& state)
+{
+    const Collider a{Mesh{k_cube_vbo, k_empty_vao}};
+
+    Mesh mesh_b{k_cube_vbo, k_empty_vao};
+    mesh_b.set_origin({0.f, 2.1f, 0.f});
+    const Collider b{mesh_b};
+
+    for ([[maybe_unused]] auto _ : state)
+        benchmark::DoNotOptimize(Gjk::is_collide(a, b));
+}
+
+// Overlapping cubes — B offset by 0.5 along X, ~1.5 units penetration depth.
+static void BM_Gjk_Overlapping(benchmark::State& state)
+{
+    const Collider a{Mesh{k_cube_vbo, k_empty_vao}};
+
+    Mesh mesh_b{k_cube_vbo, k_empty_vao};
+    mesh_b.set_origin({0.5f, 0.f, 0.f});
+    const Collider b{mesh_b};
+
+    for ([[maybe_unused]] auto _ : state)
+        benchmark::DoNotOptimize(Gjk::is_collide(a, b));
+}
+
+// Identical cubes at the same origin — deep overlap / worst case for GJK.
+static void BM_Gjk_SameOrigin(benchmark::State& state)
+{
+    const Collider a{Mesh{k_cube_vbo, k_empty_vao}};
+    const Collider b{Mesh{k_cube_vbo, k_empty_vao}};
+
+    for ([[maybe_unused]] auto _ : state)
+        benchmark::DoNotOptimize(Gjk::is_collide(a, b));
+}
+
+// ---------------------------------------------------------------------------
+// EPA benchmarks
+// ---------------------------------------------------------------------------
+
+// EPA with a pre-allocated monotonic buffer (reset each iteration).
+// Isolates algorithmic cost from allocator overhead.
+static void BM_Epa_MonotonicBuffer(benchmark::State& state)
+{
+    const Collider a{Mesh{k_cube_vbo, k_empty_vao}};
+
+    Mesh mesh_b{k_cube_vbo, k_empty_vao};
+    mesh_b.set_origin({0.5f, 0.f, 0.f});
+    const Collider b{mesh_b};
+
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    if (!hit)
+        return; // shouldn't happen, but guard for safety
+
+    constexpr Epa::Params params{.max_iterations = 64, .tolerance = 1e-4f};
+
+    // Pre-allocate a 32 KiB stack buffer — enough for typical polytope growth.
+    constexpr std::size_t k_buf_size = 32768;
+    alignas(std::max_align_t) char buf[k_buf_size];
+    std::pmr::monotonic_buffer_resource mr{buf, k_buf_size, std::pmr::null_memory_resource()};
+
+    for ([[maybe_unused]] auto _ : state)
+    {
+        mr.release(); // reset the buffer without touching the upstream resource
+        benchmark::DoNotOptimize(Epa::solve(a, b, simplex, params, mr));
+    }
+}
+
+// EPA with the default (malloc-backed) memory resource.
+// Shows total cost including allocator pressure.
+static void BM_Epa_DefaultResource(benchmark::State& state)
+{
+    const Collider a{Mesh{k_cube_vbo, k_empty_vao}};
+
+    Mesh mesh_b{k_cube_vbo, k_empty_vao};
+    mesh_b.set_origin({0.5f, 0.f, 0.f});
+    const Collider b{mesh_b};
+
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    if (!hit)
+        return;
+
+    constexpr Epa::Params params{.max_iterations = 64, .tolerance = 1e-4f};
+
+    for ([[maybe_unused]] auto _ : state)
+        benchmark::DoNotOptimize(Epa::solve(a, b, simplex, params));
+}
+
+// ---------------------------------------------------------------------------
+// Combined GJK + EPA pipeline
+// ---------------------------------------------------------------------------
+
+// Full collision pipeline: GJK detects contact, EPA resolves penetration.
+// This is the hot path in a physics engine tick.
+static void BM_GjkEpa_Pipeline(benchmark::State& state)
+{
+    const Collider a{Mesh{k_cube_vbo, k_empty_vao}};
+
+    Mesh mesh_b{k_cube_vbo, k_empty_vao};
+    mesh_b.set_origin({0.5f, 0.f, 0.f});
+    const Collider b{mesh_b};
+
+    constexpr Epa::Params params{.max_iterations = 64, .tolerance = 1e-4f};
+
+    constexpr std::size_t k_buf_size = 32768;
+    alignas(std::max_align_t) char buf[k_buf_size];
+    std::pmr::monotonic_buffer_resource mr{buf, k_buf_size, std::pmr::null_memory_resource()};
+
+    for ([[maybe_unused]] auto _ : state)
+    {
+        mr.release();
+        const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+        if (hit)
+            benchmark::DoNotOptimize(Epa::solve(a, b, simplex, params, mr));
+    }
+}
+
+BENCHMARK(BM_Gjk_Separated)->Iterations(100'000);
+BENCHMARK(BM_Gjk_Overlapping)->Iterations(100'000);
+BENCHMARK(BM_Gjk_SameOrigin)->Iterations(100'000);
+
+BENCHMARK(BM_Epa_MonotonicBuffer)->Iterations(100'000);
+BENCHMARK(BM_Epa_DefaultResource)->Iterations(100'000);
+
+BENCHMARK(BM_GjkEpa_Pipeline)->Iterations(100'000);

examples/CMakeLists.txt

@@ -2,7 +2,7 @@ add_subdirectory(example_barycentric)
 add_subdirectory(example_glfw3)
 add_subdirectory(example_proj_mat_builder)
 add_subdirectory(example_signature_scan)
-add_subdirectory(exmple_var_encryption)
+
 if(OMATH_ENABLE_VALGRIND)
     omath_setup_valgrind(example_projection_matrix_builder)
     omath_setup_valgrind(example_signature_scan)

examples/exmple_var_encryption/CMakeLists.txt

@@ -1,10 +0,0 @@
-project(example_var_encryption)
-
-add_executable(${PROJECT_NAME} main.cpp)
-set_target_properties(
-    ${PROJECT_NAME}
-    PROPERTIES CXX_STANDARD 23
-               ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/out/${CMAKE_BUILD_TYPE}"
-               LIBRARY_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/out/${CMAKE_BUILD_TYPE}"
-               RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/out/${CMAKE_BUILD_TYPE}")
-target_link_libraries(${PROJECT_NAME} PRIVATE omath::omath)

examples/exmple_var_encryption/main.cpp

@@ -1,15 +0,0 @@
-//
-// Created by orange on 24.02.2026.
-//
-#include "omath/containers/encrypted_variable.hpp"
-#include <omath/omath.hpp>
-#include <print>
-int main()
-{
-    OMATH_DEF_CRYPT_VAR(int, 64) var{5};
-    var.encrypt();
-    std::println("{}", var.value());
-    var.decrypt();
-    std::println("{}", var.value());
-    return var.value();
-}

include/omath/collision/epa_algorithm.hpp

@@ -8,6 +8,7 @@
 #include <memory>
 #include <memory_resource>
 #include <queue>
+#include <unordered_map>
 #include <utility>
 #include <vector>
 
@@ -56,83 +57,76 @@ namespace omath::collision
                                            const Simplex<VectorType>& simplex, const Params params = {},
                                            std::pmr::memory_resource& mem_resource = *std::pmr::get_default_resource())
         {
-            // --- Build initial polytope from simplex (4 points) ---
             std::pmr::vector<VectorType> vertexes = build_initial_polytope_from_simplex(simplex, mem_resource);
-
-            // Initial tetra faces (windings corrected in make_face)
             std::pmr::vector<Face> faces = create_initial_tetra_faces(mem_resource, vertexes);
 
-            auto heap = rebuild_heap(faces, mem_resource);
+            // Build initial min-heap by distance.
+            Heap heap = rebuild_heap(faces, mem_resource);
 
             Result out{};
 
+            // Hoisted outside the loop to reuse bucket allocation across iterations.
+            // Initial bucket count 16 covers a typical horizon without rehashing.
+            BoundaryMap boundary{16, &mem_resource};
+
             for (int it = 0; it < params.max_iterations; ++it)
             {
-                // If heap might be stale after face edits, rebuild lazily.
+                // Lazily discard stale (deleted or index-mismatched) heap entries.
-                if (heap.empty())
+                discard_stale_heap_entries(faces, heap);
-                    break;
-                // Rebuild when the "closest" face changed (simple cheap guard)
-                // (We could keep face handles; this is fine for small Ns.)
-
-                if (const auto top = heap.top(); faces[top.idx].d != top.d)
-                    heap = rebuild_heap(faces, mem_resource);
 
                 if (heap.empty())
                     break;
 
-                // FIXME: STORE REF VALUE, DO NOT USE
-                //  AFTER IF STATEMENT BLOCK
                 const Face& face = faces[heap.top().idx];
 
-                // Get the furthest point in face normal direction
                 const VectorType p = support_point(a, b, face.n);
                 const auto p_dist = face.n.dot(p);
 
-                // Converged if we can’t push the face closer than tolerance
+                // Converged: new support can't push the face closer than tolerance.
                 if (p_dist - face.d <= params.tolerance)
                 {
                     out.normal = face.n;
-                    out.depth = face.d; // along unit normal
+                    out.depth = face.d;
                     out.iterations = it + 1;
                     out.num_vertices = static_cast<int>(vertexes.size());
                     out.num_faces = static_cast<int>(faces.size());
-
                     out.penetration_vector = out.normal * out.depth;
                     return out;
                 }
 
-                // Add new vertex
                 const int new_idx = static_cast<int>(vertexes.size());
                 vertexes.emplace_back(p);
 
-                const auto [to_delete, boundary] = mark_visible_and_collect_horizon(faces, p);
+                // Tombstone visible faces and collect the horizon boundary.
+                // This avoids copying the faces array (O(n)) each iteration.
+                tombstone_visible_faces(faces, boundary, p);
 
-                erase_marked(faces, to_delete);
+                // Stitch new faces around the horizon and push them directly onto the
-
+                // heap — no full O(n log n) rebuild needed.
-                // Stitch new faces around the horizon
+                for (const auto& [key, e] : boundary)
-                for (const auto& e : boundary)
+                {
+                    const int fi = static_cast<int>(faces.size());
                     faces.emplace_back(make_face(vertexes, e.a, e.b, new_idx));
-
+                    heap.emplace(faces.back().d, fi);
-                // Rebuild heap after topology change
+                }
-                heap = rebuild_heap(faces, mem_resource);
 
                 if (!std::isfinite(vertexes.back().dot(vertexes.back())))
                     break; // safety
+
                 out.iterations = it + 1;
             }
 
-            if (faces.empty())
+            // Find the best surviving (non-deleted) face.
+            const Face* best = find_best_surviving_face(faces);
+
+            if (!best)
                 return std::nullopt;
 
-            const auto best = *std::ranges::min_element(faces, [](const auto& first, const auto& second)
+            out.normal = best->n;
-                                                        { return first.d < second.d; });
+            out.depth = best->d;
-            out.normal = best.n;
-            out.depth = best.d;
             out.num_vertices = static_cast<int>(vertexes.size());
             out.num_faces = static_cast<int>(faces.size());
-
             out.penetration_vector = out.normal * out.depth;
-
             return out;
         }
 
@@ -141,7 +135,8 @@ namespace omath::collision
         {
             int i0, i1, i2;
             VectorType n; // unit outward normal
-            FloatingType d; // n · v0  (>=0 ideally because origin is inside)
+            FloatingType d; // n · v0  (>= 0 ideally because origin is inside)
+            bool deleted{false}; // tombstone flag — avoids O(n) compaction per iteration
         };
 
         struct Edge final
@@ -154,6 +149,7 @@ namespace omath::collision
             FloatingType d;
             int idx;
         };
+
         struct HeapCmp final
         {
             [[nodiscard]]
@@ -165,35 +161,44 @@ namespace omath::collision
 
         using Heap = std::priority_queue<HeapItem, std::pmr::vector<HeapItem>, HeapCmp>;
 
+        // Horizon boundary: maps packed(a,b) → Edge.
+        // Opposite edges cancel in O(1) via hash lookup instead of O(h) linear scan.
+        using BoundaryMap = std::pmr::unordered_map<std::int64_t, Edge>;
+
+        [[nodiscard]]
+        static constexpr std::int64_t pack_edge(const int a, const int b) noexcept
+        {
+            return (static_cast<std::int64_t>(a) << 32) | static_cast<std::uint32_t>(b);
+        }
+
         [[nodiscard]]
         static Heap rebuild_heap(const std::pmr::vector<Face>& faces, auto& memory_resource)
         {
             std::pmr::vector<HeapItem> storage{&memory_resource};
-            storage.reserve(faces.size()); // optional but recommended
+            storage.reserve(faces.size());
-
             Heap h{HeapCmp{}, std::move(storage)};
-
             for (int i = 0; i < static_cast<int>(faces.size()); ++i)
-                h.emplace(faces[i].d, i);
+                if (!faces[i].deleted)
-
+                    h.emplace(faces[i].d, i);
-            return h; // allocator is preserved
+            return h;
         }
 
         [[nodiscard]]
         static bool visible_from(const Face& f, const VectorType& p)
         {
-            // positive if p is in front of the face
             return f.n.dot(p) - f.d > static_cast<FloatingType>(1e-7);
         }
 
-        static void add_edge_boundary(std::pmr::vector<Edge>& boundary, int a, int b)
+        static void add_edge_boundary(BoundaryMap& boundary, int a, int b)
         {
-            // Keep edges that appear only once; erase if opposite already present
+            // O(1) cancel: if the opposite edge (b→a) is already in the map it is an
-            auto itb = std::ranges::find_if(boundary, [&](const Edge& e) { return e.a == b && e.b == a; });
+            // internal edge shared by two visible faces and must be removed.
-            if (itb != boundary.end())
+            // Otherwise this is a horizon edge and we insert it.
-                boundary.erase(itb); // internal edge cancels out
+            const std::int64_t rev = pack_edge(b, a);
+            if (const auto it = boundary.find(rev); it != boundary.end())
+                boundary.erase(it);
             else
-                boundary.emplace_back(a, b); // horizon edge (directed)
+                boundary.emplace(pack_edge(a, b), Edge{a, b});
         }
 
         [[nodiscard]]
@@ -204,9 +209,7 @@ namespace omath::collision
             const VectorType& a2 = vertexes[i2];
             VectorType n = (a1 - a0).cross(a2 - a0);
             if (n.dot(n) <= static_cast<FloatingType>(1e-30))
-            {
                 n = any_perp_vec(a1 - a0); // degenerate guard
-            }
             // Ensure normal points outward (away from origin): require n·a0 >= 0
             if (n.dot(a0) < static_cast<FloatingType>(0.0))
             {
@@ -243,6 +246,7 @@ namespace omath::collision
                     return d;
             return V{1, 0, 0};
         }
+
         [[nodiscard]]
         static std::pmr::vector<Face> create_initial_tetra_faces(std::pmr::memory_resource& mem_resource,
                                                                  const std::pmr::vector<VectorType>& vertexes)
@@ -262,48 +266,45 @@ namespace omath::collision
         {
             std::pmr::vector<VectorType> vertexes{&mem_resource};
             vertexes.reserve(simplex.size());
-
             for (std::size_t i = 0; i < simplex.size(); ++i)
                 vertexes.emplace_back(simplex[i]);
-
             return vertexes;
         }
-        static void erase_marked(std::pmr::vector<Face>& faces, const std::pmr::vector<bool>& to_delete)
+
+        static const Face* find_best_surviving_face(const std::pmr::vector<Face>& faces)
         {
-            auto* mr = faces.get_allocator().resource(); // keep same resource
+            const Face* best = nullptr;
-            std::pmr::vector<Face> kept{mr};
+            for (const auto& f : faces)
-            kept.reserve(faces.size());
+                if (!f.deleted && (best == nullptr || f.d < best->d))
-
+                    best = &f;
-            for (std::size_t i = 0; i < faces.size(); ++i)
+            return best;
-                if (!to_delete[i])
-                    kept.emplace_back(faces[i]);
-
-            faces.swap(kept);
         }
-        struct Horizon
+        static void tombstone_visible_faces(std::pmr::vector<Face>& faces, BoundaryMap& boundary,
+                                            const VectorType& p)
         {
-            std::pmr::vector<bool> to_delete;
+            boundary.clear();
-            std::pmr::vector<Edge> boundary;
+            for (auto& f : faces)
-        };
+            {
-
+                if (!f.deleted && visible_from(f, p))
-        static Horizon mark_visible_and_collect_horizon(const std::pmr::vector<Face>& faces, const VectorType& p)
-        {
-            auto* mr = faces.get_allocator().resource();
-
-            Horizon horizon{std::pmr::vector<bool>(faces.size(), false, mr), std::pmr::vector<Edge>(mr)};
-            horizon.boundary.reserve(faces.size());
-
-            for (std::size_t i = 0; i < faces.size(); ++i)
-                if (visible_from(faces[i], p))
                 {
-                    const auto& rf = faces[i];
+                    f.deleted = true;
-                    horizon.to_delete[i] = true;
+                    add_edge_boundary(boundary, f.i0, f.i1);
-                    add_edge_boundary(horizon.boundary, rf.i0, rf.i1);
+                    add_edge_boundary(boundary, f.i1, f.i2);
-                    add_edge_boundary(horizon.boundary, rf.i1, rf.i2);
+                    add_edge_boundary(boundary, f.i2, f.i0);
-                    add_edge_boundary(horizon.boundary, rf.i2, rf.i0);
                 }
+            }
+        }
 
-            return horizon;
+        static void discard_stale_heap_entries(const std::pmr::vector<Face>& faces,
+                                               std::priority_queue<HeapItem, std::pmr::vector<HeapItem>, HeapCmp>& heap)
+        {
+            while (!heap.empty())
+            {
+                const auto& top = heap.top();
+                if (!faces[top.idx].deleted && faces[top.idx].d == top.d)
+                    break;
+                heap.pop();
+            }
         }
     };
 } // namespace omath::collision

include/omath/collision/gjk_algorithm.hpp

@@ -14,11 +14,15 @@ namespace omath::collision
         Simplex<VertexType> simplex; // valid only if hit == true and size==4
     };
 
+    struct GjkSettings final
+    {
+        float epsilon = 1e-6f;
+        std::size_t max_iterations = 64;
+    };
     template<class ColliderInterfaceType>
     class GjkAlgorithm final
     {
         using VectorType = ColliderInterfaceType::VectorType;
-
     public:
         [[nodiscard]]
         static VectorType find_support_vertex(const ColliderInterfaceType& collider_a,
@@ -36,20 +40,34 @@ namespace omath::collision
 
         [[nodiscard]]
         static GjkHitInfo<VectorType> is_collide_with_simplex_info(const ColliderInterfaceType& collider_a,
-                                                                   const ColliderInterfaceType& collider_b)
+                                                                   const ColliderInterfaceType& collider_b,
+                                                                   const GjkSettings& settings = {})
         {
-            auto support = find_support_vertex(collider_a, collider_b, VectorType{1, 0, 0});
+            // Use centroid difference as initial direction — greatly reduces iterations for separated shapes.
+            VectorType initial_dir;
+            if constexpr (requires { collider_b.get_origin() - collider_a.get_origin(); })
+            {
+                initial_dir = collider_b.get_origin() - collider_a.get_origin();
+                if (initial_dir.dot(initial_dir) < settings.epsilon * settings.epsilon)
+                    initial_dir = VectorType{1, 0, 0};
+            }
+            else
+            {
+                initial_dir = VectorType{1, 0, 0};
+            }
+
+            auto support = find_support_vertex(collider_a, collider_b, initial_dir);
 
             Simplex<VectorType> simplex;
             simplex.push_front(support);
 
             auto direction = -support;
 
-            while (true)
+            for (std::size_t iteration = 0; iteration < settings.max_iterations; ++iteration)
             {
                 support = find_support_vertex(collider_a, collider_b, direction);
 
-                if (support.dot(direction) <= 0.f)
+                if (support.dot(direction) <= settings.epsilon)
                     return {false, simplex};
 
                 simplex.push_front(support);
@@ -57,6 +75,7 @@ namespace omath::collision
                 if (simplex.handle(direction))
                     return {true, simplex};
             }
+            return {false, simplex};
         }
     };
 } // namespace omath::collision

include/omath/collision/mesh_collider.hpp

@@ -42,13 +42,40 @@ namespace omath::collision
             m_mesh.set_origin(new_origin);
         }
 
+        [[nodiscard]]
+        const MeshType& get_mesh() const
+        {
+            return m_mesh;
+        }
+        [[nodiscard]]
+        MeshType& get_mesh()
+        {
+            return m_mesh;
+        }
     private:
         [[nodiscard]]
         const VertexType& find_furthest_vertex(const VectorType& direction) const
         {
-            return *std::ranges::max_element(
+            // The support query arrives in world space, but vertex positions are stored
-                    m_mesh.m_vertex_buffer, [&direction](const auto& first, const auto& second)
+            // in local space.  We need argmax_v { world(v) · d }.
-                    { return first.position.dot(direction) < second.position.dot(direction); });
+            //
+            // world(v) = M·v  (ignoring translation, which is constant across vertices)
+            // world(v) · d = v · Mᵀ·d
+            //
+            // So we transform the direction to local space once — O(1) — then compare
+            // raw local positions, which is far cheaper than calling
+            // vertex_position_to_world_space (full 4×4 multiply) for every vertex.
+            //
+            // d_local = upper-left 3×3 of M, transposed, times d_world:
+            //   d_local[j] = sum_i  M.at(i,j) * d[i]   (i.e. column j of M dotted with d)
+            const auto& m = m_mesh.get_to_world_matrix();
+            const VectorType d_local = {
+                    m[0, 0] * direction.x + m[1, 0] * direction.y + m[2, 0] * direction.z,
+                    m[0, 1] * direction.x + m[1, 1] * direction.y + m[2, 1] * direction.z,
+                    m[0, 2] * direction.x + m[1, 2] * direction.y + m[2, 2] * direction.z,
+            };
+            return *std::ranges::max_element(m_mesh.m_vertex_buffer, [&d_local](const auto& first, const auto& second)
+                                             { return first.position.dot(d_local) < second.position.dot(d_local); });
         }
         MeshType m_mesh;
     };

include/omath/containers/encrypted_variable.hpp

@@ -2,14 +2,10 @@
 // Created by Vladislav on 04.01.2026.
 //
 #pragma once
-#include <VMProtectSDK.h>
 #include <array>
 #include <cstddef>
 #include <cstdint>
 #include <span>
-#include <source_location>
-
-
 #ifdef OMATH_ENABLE_FORCE_INLINE
 #ifdef _MSC_VER
 #define OMATH_FORCE_INLINE __forceinline
@@ -114,18 +110,16 @@ namespace omath
         bool m_is_encrypted{};
         value_type m_data{};
 
-        OMATH_FORCE_INLINE void xor_contained_var_by_key()
+        OMATH_FORCE_INLINE constexpr void xor_contained_var_by_key()
         {
-            VMProtectBeginVirtualization(nullptr);
             // Safe, keeps const-correctness, and avoids reinterpret_cast issues
             auto bytes = std::as_writable_bytes(std::span<value_type, 1>{&m_data, 1});
 
             for (std::size_t i = 0; i < bytes.size(); ++i)
             {
-                const auto k = static_cast<std::uint8_t>(key[i % key_size] + (i * key_size));
+                const std::uint8_t k = static_cast<std::uint8_t>(key[i % key_size] + (i * key_size));
                 bytes[i] ^= static_cast<std::byte>(k);
             }
-            VMProtectEnd();
         }
 
     public:
@@ -140,7 +134,7 @@ namespace omath
             return m_is_encrypted;
         }
 
-        OMATH_FORCE_INLINE void decrypt()
+        OMATH_FORCE_INLINE constexpr void decrypt()
         {
             if (!m_is_encrypted)
                 return;
@@ -148,7 +142,7 @@ namespace omath
             m_is_encrypted = false;
         }
 
-        OMATH_FORCE_INLINE void encrypt()
+        OMATH_FORCE_INLINE constexpr void encrypt()
         {
             if (m_is_encrypted)
                 return;

include/omath/linear_algebra/quaternion.hpp

@@ -0,0 +1,219 @@
+//
+// Created by vlad on 3/1/2026.
+//
+#pragma once
+
+#include "omath/linear_algebra/mat.hpp"
+#include "omath/linear_algebra/vector3.hpp"
+#include <array>
+#include <cmath>
+#include <format>
+
+namespace omath
+{
+    template<class Type>
+    requires std::is_arithmetic_v<Type>
+    class Quaternion
+    {
+    public:
+        using ContainedType = Type;
+
+        Type x = static_cast<Type>(0);
+        Type y = static_cast<Type>(0);
+        Type z = static_cast<Type>(0);
+        Type w = static_cast<Type>(1);  // identity quaternion
+
+        constexpr Quaternion() noexcept = default;
+
+        constexpr Quaternion(const Type& x, const Type& y, const Type& z, const Type& w) noexcept
+            : x(x), y(y), z(z), w(w)
+        {
+        }
+
+        // Factory: build from a normalized axis and an angle in radians
+        [[nodiscard]]
+        static Quaternion from_axis_angle(const Vector3<Type>& axis, const Type& angle_rad) noexcept
+        {
+            const Type half = angle_rad / static_cast<Type>(2);
+            const Type s = std::sin(half);
+            return {axis.x * s, axis.y * s, axis.z * s, std::cos(half)};
+        }
+
+        [[nodiscard]] constexpr bool operator==(const Quaternion& other) const noexcept
+        {
+            return x == other.x && y == other.y && z == other.z && w == other.w;
+        }
+
+        [[nodiscard]] constexpr bool operator!=(const Quaternion& other) const noexcept
+        {
+            return !(*this == other);
+        }
+
+        // Hamilton product: this * other
+        [[nodiscard]] constexpr Quaternion operator*(const Quaternion& other) const noexcept
+        {
+            return {
+                w * other.x + x * other.w + y * other.z - z * other.y,
+                w * other.y - x * other.z + y * other.w + z * other.x,
+                w * other.z + x * other.y - y * other.x + z * other.w,
+                w * other.w - x * other.x - y * other.y - z * other.z,
+            };
+        }
+
+        constexpr Quaternion& operator*=(const Quaternion& other) noexcept
+        {
+            return *this = *this * other;
+        }
+
+        [[nodiscard]] constexpr Quaternion operator*(const Type& scalar) const noexcept
+        {
+            return {x * scalar, y * scalar, z * scalar, w * scalar};
+        }
+
+        constexpr Quaternion& operator*=(const Type& scalar) noexcept
+        {
+            x *= scalar;
+            y *= scalar;
+            z *= scalar;
+            w *= scalar;
+            return *this;
+        }
+
+        [[nodiscard]] constexpr Quaternion operator+(const Quaternion& other) const noexcept
+        {
+            return {x + other.x, y + other.y, z + other.z, w + other.w};
+        }
+
+        constexpr Quaternion& operator+=(const Quaternion& other) noexcept
+        {
+            x += other.x;
+            y += other.y;
+            z += other.z;
+            w += other.w;
+            return *this;
+        }
+
+        [[nodiscard]] constexpr Quaternion operator-() const noexcept
+        {
+            return {-x, -y, -z, -w};
+        }
+
+        // Conjugate: negates the vector part (x, y, z)
+        [[nodiscard]] constexpr Quaternion conjugate() const noexcept
+        {
+            return {-x, -y, -z, w};
+        }
+
+        [[nodiscard]] constexpr Type dot(const Quaternion& other) const noexcept
+        {
+            return x * other.x + y * other.y + z * other.z + w * other.w;
+        }
+
+        [[nodiscard]] constexpr Type length_sqr() const noexcept
+        {
+            return x * x + y * y + z * z + w * w;
+        }
+
+#ifndef _MSC_VER
+        [[nodiscard]] constexpr Type length() const noexcept
+        {
+            return std::sqrt(length_sqr());
+        }
+
+        [[nodiscard]] constexpr Quaternion normalized() const noexcept
+        {
+            const Type len = length();
+            return len != static_cast<Type>(0) ? *this * (static_cast<Type>(1) / len) : *this;
+        }
+#else
+        [[nodiscard]] Type length() const noexcept
+        {
+            return std::sqrt(length_sqr());
+        }
+
+        [[nodiscard]] Quaternion normalized() const noexcept
+        {
+            const Type len = length();
+            return len != static_cast<Type>(0) ? *this * (static_cast<Type>(1) / len) : *this;
+        }
+#endif
+
+        // Inverse: q* / |q|^2  (for unit quaternions inverse == conjugate)
+        [[nodiscard]] constexpr Quaternion inverse() const noexcept
+        {
+            return conjugate() * (static_cast<Type>(1) / length_sqr());
+        }
+
+        // Rotate a 3D vector: v' = q * pure(v) * q^-1
+        // Computed via Rodrigues' formula to avoid full quaternion product overhead
+        [[nodiscard]] constexpr Vector3<Type> rotate(const Vector3<Type>& v) const noexcept
+        {
+            const Vector3<Type> q_vec{x, y, z};
+            const Vector3<Type> cross = q_vec.cross(v);
+            return v + cross * (static_cast<Type>(2) * w) + q_vec.cross(cross) * static_cast<Type>(2);
+        }
+
+        // 3x3 rotation matrix from this (unit) quaternion
+        [[nodiscard]] constexpr Mat<3, 3, Type> to_rotation_matrix3() const noexcept
+        {
+            const Type xx = x * x, yy = y * y, zz = z * z;
+            const Type xy = x * y, xz = x * z, yz = y * z;
+            const Type wx = w * x, wy = w * y, wz = w * z;
+            const Type one = static_cast<Type>(1);
+            const Type two = static_cast<Type>(2);
+
+            return {
+                {one - two * (yy + zz), two * (xy - wz),       two * (xz + wy)      },
+                {two * (xy + wz),       one - two * (xx + zz), two * (yz - wx)      },
+                {two * (xz - wy),       two * (yz + wx),       one - two * (xx + yy)},
+            };
+        }
+
+        // 4x4 rotation matrix (with homogeneous row/column)
+        [[nodiscard]] constexpr Mat<4, 4, Type> to_rotation_matrix4() const noexcept
+        {
+            const Type xx = x * x, yy = y * y, zz = z * z;
+            const Type xy = x * y, xz = x * z, yz = y * z;
+            const Type wx = w * x, wy = w * y, wz = w * z;
+            const Type one = static_cast<Type>(1);
+            const Type two = static_cast<Type>(2);
+            const Type zero = static_cast<Type>(0);
+
+            return {
+                {one - two * (yy + zz), two * (xy - wz),       two * (xz + wy),       zero},
+                {two * (xy + wz),       one - two * (xx + zz), two * (yz - wx),       zero},
+                {two * (xz - wy),       two * (yz + wx),       one - two * (xx + yy), zero},
+                {zero,                  zero,                   zero,                  one },
+            };
+        }
+
+        [[nodiscard]] constexpr std::array<Type, 4> as_array() const noexcept
+        {
+            return {x, y, z, w};
+        }
+    };
+} // namespace omath
+
+template<class Type>
+struct std::formatter<omath::Quaternion<Type>> // NOLINT(*-dcl58-cpp)
+{
+    [[nodiscard]]
+    static constexpr auto parse(std::format_parse_context& ctx)
+    {
+        return ctx.begin();
+    }
+
+    template<class FormatContext>
+    [[nodiscard]]
+    static auto format(const omath::Quaternion<Type>& q, FormatContext& ctx)
+    {
+        if constexpr (std::is_same_v<typename FormatContext::char_type, char>)
+            return std::format_to(ctx.out(), "[{}, {}, {}, {}]", q.x, q.y, q.z, q.w);
+
+        if constexpr (std::is_same_v<typename FormatContext::char_type, wchar_t>)
+            return std::format_to(ctx.out(), L"[{}, {}, {}, {}]", q.x, q.y, q.z, q.w);
+
+        if constexpr (std::is_same_v<typename FormatContext::char_type, char8_t>)
+            return std::format_to(ctx.out(), u8"[{}, {}, {}, {}]", q.x, q.y, q.z, q.w);
+    }
+};

include/omath/omath.hpp

@@ -17,6 +17,9 @@
 // Matrix classes
 #include "omath/linear_algebra/mat.hpp"
 
+// Quaternion
+#include "omath/linear_algebra/quaternion.hpp"
+
 // Color functionality
 #include "omath/utility/color.hpp"
 

include/omath/utility/color.hpp

@@ -16,19 +16,28 @@ namespace omath
         float value{};
     };
 
-    class Color final : public Vector4<float>
+    class Color final
     {
+        Vector4<float> m_value;
     public:
-        constexpr Color(const float r, const float g, const float b, const float a) noexcept: Vector4(r, g, b, a)
+        constexpr const Vector4<float>& value() const
         {
-            clamp(0.f, 1.f);
+            return m_value;
+        }
+        constexpr Color(const float r, const float g, const float b, const float a) noexcept: m_value(r, g, b, a)
+        {
+            m_value.clamp(0.f, 1.f);
         }
 
+        constexpr explicit Color(const Vector4<float>& value) : m_value(value)
+        {
+            m_value.clamp(0.f, 1.f);
+        }
         constexpr explicit Color() noexcept = default;
         [[nodiscard]]
         constexpr static Color from_rgba(const uint8_t r, const uint8_t g, const uint8_t b, const uint8_t a) noexcept
         {
-            return Color{Vector4(r, g, b, a) / 255.f};
+            return Color(Vector4<float>(r, g, b, a) / 255.f);
         }
 
         [[nodiscard]]
@@ -82,9 +91,9 @@ namespace omath
         {
             Hsv hsv_data;
 
-            const float& red = x;
+            const float& red = m_value.x;
-            const float& green = y;
+            const float& green = m_value.y;
-            const float& blue = z;
+            const float& blue = m_value.z;
 
             const float max = std::max({red, green, blue});
             const float min = std::min({red, green, blue});
@@ -109,11 +118,6 @@ namespace omath
 
             return hsv_data;
         }
-
-        constexpr explicit Color(const Vector4& vec) noexcept: Vector4(vec)
-        {
-            clamp(0.f, 1.f);
-        }
         constexpr void set_hue(const float hue) noexcept
         {
             auto hsv = to_hsv();
@@ -141,7 +145,7 @@ namespace omath
         constexpr Color blend(const Color& other, float ratio) const noexcept
         {
             ratio = std::clamp(ratio, 0.f, 1.f);
-            return Color(*this * (1.f - ratio) + other * ratio);
+            return Color(this->m_value * (1.f - ratio) + other.m_value * ratio);
         }
 
         [[nodiscard]] static constexpr Color red()
@@ -160,16 +164,26 @@ namespace omath
         [[nodiscard]]
         ImColor to_im_color() const noexcept
         {
-            return {to_im_vec4()};
+            return {m_value.to_im_vec4()};
         }
 #endif
         [[nodiscard]] std::string to_string() const noexcept
         {
             return std::format("[r:{}, g:{}, b:{}, a:{}]",
-                            static_cast<int>(x * 255.f),
+                            static_cast<int>(m_value.x * 255.f),
-                            static_cast<int>(y * 255.f),
+                            static_cast<int>(m_value.y * 255.f),
-                            static_cast<int>(z * 255.f),
+                            static_cast<int>(m_value.z * 255.f),
-                            static_cast<int>(w * 255.f));
+                            static_cast<int>(m_value.w * 255.f));
+        }
+        [[nodiscard]] std::string to_rgbf_string() const noexcept
+        {
+            return std::format("[r:{}, g:{}, b:{}, a:{}]",
+                            m_value.x, m_value.y, m_value.z, m_value.w);
+        }
+        [[nodiscard]] std::string to_hsv_string() const noexcept
+        {
+            const auto [hue, saturation, value] = to_hsv();
+            return std::format("[h:{}, s:{}, v:{}]", hue, saturation, value);
         }
         [[nodiscard]] std::wstring to_wstring() const noexcept
         {
@@ -188,23 +202,55 @@ namespace omath
 template<>
 struct std::formatter<omath::Color> // NOLINT(*-dcl58-cpp)
 {
-    [[nodiscard]]
+    enum class ColorFormat { rgb, rgbf, hsv };
-    static constexpr auto parse(const std::format_parse_context& ctx)
+    ColorFormat color_format = ColorFormat::rgb;
+
+    constexpr auto parse(std::format_parse_context& ctx)
     {
-        return ctx.begin();
+        const auto it = ctx.begin();
+        const auto end = ctx.end();
+
+        if (it == end || *it == '}')
+            return it;
+
+        const std::string_view spec(it, end);
+
+        if (spec.starts_with("rgbf"))
+        {
+            color_format = ColorFormat::rgbf;
+            return it + 4;
+        }
+        if (spec.starts_with("rgb"))
+        {
+            color_format = ColorFormat::rgb;
+            return it + 3;
+        }
+        if (spec.starts_with("hsv"))
+        {
+            color_format = ColorFormat::hsv;
+            return it + 3;
+        }
+
+        throw std::format_error("Invalid format specifier for omath::Color. Use rgb, rgbf, or hsv.");
     }
 
     template<class FormatContext>
-    [[nodiscard]]
+    auto format(const omath::Color& col, FormatContext& ctx) const
-    static auto format(const omath::Color& col, FormatContext& ctx)
     {
-        if constexpr (std::is_same_v<typename FormatContext::char_type, char>)
+        std::string str;
-            return std::format_to(ctx.out(), "{}", col.to_string());
+        switch (color_format)
-        if constexpr (std::is_same_v<typename FormatContext::char_type, wchar_t>)
+        {
-            return std::format_to(ctx.out(), L"{}", col.to_wstring());
+            case ColorFormat::rgb:  str = col.to_string();      break;
+            case ColorFormat::rgbf: str = col.to_rgbf_string(); break;
+            case ColorFormat::hsv:  str = col.to_hsv_string();  break;
+        }
 
+        if constexpr (std::is_same_v<typename FormatContext::char_type, char>)
+            return std::format_to(ctx.out(), "{}", str);
+        if constexpr (std::is_same_v<typename FormatContext::char_type, wchar_t>)
+            return std::format_to(ctx.out(), L"{}", std::wstring(str.cbegin(), str.cend()));
         if constexpr (std::is_same_v<typename FormatContext::char_type, char8_t>)
-            return std::format_to(ctx.out(), u8"{}", col.to_u8string());
+            return std::format_to(ctx.out(), u8"{}", std::u8string(str.cbegin(), str.cend()));
 
         std::unreachable();
     }

source/utility/elf_pattern_scan.cpp

@@ -8,7 +8,6 @@
 #include <utility>
 #include <variant>
 #include <vector>
-#include <VMProtectSDK.h>
 
 #pragma pack(push, 1)
 

tests/general/unit_test_color_grouped.cpp

@@ -26,38 +26,38 @@ protected:
 TEST_F(UnitTestColorGrouped, Constructor_Float)
 {
     constexpr Color color(0.5f, 0.5f, 0.5f, 1.0f);
-    EXPECT_FLOAT_EQ(color.x, 0.5f);
+    EXPECT_FLOAT_EQ(color.value().x, 0.5f);
-    EXPECT_FLOAT_EQ(color.y, 0.5f);
+    EXPECT_FLOAT_EQ(color.value().y, 0.5f);
-    EXPECT_FLOAT_EQ(color.z, 0.5f);
+    EXPECT_FLOAT_EQ(color.value().z, 0.5f);
-    EXPECT_FLOAT_EQ(color.w, 1.0f);
+    EXPECT_FLOAT_EQ(color.value().w, 1.0f);
 }
 
 TEST_F(UnitTestColorGrouped, Constructor_Vector4)
 {
     constexpr omath::Vector4 vec(0.2f, 0.4f, 0.6f, 0.8f);
     constexpr Color color(vec);
-    EXPECT_FLOAT_EQ(color.x, 0.2f);
+    EXPECT_FLOAT_EQ(color.value().x, 0.2f);
-    EXPECT_FLOAT_EQ(color.y, 0.4f);
+    EXPECT_FLOAT_EQ(color.value().y, 0.4f);
-    EXPECT_FLOAT_EQ(color.z, 0.6f);
+    EXPECT_FLOAT_EQ(color.value().z, 0.6f);
-    EXPECT_FLOAT_EQ(color.w, 0.8f);
+    EXPECT_FLOAT_EQ(color.value().w, 0.8f);
 }
 
 TEST_F(UnitTestColorGrouped, FromRGBA)
 {
     constexpr Color color = Color::from_rgba(128, 64, 32, 255);
-    EXPECT_FLOAT_EQ(color.x, 128.0f / 255.0f);
+    EXPECT_FLOAT_EQ(color.value().x, 128.0f / 255.0f);
-    EXPECT_FLOAT_EQ(color.y, 64.0f / 255.0f);
+    EXPECT_FLOAT_EQ(color.value().y, 64.0f / 255.0f);
-    EXPECT_FLOAT_EQ(color.z, 32.0f / 255.0f);
+    EXPECT_FLOAT_EQ(color.value().z, 32.0f / 255.0f);
-    EXPECT_FLOAT_EQ(color.w, 1.0f);
+    EXPECT_FLOAT_EQ(color.value().w, 1.0f);
 }
 
 TEST_F(UnitTestColorGrouped, FromHSV)
 {
     constexpr Color color = Color::from_hsv(0.0f, 1.0f, 1.0f); // Red in HSV
-    EXPECT_FLOAT_EQ(color.x, 1.0f);
+    EXPECT_FLOAT_EQ(color.value().x, 1.0f);
-    EXPECT_FLOAT_EQ(color.y, 0.0f);
+    EXPECT_FLOAT_EQ(color.value().y, 0.0f);
-    EXPECT_FLOAT_EQ(color.z, 0.0f);
+    EXPECT_FLOAT_EQ(color.value().z, 0.0f);
-    EXPECT_FLOAT_EQ(color.w, 1.0f);
+    EXPECT_FLOAT_EQ(color.value().w, 1.0f);
 }
 
 TEST_F(UnitTestColorGrouped, ToHSV)
@@ -71,10 +71,10 @@ TEST_F(UnitTestColorGrouped, ToHSV)
 TEST_F(UnitTestColorGrouped, Blend)
 {
     const Color blended = color1.blend(color2, 0.5f);
-    EXPECT_FLOAT_EQ(blended.x, 0.5f);
+    EXPECT_FLOAT_EQ(blended.value().x, 0.5f);
-    EXPECT_FLOAT_EQ(blended.y, 0.5f);
+    EXPECT_FLOAT_EQ(blended.value().y, 0.5f);
-    EXPECT_FLOAT_EQ(blended.z, 0.0f);
+    EXPECT_FLOAT_EQ(blended.value().z, 0.0f);
-    EXPECT_FLOAT_EQ(blended.w, 1.0f);
+    EXPECT_FLOAT_EQ(blended.value().w, 1.0f);
 }
 
 TEST_F(UnitTestColorGrouped, PredefinedColors)
@@ -83,20 +83,20 @@ TEST_F(UnitTestColorGrouped, PredefinedColors)
     constexpr Color green = Color::green();
     constexpr Color blue = Color::blue();
 
-    EXPECT_FLOAT_EQ(red.x, 1.0f);
+    EXPECT_FLOAT_EQ(red.value().x, 1.0f);
-    EXPECT_FLOAT_EQ(red.y, 0.0f);
+    EXPECT_FLOAT_EQ(red.value().y, 0.0f);
-    EXPECT_FLOAT_EQ(red.z, 0.0f);
+    EXPECT_FLOAT_EQ(red.value().z, 0.0f);
-    EXPECT_FLOAT_EQ(red.w, 1.0f);
+    EXPECT_FLOAT_EQ(red.value().w, 1.0f);
 
-    EXPECT_FLOAT_EQ(green.x, 0.0f);
+    EXPECT_FLOAT_EQ(green.value().x, 0.0f);
-    EXPECT_FLOAT_EQ(green.y, 1.0f);
+    EXPECT_FLOAT_EQ(green.value().y, 1.0f);
-    EXPECT_FLOAT_EQ(green.z, 0.0f);
+    EXPECT_FLOAT_EQ(green.value().z, 0.0f);
-    EXPECT_FLOAT_EQ(green.w, 1.0f);
+    EXPECT_FLOAT_EQ(green.value().w, 1.0f);
 
-    EXPECT_FLOAT_EQ(blue.x, 0.0f);
+    EXPECT_FLOAT_EQ(blue.value().x, 0.0f);
-    EXPECT_FLOAT_EQ(blue.y, 0.0f);
+    EXPECT_FLOAT_EQ(blue.value().y, 0.0f);
-    EXPECT_FLOAT_EQ(blue.z, 1.0f);
+    EXPECT_FLOAT_EQ(blue.value().z, 1.0f);
-    EXPECT_FLOAT_EQ(blue.w, 1.0f);
+    EXPECT_FLOAT_EQ(blue.value().w, 1.0f);
 }
 
 TEST_F(UnitTestColorGrouped, BlendVector3)
@@ -104,9 +104,9 @@ TEST_F(UnitTestColorGrouped, BlendVector3)
     constexpr Color v1(1.0f, 0.0f, 0.0f, 1.f); // Red
     constexpr Color v2(0.0f, 1.0f, 0.0f, 1.f); // Green
     constexpr Color blended = v1.blend(v2, 0.5f);
-    EXPECT_FLOAT_EQ(blended.x, 0.5f);
+    EXPECT_FLOAT_EQ(blended.value().x, 0.5f);
-    EXPECT_FLOAT_EQ(blended.y, 0.5f);
+    EXPECT_FLOAT_EQ(blended.value().y, 0.5f);
-    EXPECT_FLOAT_EQ(blended.z, 0.0f);
+    EXPECT_FLOAT_EQ(blended.value().z, 0.0f);
 }
 
 // From unit_test_color_extra.cpp
@@ -148,37 +148,37 @@ TEST(UnitTestColorGrouped_Extra, BlendEdgeCases)
     constexpr Color a = Color::red();
     constexpr Color b = Color::blue();
     constexpr auto r0 = a.blend(b, 0.f);
-    EXPECT_FLOAT_EQ(r0.x, a.x);
+    EXPECT_FLOAT_EQ(r0.value().x, a.value().x);
     constexpr auto r1 = a.blend(b, 1.f);
-    EXPECT_FLOAT_EQ(r1.x, b.x);
+    EXPECT_FLOAT_EQ(r1.value().x, b.value().x);
 }
 
 // From unit_test_color_more.cpp
 TEST(UnitTestColorGrouped_More, DefaultCtorIsZero)
 {
     constexpr Color c;
-    EXPECT_FLOAT_EQ(c.x, 0.0f);
+    EXPECT_FLOAT_EQ(c.value().x, 0.0f);
-    EXPECT_FLOAT_EQ(c.y, 0.0f);
+    EXPECT_FLOAT_EQ(c.value().y, 0.0f);
-    EXPECT_FLOAT_EQ(c.z, 0.0f);
+    EXPECT_FLOAT_EQ(c.value().z, 0.0f);
-    EXPECT_FLOAT_EQ(c.w, 0.0f);
+    EXPECT_FLOAT_EQ(c.value().w, 0.0f);
 }
 
 TEST(UnitTestColorGrouped_More, FloatCtorAndClampForRGB)
 {
     constexpr Color c(1.2f, -0.5f, 0.5f, 2.0f);
-    EXPECT_FLOAT_EQ(c.x, 1.0f);
+    EXPECT_FLOAT_EQ(c.value().x, 1.0f);
-    EXPECT_FLOAT_EQ(c.y, 0.0f);
+    EXPECT_FLOAT_EQ(c.value().y, 0.0f);
-    EXPECT_FLOAT_EQ(c.z, 0.5f);
+    EXPECT_FLOAT_EQ(c.value().z, 0.5f);
-    EXPECT_FLOAT_EQ(c.w, 2.0f);
+    EXPECT_FLOAT_EQ(c.value().w, 2.0f);
 }
 
 TEST(UnitTestColorGrouped_More, FromRgbaProducesScaledComponents)
 {
     constexpr Color c = Color::from_rgba(25u, 128u, 230u, 64u);
-    EXPECT_NEAR(c.x, 25.0f/255.0f, 1e-6f);
+    EXPECT_NEAR(c.value().x, 25.0f/255.0f, 1e-6f);
-    EXPECT_NEAR(c.y, 128.0f/255.0f, 1e-6f);
+    EXPECT_NEAR(c.value().y, 128.0f/255.0f, 1e-6f);
-    EXPECT_NEAR(c.z, 230.0f/255.0f, 1e-6f);
+    EXPECT_NEAR(c.value().z, 230.0f/255.0f, 1e-6f);
-    EXPECT_NEAR(c.w, 64.0f/255.0f, 1e-6f);
+    EXPECT_NEAR(c.value().w, 64.0f/255.0f, 1e-6f);
 }
 
 TEST(UnitTestColorGrouped_More, BlendProducesIntermediate)
@@ -186,10 +186,10 @@ TEST(UnitTestColorGrouped_More, BlendProducesIntermediate)
     constexpr Color c0(0.0f, 0.0f, 0.0f, 1.0f);
     constexpr Color c1(1.0f, 1.0f, 1.0f, 0.0f);
     constexpr Color mid = c0.blend(c1, 0.5f);
-    EXPECT_FLOAT_EQ(mid.x, 0.5f);
+    EXPECT_FLOAT_EQ(mid.value().x, 0.5f);
-    EXPECT_FLOAT_EQ(mid.y, 0.5f);
+    EXPECT_FLOAT_EQ(mid.value().y, 0.5f);
-    EXPECT_FLOAT_EQ(mid.z, 0.5f);
+    EXPECT_FLOAT_EQ(mid.value().z, 0.5f);
-    EXPECT_FLOAT_EQ(mid.w, 0.5f);
+    EXPECT_FLOAT_EQ(mid.value().w, 0.5f);
 }
 
 TEST(UnitTestColorGrouped_More, HsvRoundTrip)
@@ -197,9 +197,9 @@ TEST(UnitTestColorGrouped_More, HsvRoundTrip)
     constexpr Color red = Color::red();
     const auto hsv = red.to_hsv();
     const Color back = Color::from_hsv(hsv);
-    EXPECT_NEAR(back.x, 1.0f, 1e-6f);
+    EXPECT_NEAR(back.value().x, 1.0f, 1e-6f);
-    EXPECT_NEAR(back.y, 0.0f, 1e-6f);
+    EXPECT_NEAR(back.value().y, 0.0f, 1e-6f);
-    EXPECT_NEAR(back.z, 0.0f, 1e-6f);
+    EXPECT_NEAR(back.value().z, 0.0f, 1e-6f);
 }
 
 TEST(UnitTestColorGrouped_More, ToStringContainsComponents)
@@ -230,18 +230,18 @@ TEST(UnitTestColorGrouped_More2, FromHsvCases)
     auto check_hue = [&](float h) {
         SCOPED_TRACE(::testing::Message() << "h=" << h);
         Color c = Color::from_hsv(h, 1.f, 1.f);
-        EXPECT_TRUE(std::isfinite(c.x));
+        EXPECT_TRUE(std::isfinite(c.value().x));
-        EXPECT_TRUE(std::isfinite(c.y));
+        EXPECT_TRUE(std::isfinite(c.value().y));
-        EXPECT_TRUE(std::isfinite(c.z));
+        EXPECT_TRUE(std::isfinite(c.value().z));
-        EXPECT_GE(c.x, -eps);
+        EXPECT_GE(c.value().x, -eps);
-        EXPECT_LE(c.x, 1.f + eps);
+        EXPECT_LE(c.value().x, 1.f + eps);
-        EXPECT_GE(c.y, -eps);
+        EXPECT_GE(c.value().y, -eps);
-        EXPECT_LE(c.y, 1.f + eps);
+        EXPECT_LE(c.value().y, 1.f + eps);
-        EXPECT_GE(c.z, -eps);
+        EXPECT_GE(c.value().z, -eps);
-        EXPECT_LE(c.z, 1.f + eps);
+        EXPECT_LE(c.value().z, 1.f + eps);
 
-        float mx = std::max({c.x, c.y, c.z});
+        float mx = std::max({c.value().x, c.value().y, c.value().z});
-        float mn = std::min({c.x, c.y, c.z});
+        float mn = std::min({c.value().x, c.value().y, c.value().z});
         EXPECT_GE(mx, 0.999f);
         EXPECT_LE(mn, 1e-3f + 1e-4f);
     };
@@ -261,13 +261,13 @@ TEST(UnitTestColorGrouped_More2, ToHsvAndSetters)
     EXPECT_NEAR(hsv.value, 0.6f, 1e-6f);
 
     c.set_hue(0.0f);
-    EXPECT_TRUE(std::isfinite(c.x));
+    EXPECT_TRUE(std::isfinite(c.value().x));
 
     c.set_saturation(0.0f);
-    EXPECT_TRUE(std::isfinite(c.y));
+    EXPECT_TRUE(std::isfinite(c.value().y));
 
     c.set_value(0.5f);
-    EXPECT_TRUE(std::isfinite(c.z));
+    EXPECT_TRUE(std::isfinite(c.value().z));
 }
 
 TEST(UnitTestColorGrouped_More2, BlendAndStaticColors)
@@ -275,14 +275,14 @@ TEST(UnitTestColorGrouped_More2, BlendAndStaticColors)
     constexpr Color a = Color::red();
     constexpr Color b = Color::blue();
     constexpr auto mid = a.blend(b, 0.5f);
-    EXPECT_GT(mid.x, 0.f);
+    EXPECT_GT(mid.value().x, 0.f);
-    EXPECT_GT(mid.z, 0.f);
+    EXPECT_GT(mid.value().z, 0.f);
 
     constexpr auto all_a = a.blend(b, -1.f);
-    EXPECT_NEAR(all_a.x, a.x, 1e-6f);
+    EXPECT_NEAR(all_a.value().x, a.value().x, 1e-6f);
 
     constexpr auto all_b = a.blend(b, 2.f);
-    EXPECT_NEAR(all_b.z, b.z, 1e-6f);
+    EXPECT_NEAR(all_b.value().z, b.value().z, 1e-6f);
 }
 
 TEST(UnitTestColorGrouped_More2, FormatterUsesToString)
@@ -291,3 +291,35 @@ TEST(UnitTestColorGrouped_More2, FormatterUsesToString)
     const auto formatted = std::format("{}", c);
     EXPECT_NE(formatted.find("r:10"), std::string::npos);
 }
+
+TEST(UnitTestColorGrouped_More2, FormatterRgb)
+{
+    constexpr Color c = Color::from_rgba(255, 128, 0, 64);
+    const auto s = std::format("{:rgb}", c);
+    EXPECT_NE(s.find("r:255"), std::string::npos);
+    EXPECT_NE(s.find("g:128"), std::string::npos);
+    EXPECT_NE(s.find("b:0"), std::string::npos);
+    EXPECT_NE(s.find("a:64"), std::string::npos);
+}
+
+TEST(UnitTestColorGrouped_More2, FormatterRgbf)
+{
+    constexpr Color c(0.5f, 0.25f, 1.0f, 0.75f);
+    const auto s = std::format("{:rgbf}", c);
+    EXPECT_NE(s.find("r:"), std::string::npos);
+    EXPECT_NE(s.find("g:"), std::string::npos);
+    EXPECT_NE(s.find("b:"), std::string::npos);
+    EXPECT_NE(s.find("a:"), std::string::npos);
+    // Values should be in [0,1] float range, not 0-255
+    EXPECT_EQ(s.find("r:127"), std::string::npos);
+    EXPECT_EQ(s.find("r:255"), std::string::npos);
+}
+
+TEST(UnitTestColorGrouped_More2, FormatterHsv)
+{
+    const Color c = Color::red();
+    const auto s = std::format("{:hsv}", c);
+    EXPECT_NE(s.find("h:"), std::string::npos);
+    EXPECT_NE(s.find("s:"), std::string::npos);
+    EXPECT_NE(s.find("v:"), std::string::npos);
+}

tests/general/unit_test_epa_comprehensive.cpp

@@ -0,0 +1,471 @@
+//
+// Comprehensive EPA tests.
+// Covers: all 3 axis directions, multiple depth levels, penetration-vector
+// round-trips, depth monotonicity, symmetry, asymmetric sizes, memory
+// resource variants, tolerance sensitivity, and iteration bookkeeping.
+//
+#include <cmath>
+#include <gtest/gtest.h>
+#include <memory_resource>
+#include <omath/collision/epa_algorithm.hpp>
+#include <omath/collision/gjk_algorithm.hpp>
+#include <omath/engines/source_engine/collider.hpp>
+#include <omath/engines/source_engine/mesh.hpp>
+
+using Mesh     = omath::source_engine::Mesh;
+using Collider = omath::source_engine::MeshCollider;
+using Gjk      = omath::collision::GjkAlgorithm<Collider>;
+using Epa      = omath::collision::Epa<Collider>;
+using Vec3     = omath::Vector3<float>;
+
+namespace
+{
+    const std::vector<omath::primitives::Vertex<>> k_cube_vbo = {
+        { { -1.f, -1.f, -1.f }, {}, {} },
+        { { -1.f, -1.f,  1.f }, {}, {} },
+        { { -1.f,  1.f, -1.f }, {}, {} },
+        { { -1.f,  1.f,  1.f }, {}, {} },
+        { {  1.f,  1.f,  1.f }, {}, {} },
+        { {  1.f,  1.f, -1.f }, {}, {} },
+        { {  1.f, -1.f,  1.f }, {}, {} },
+        { {  1.f, -1.f, -1.f }, {}, {} },
+    };
+    const std::vector<omath::Vector3<std::uint32_t>> k_empty_ebo{};
+
+    constexpr Epa::Params k_default_params{ .max_iterations = 64, .tolerance = 1e-4f };
+
+    Collider make_cube(const Vec3& origin = {}, const Vec3& scale = { 1, 1, 1 })
+    {
+        Mesh m{ k_cube_vbo, k_empty_ebo, scale };
+        m.set_origin(origin);
+        return Collider{ m };
+    }
+
+    // Run GJK then EPA; asserts GJK hit and EPA converged.
+    Epa::Result solve(const Collider& a, const Collider& b,
+                      const Epa::Params& params = k_default_params)
+    {
+        const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+        EXPECT_TRUE(hit) << "GJK must detect collision before EPA can run";
+        auto result = Epa::solve(a, b, simplex, params);
+        EXPECT_TRUE(result.has_value()) << "EPA must converge";
+        return *result;
+    }
+} // namespace
+
+// ---------------------------------------------------------------------------
+// Normal direction per axis
+// ---------------------------------------------------------------------------
+
+// For two unit cubes (half-extent 1) with B offset by d along an axis:
+//   depth = 2 - d   (distance from origin to nearest face of Minkowski diff)
+//   normal component along that axis ≈ ±1
+
+TEST(EpaComprehensive, NormalAlongX_Positive)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0.5f, 0, 0 }));
+    EXPECT_NEAR(std::abs(r.normal.x), 1.f, 1e-3f);
+    EXPECT_NEAR(r.normal.y, 0.f, 1e-3f);
+    EXPECT_NEAR(r.normal.z, 0.f, 1e-3f);
+}
+
+TEST(EpaComprehensive, NormalAlongX_Negative)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ -0.5f, 0, 0 }));
+    EXPECT_NEAR(std::abs(r.normal.x), 1.f, 1e-3f);
+    EXPECT_NEAR(r.normal.y, 0.f, 1e-3f);
+    EXPECT_NEAR(r.normal.z, 0.f, 1e-3f);
+}
+
+TEST(EpaComprehensive, NormalAlongY_Positive)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 0.5f, 0 }));
+    EXPECT_NEAR(r.normal.x, 0.f, 1e-3f);
+    EXPECT_NEAR(std::abs(r.normal.y), 1.f, 1e-3f);
+    EXPECT_NEAR(r.normal.z, 0.f, 1e-3f);
+}
+
+TEST(EpaComprehensive, NormalAlongY_Negative)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, -0.5f, 0 }));
+    EXPECT_NEAR(r.normal.x, 0.f, 1e-3f);
+    EXPECT_NEAR(std::abs(r.normal.y), 1.f, 1e-3f);
+    EXPECT_NEAR(r.normal.z, 0.f, 1e-3f);
+}
+
+TEST(EpaComprehensive, NormalAlongZ_Positive)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, 0.5f }));
+    EXPECT_NEAR(r.normal.x, 0.f, 1e-3f);
+    EXPECT_NEAR(r.normal.y, 0.f, 1e-3f);
+    EXPECT_NEAR(std::abs(r.normal.z), 1.f, 1e-3f);
+}
+
+TEST(EpaComprehensive, NormalAlongZ_Negative)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, -0.5f }));
+    EXPECT_NEAR(r.normal.x, 0.f, 1e-3f);
+    EXPECT_NEAR(r.normal.y, 0.f, 1e-3f);
+    EXPECT_NEAR(std::abs(r.normal.z), 1.f, 1e-3f);
+}
+
+// ---------------------------------------------------------------------------
+// Depth correctness  (depth = 2 - offset for unit cubes)
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, Depth_ShallowOverlap)
+{
+    // offset 1.9 → depth 0.1
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 1.9f, 0, 0 }));
+    EXPECT_NEAR(r.depth, 0.1f, 1e-2f);
+}
+
+TEST(EpaComprehensive, Depth_QuarterOverlap)
+{
+    // offset 1.5 → depth 0.5
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 1.5f, 0, 0 }));
+    EXPECT_NEAR(r.depth, 0.5f, 1e-2f);
+}
+
+TEST(EpaComprehensive, Depth_HalfOverlap)
+{
+    // offset 1.0 → depth 1.0
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 1.0f, 0, 0 }));
+    EXPECT_NEAR(r.depth, 1.0f, 1e-2f);
+}
+
+TEST(EpaComprehensive, Depth_ThreeQuarterOverlap)
+{
+    // offset 0.5 → depth 1.5
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0.5f, 0, 0 }));
+    EXPECT_NEAR(r.depth, 1.5f, 1e-2f);
+}
+
+TEST(EpaComprehensive, Depth_AlongY_HalfOverlap)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 1.0f, 0 }));
+    EXPECT_NEAR(r.depth, 1.0f, 1e-2f);
+}
+
+TEST(EpaComprehensive, Depth_AlongZ_HalfOverlap)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, 1.0f }));
+    EXPECT_NEAR(r.depth, 1.0f, 1e-2f);
+}
+
+// ---------------------------------------------------------------------------
+// Depth monotonicity — deeper overlap → larger depth
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, DepthMonotonic_AlongX)
+{
+    const float d1 = solve(make_cube({ 0, 0, 0 }), make_cube({ 1.9f, 0, 0 })).depth; // ~0.1
+    const float d2 = solve(make_cube({ 0, 0, 0 }), make_cube({ 1.5f, 0, 0 })).depth; // ~0.5
+    const float d3 = solve(make_cube({ 0, 0, 0 }), make_cube({ 1.0f, 0, 0 })).depth; // ~1.0
+    const float d4 = solve(make_cube({ 0, 0, 0 }), make_cube({ 0.5f, 0, 0 })).depth; // ~1.5
+
+    EXPECT_LT(d1, d2);
+    EXPECT_LT(d2, d3);
+    EXPECT_LT(d3, d4);
+}
+
+// ---------------------------------------------------------------------------
+// Normal is a unit vector
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, NormalIsUnit_AlongX)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0.5f, 0, 0 }));
+    EXPECT_NEAR(r.normal.dot(r.normal), 1.f, 1e-5f);
+}
+
+TEST(EpaComprehensive, NormalIsUnit_AlongY)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 1.2f, 0 }));
+    EXPECT_NEAR(r.normal.dot(r.normal), 1.f, 1e-5f);
+}
+
+TEST(EpaComprehensive, NormalIsUnit_AlongZ)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, 0.8f }));
+    EXPECT_NEAR(r.normal.dot(r.normal), 1.f, 1e-5f);
+}
+
+// ---------------------------------------------------------------------------
+// Penetration vector = normal * depth
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, PenetrationVectorLength_EqualsDepth)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0.5f, 0, 0 }));
+    const float pen_len = std::sqrt(r.penetration_vector.dot(r.penetration_vector));
+    EXPECT_NEAR(pen_len, r.depth, 1e-5f);
+}
+
+TEST(EpaComprehensive, PenetrationVectorDirection_ParallelToNormal)
+{
+    const auto r = solve(make_cube({ 0, 0, 0 }), make_cube({ 0, 1.0f, 0 }));
+    // penetration_vector = normal * depth → cross product must be ~zero
+    const auto cross = r.penetration_vector.cross(r.normal);
+    EXPECT_NEAR(cross.dot(cross), 0.f, 1e-8f);
+}
+
+// ---------------------------------------------------------------------------
+// Round-trip: applying penetration_vector separates the shapes
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, RoundTrip_AlongX)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    Mesh mesh_b{ k_cube_vbo, k_empty_ebo };
+    mesh_b.set_origin({ 0.5f, 0, 0 });
+    const auto b = Collider{ mesh_b };
+
+    const auto r = solve(a, b);
+    constexpr float margin = 1.f + 1e-3f;
+
+    // Move B along the penetration vector; it should separate from A
+    Mesh mesh_sep{ k_cube_vbo, k_empty_ebo };
+    mesh_sep.set_origin(mesh_b.get_origin() + r.penetration_vector * margin);
+    EXPECT_FALSE(Gjk::is_collide(a, Collider{ mesh_sep })) << "Applying pen vector must separate";
+
+    // Moving the wrong way must still collide
+    Mesh mesh_wrong{ k_cube_vbo, k_empty_ebo };
+    mesh_wrong.set_origin(mesh_b.get_origin() - r.penetration_vector * margin);
+    EXPECT_TRUE(Gjk::is_collide(a, Collider{ mesh_wrong })) << "Opposite direction must still collide";
+}
+
+TEST(EpaComprehensive, RoundTrip_AlongY)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    Mesh mesh_b{ k_cube_vbo, k_empty_ebo };
+    mesh_b.set_origin({ 0, 0.8f, 0 });
+    const auto b = Collider{ mesh_b };
+
+    const auto r = solve(a, b);
+    constexpr float margin = 1.f + 1e-3f;
+
+    Mesh mesh_sep{ k_cube_vbo, k_empty_ebo };
+    mesh_sep.set_origin(mesh_b.get_origin() + r.penetration_vector * margin);
+    EXPECT_FALSE(Gjk::is_collide(a, Collider{ mesh_sep }));
+
+    Mesh mesh_wrong{ k_cube_vbo, k_empty_ebo };
+    mesh_wrong.set_origin(mesh_b.get_origin() - r.penetration_vector * margin);
+    EXPECT_TRUE(Gjk::is_collide(a, Collider{ mesh_wrong }));
+}
+
+TEST(EpaComprehensive, RoundTrip_AlongZ)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    Mesh mesh_b{ k_cube_vbo, k_empty_ebo };
+    mesh_b.set_origin({ 0, 0, 1.2f });
+    const auto b = Collider{ mesh_b };
+
+    const auto r = solve(a, b);
+    constexpr float margin = 1.f + 1e-3f;
+
+    Mesh mesh_sep{ k_cube_vbo, k_empty_ebo };
+    mesh_sep.set_origin(mesh_b.get_origin() + r.penetration_vector * margin);
+    EXPECT_FALSE(Gjk::is_collide(a, Collider{ mesh_sep }));
+}
+
+// ---------------------------------------------------------------------------
+// Symmetry — swapping A and B preserves depth
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, Symmetry_DepthIsIndependentOfOrder)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+
+    const float depth_ab = solve(a, b).depth;
+    const float depth_ba = solve(b, a).depth;
+
+    EXPECT_NEAR(depth_ab, depth_ba, 1e-2f);
+}
+
+TEST(EpaComprehensive, Symmetry_NormalsAreOpposite)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+
+    const Vec3 n_ab = solve(a, b).normal;
+    const Vec3 n_ba = solve(b, a).normal;
+
+    // The normals should be anti-parallel: n_ab · n_ba ≈ -1
+    EXPECT_NEAR(n_ab.dot(n_ba), -1.f, 1e-3f);
+}
+
+// ---------------------------------------------------------------------------
+// Asymmetric sizes
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, LargeVsSmall_DepthCorrect)
+{
+    // Big (half-ext 2) at origin, small (half-ext 0.5) at (2.0, 0, 0)
+    // Minkowski diff closest face in X at distance 0.5
+    const auto r = solve(make_cube({ 0, 0, 0 }, { 2, 2, 2 }), make_cube({ 2.0f, 0, 0 }, { 0.5f, 0.5f, 0.5f }));
+    EXPECT_NEAR(r.depth, 0.5f, 1e-2f);
+    EXPECT_NEAR(std::abs(r.normal.x), 1.f, 1e-3f);
+}
+
+TEST(EpaComprehensive, LargeVsSmall_RoundTrip)
+{
+    const auto a = make_cube({ 0, 0, 0 }, { 2, 2, 2 });
+
+    Mesh mesh_b{ k_cube_vbo, k_empty_ebo, { 0.5f, 0.5f, 0.5f } };
+    mesh_b.set_origin({ 2.0f, 0, 0 });
+    const auto b = Collider{ mesh_b };
+
+    const auto r = solve(a, b);
+    constexpr float margin = 1.f + 1e-3f;
+
+    Mesh mesh_sep{ k_cube_vbo, k_empty_ebo, { 0.5f, 0.5f, 0.5f } };
+    mesh_sep.set_origin(mesh_b.get_origin() + r.penetration_vector * margin);
+    EXPECT_FALSE(Gjk::is_collide(a, Collider{ mesh_sep }));
+}
+
+// ---------------------------------------------------------------------------
+// Memory resource variants
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, MonotonicBuffer_ConvergesCorrectly)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    ASSERT_TRUE(hit);
+
+    constexpr std::size_t k_buf = 32768;
+    alignas(std::max_align_t) char buf[k_buf];
+    std::pmr::monotonic_buffer_resource mr{ buf, k_buf, std::pmr::null_memory_resource() };
+
+    const auto r = Epa::solve(a, b, simplex, k_default_params, mr);
+    ASSERT_TRUE(r.has_value());
+    EXPECT_NEAR(r->depth, 1.5f, 1e-2f);
+}
+
+TEST(EpaComprehensive, MonotonicBuffer_MultipleReleaseCycles)
+{
+    // Verify mr.release() correctly resets the buffer across multiple calls
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    ASSERT_TRUE(hit);
+
+    constexpr std::size_t k_buf = 32768;
+    alignas(std::max_align_t) char buf[k_buf];
+    std::pmr::monotonic_buffer_resource mr{ buf, k_buf, std::pmr::null_memory_resource() };
+
+    float first_depth = 0.f;
+    for (int i = 0; i < 5; ++i)
+    {
+        mr.release();
+        const auto r = Epa::solve(a, b, simplex, k_default_params, mr);
+        ASSERT_TRUE(r.has_value()) << "solve must converge on iteration " << i;
+        if (i == 0)
+            first_depth = r->depth;
+        else
+            EXPECT_NEAR(r->depth, first_depth, 1e-6f) << "depth must be deterministic";
+    }
+}
+
+TEST(EpaComprehensive, DefaultResource_ConvergesCorrectly)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 1.0f, 0, 0 });
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    ASSERT_TRUE(hit);
+
+    const auto r = Epa::solve(a, b, simplex);
+    ASSERT_TRUE(r.has_value());
+    EXPECT_NEAR(r->depth, 1.0f, 1e-2f);
+}
+
+// ---------------------------------------------------------------------------
+// Tolerance sensitivity
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, TighterTolerance_MoreAccurateDepth)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 1.0f, 0, 0 });
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    ASSERT_TRUE(hit);
+
+    const Epa::Params loose{ .max_iterations = 64, .tolerance = 1e-2f };
+    const Epa::Params tight{ .max_iterations = 64, .tolerance = 1e-5f };
+
+    const auto r_loose = Epa::solve(a, b, simplex, loose);
+    const auto r_tight = Epa::solve(a, b, simplex, tight);
+    ASSERT_TRUE(r_loose.has_value());
+    ASSERT_TRUE(r_tight.has_value());
+
+    // Tighter tolerance must yield a result at least as accurate
+    EXPECT_LE(std::abs(r_tight->depth - 1.0f), std::abs(r_loose->depth - 1.0f) + 1e-4f);
+}
+
+// ---------------------------------------------------------------------------
+// Bookkeeping fields
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, Bookkeeping_IterationsInBounds)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+    const auto r = solve(a, b);
+
+    EXPECT_GT(r.iterations, 0);
+    EXPECT_LE(r.iterations, k_default_params.max_iterations);
+}
+
+TEST(EpaComprehensive, Bookkeeping_FacesAndVerticesGrow)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+    const auto r = solve(a, b);
+
+    // Started with a tetrahedron (4 faces, 4 vertices); EPA must have expanded it
+    EXPECT_GE(r.num_faces, 4);
+    EXPECT_GE(r.num_vertices, 4);
+}
+
+TEST(EpaComprehensive, Bookkeeping_MaxIterationsRespected)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.5f, 0, 0 });
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    ASSERT_TRUE(hit);
+
+    constexpr Epa::Params tight{ .max_iterations = 3, .tolerance = 1e-10f };
+    const auto r = Epa::solve(a, b, simplex, tight);
+
+    // Must return something (fallback best-face path) and respect the cap
+    if (r.has_value())
+        EXPECT_LE(r->iterations, tight.max_iterations);
+}
+
+// ---------------------------------------------------------------------------
+// Determinism
+// ---------------------------------------------------------------------------
+
+TEST(EpaComprehensive, Deterministic_SameResultOnRepeatedCalls)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 0.7f, 0, 0 });
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(a, b);
+    ASSERT_TRUE(hit);
+
+    const auto first = Epa::solve(a, b, simplex);
+    ASSERT_TRUE(first.has_value());
+
+    for (int i = 0; i < 5; ++i)
+    {
+        const auto r = Epa::solve(a, b, simplex);
+        ASSERT_TRUE(r.has_value());
+        EXPECT_NEAR(r->depth, first->depth, 1e-6f);
+        EXPECT_NEAR(r->normal.x, first->normal.x, 1e-6f);
+        EXPECT_NEAR(r->normal.y, first->normal.y, 1e-6f);
+        EXPECT_NEAR(r->normal.z, first->normal.z, 1e-6f);
+    }
+}

tests/general/unit_test_gjk_comprehensive.cpp

@@ -0,0 +1,277 @@
+//
+// Comprehensive GJK tests.
+// Covers: all 6 axis directions, diagonal cases, boundary touching,
+// asymmetric sizes, nesting, symmetry, simplex info, far separation.
+//
+#include <gtest/gtest.h>
+#include <omath/collision/gjk_algorithm.hpp>
+#include <omath/engines/source_engine/collider.hpp>
+#include <omath/engines/source_engine/mesh.hpp>
+
+using Mesh     = omath::source_engine::Mesh;
+using Collider = omath::source_engine::MeshCollider;
+using Gjk      = omath::collision::GjkAlgorithm<Collider>;
+using Vec3     = omath::Vector3<float>;
+
+namespace
+{
+    // Unit cube [-1, 1]^3 in local space.
+    const std::vector<omath::primitives::Vertex<>> k_cube_vbo = {
+        { { -1.f, -1.f, -1.f }, {}, {} },
+        { { -1.f, -1.f,  1.f }, {}, {} },
+        { { -1.f,  1.f, -1.f }, {}, {} },
+        { { -1.f,  1.f,  1.f }, {}, {} },
+        { {  1.f,  1.f,  1.f }, {}, {} },
+        { {  1.f,  1.f, -1.f }, {}, {} },
+        { {  1.f, -1.f,  1.f }, {}, {} },
+        { {  1.f, -1.f, -1.f }, {}, {} },
+    };
+    const std::vector<omath::Vector3<std::uint32_t>> k_empty_ebo{};
+
+    Collider make_cube(const Vec3& origin = {}, const Vec3& scale = { 1, 1, 1 })
+    {
+        Mesh m{ k_cube_vbo, k_empty_ebo, scale };
+        m.set_origin(origin);
+        return Collider{ m };
+    }
+} // namespace
+
+// ---------------------------------------------------------------------------
+// Separation — expect false
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, Separated_AlongPosX)
+{
+    // A extends to x=1, B starts at x=1.1 → clear gap
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 2.1f, 0, 0 })));
+}
+
+TEST(GjkComprehensive, Separated_AlongNegX)
+{
+    // B to the left of A
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ -2.1f, 0, 0 })));
+}
+
+TEST(GjkComprehensive, Separated_AlongPosY)
+{
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 0, 2.1f, 0 })));
+}
+
+TEST(GjkComprehensive, Separated_AlongNegY)
+{
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 0, -2.1f, 0 })));
+}
+
+TEST(GjkComprehensive, Separated_AlongPosZ)
+{
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, 2.1f })));
+}
+
+TEST(GjkComprehensive, Separated_AlongNegZ)
+{
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, -2.1f })));
+}
+
+TEST(GjkComprehensive, Separated_AlongDiagonal)
+{
+    // All components exceed 2.0 — no overlap on any axis
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 2.1f, 2.1f, 2.1f })));
+}
+
+TEST(GjkComprehensive, Separated_LargeDistance)
+{
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 100.f, 0, 0 })));
+}
+
+TEST(GjkComprehensive, Separated_AsymmetricSizes)
+{
+    // Big (scale 2, half-ext 2), small (scale 0.5, half-ext 0.5) at 2.6 → gap of 0.1
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }, { 2, 2, 2 }), make_cube({ 2.6f, 0, 0 }, { 0.5f, 0.5f, 0.5f })));
+}
+
+// ---------------------------------------------------------------------------
+// Overlap — expect true
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, Overlapping_AlongPosX)
+{
+    // B offset 1.5 → overlap depth 0.5 in X
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 1.5f, 0, 0 })));
+}
+
+TEST(GjkComprehensive, Overlapping_AlongNegX)
+{
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ -1.5f, 0, 0 })));
+}
+
+TEST(GjkComprehensive, Overlapping_AlongPosZ)
+{
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, 1.5f })));
+}
+
+TEST(GjkComprehensive, Overlapping_AlongNegZ)
+{
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, -1.5f })));
+}
+
+TEST(GjkComprehensive, Overlapping_AlongDiagonalXY)
+{
+    // Minkowski sum extends ±2 on each axis; offset (1,1,0) is inside
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 1.f, 1.f, 0.f })));
+}
+
+TEST(GjkComprehensive, Overlapping_AlongDiagonalXYZ)
+{
+    // All three axes overlap: (1,1,1) is inside the Minkowski sum
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 1.f, 1.f, 1.f })));
+}
+
+TEST(GjkComprehensive, FullyNested_SmallInsideBig)
+{
+    // Small cube (half-ext 0.5) fully inside big cube (half-ext 2)
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }, { 2, 2, 2 }), make_cube({ 0, 0, 0 }, { 0.5f, 0.5f, 0.5f })));
+}
+
+TEST(GjkComprehensive, FullyNested_OffCenter)
+{
+    // Small at (0.5, 0, 0) still fully inside big (half-ext 2)
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }, { 2, 2, 2 }), make_cube({ 0.5f, 0, 0 }, { 0.5f, 0.5f, 0.5f })));
+}
+
+TEST(GjkComprehensive, Overlapping_AsymmetricSizes)
+{
+    // Big (scale 2, half-ext 2) and small (scale 0.5, half-ext 0.5) at 2.0 → overlap 0.5 in X
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }, { 2, 2, 2 }), make_cube({ 2.0f, 0, 0 }, { 0.5f, 0.5f, 0.5f })));
+}
+
+// ---------------------------------------------------------------------------
+// Boundary cases
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, BoundaryCase_JustColliding)
+{
+    // B at 1.999 — 0.001 overlap in X
+    EXPECT_TRUE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 1.999f, 0, 0 })));
+}
+
+TEST(GjkComprehensive, BoundaryCase_JustSeparated)
+{
+    // B at 2.001 — 0.001 gap in X
+    EXPECT_FALSE(Gjk::is_collide(make_cube({ 0, 0, 0 }), make_cube({ 2.001f, 0, 0 })));
+}
+
+// ---------------------------------------------------------------------------
+// Symmetry
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, Symmetry_WhenColliding)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 1.5f, 0, 0 });
+    EXPECT_EQ(Gjk::is_collide(a, b), Gjk::is_collide(b, a));
+}
+
+TEST(GjkComprehensive, Symmetry_WhenSeparated)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 2.1f, 0.5f, 0 });
+    EXPECT_EQ(Gjk::is_collide(a, b), Gjk::is_collide(b, a));
+}
+
+TEST(GjkComprehensive, Symmetry_DiagonalSeparation)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 1.5f, 1.5f, 1.5f });
+    EXPECT_EQ(Gjk::is_collide(a, b), Gjk::is_collide(b, a));
+}
+
+// ---------------------------------------------------------------------------
+// Simplex info
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, SimplexInfo_HitProducesSimplex4)
+{
+    // On collision the simplex must be a full tetrahedron (4 points)
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(make_cube({ 0, 0, 0 }), make_cube({ 0.5f, 0, 0 }));
+    EXPECT_TRUE(hit);
+    EXPECT_EQ(simplex.size(), 4u);
+}
+
+TEST(GjkComprehensive, SimplexInfo_MissProducesLessThan4)
+{
+    // On non-collision the simplex can never be a full tetrahedron
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(make_cube({ 0, 0, 0 }), make_cube({ 2.1f, 0, 0 }));
+    EXPECT_FALSE(hit);
+    EXPECT_LT(simplex.size(), 4u);
+}
+
+TEST(GjkComprehensive, SimplexInfo_HitAlongY)
+{
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(make_cube({ 0, 0, 0 }), make_cube({ 0, 1.5f, 0 }));
+    EXPECT_TRUE(hit);
+    EXPECT_EQ(simplex.size(), 4u);
+}
+
+TEST(GjkComprehensive, SimplexInfo_HitAlongZ)
+{
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(make_cube({ 0, 0, 0 }), make_cube({ 0, 0, 1.5f }));
+    EXPECT_TRUE(hit);
+    EXPECT_EQ(simplex.size(), 4u);
+}
+
+TEST(GjkComprehensive, SimplexInfo_MissAlongDiagonal)
+{
+    const auto [hit, simplex] = Gjk::is_collide_with_simplex_info(make_cube({ 0, 0, 0 }), make_cube({ 2.1f, 2.1f, 2.1f }));
+    EXPECT_FALSE(hit);
+    EXPECT_LT(simplex.size(), 4u);
+}
+
+// ---------------------------------------------------------------------------
+// Non-trivial geometry — tetrahedron shaped colliders
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, TetrahedronShapes_Overlapping)
+{
+    // A rough tetrahedron mesh; two of them close enough to overlap
+    const std::vector<omath::primitives::Vertex<>> tet_vbo = {
+        { {  0.f,  1.f,  0.f }, {}, {} },
+        { { -1.f, -1.f,  1.f }, {}, {} },
+        { {  1.f, -1.f,  1.f }, {}, {} },
+        { {  0.f, -1.f, -1.f }, {}, {} },
+    };
+
+    Mesh m_a{ tet_vbo, k_empty_ebo };
+    Mesh m_b{ tet_vbo, k_empty_ebo };
+    m_b.set_origin({ 0.5f, 0.f, 0.f });
+
+    EXPECT_TRUE(Gjk::is_collide(Collider{ m_a }, Collider{ m_b }));
+}
+
+TEST(GjkComprehensive, TetrahedronShapes_Separated)
+{
+    const std::vector<omath::primitives::Vertex<>> tet_vbo = {
+        { {  0.f,  1.f,  0.f }, {}, {} },
+        { { -1.f, -1.f,  1.f }, {}, {} },
+        { {  1.f, -1.f,  1.f }, {}, {} },
+        { {  0.f, -1.f, -1.f }, {}, {} },
+    };
+
+    Mesh m_a{ tet_vbo, k_empty_ebo };
+    Mesh m_b{ tet_vbo, k_empty_ebo };
+    m_b.set_origin({ 3.f, 0.f, 0.f });
+
+    EXPECT_FALSE(Gjk::is_collide(Collider{ m_a }, Collider{ m_b }));
+}
+
+// ---------------------------------------------------------------------------
+// Determinism
+// ---------------------------------------------------------------------------
+
+TEST(GjkComprehensive, Deterministic_SameResultOnRepeatedCalls)
+{
+    const auto a = make_cube({ 0, 0, 0 });
+    const auto b = make_cube({ 1.2f, 0.3f, 0.1f });
+    const bool first = Gjk::is_collide(a, b);
+    for (int i = 0; i < 10; ++i)
+        EXPECT_EQ(Gjk::is_collide(a, b), first);
+}

tests/general/unit_test_quaternion.cpp

@@ -0,0 +1,402 @@
+//
+// Created by vlad on 3/1/2026.
+//
+#include <omath/linear_algebra/quaternion.hpp>
+#include <cmath>
+#include <gtest/gtest.h>
+#include <numbers>
+
+using namespace omath;
+
+static constexpr float kEps = 1e-5f;
+
+// ── Helpers ──────────────────────────────────────────────────────────────────
+
+static void expect_quat_near(const Quaternion<float>& a, const Quaternion<float>& b, float eps = kEps)
+{
+    EXPECT_NEAR(a.x, b.x, eps);
+    EXPECT_NEAR(a.y, b.y, eps);
+    EXPECT_NEAR(a.z, b.z, eps);
+    EXPECT_NEAR(a.w, b.w, eps);
+}
+
+static void expect_vec3_near(const Vector3<float>& a, const Vector3<float>& b, float eps = kEps)
+{
+    EXPECT_NEAR(a.x, b.x, eps);
+    EXPECT_NEAR(a.y, b.y, eps);
+    EXPECT_NEAR(a.z, b.z, eps);
+}
+
+// ── Constructors ─────────────────────────────────────────────────────────────
+
+TEST(Quaternion, DefaultConstructorIsIdentity)
+{
+    constexpr Quaternion<float> q;
+    EXPECT_FLOAT_EQ(q.x, 0.f);
+    EXPECT_FLOAT_EQ(q.y, 0.f);
+    EXPECT_FLOAT_EQ(q.z, 0.f);
+    EXPECT_FLOAT_EQ(q.w, 1.f);
+}
+
+TEST(Quaternion, ValueConstructor)
+{
+    constexpr Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    EXPECT_FLOAT_EQ(q.x, 1.f);
+    EXPECT_FLOAT_EQ(q.y, 2.f);
+    EXPECT_FLOAT_EQ(q.z, 3.f);
+    EXPECT_FLOAT_EQ(q.w, 4.f);
+}
+
+TEST(Quaternion, DoubleInstantiation)
+{
+    constexpr Quaternion<double> q{0.0, 0.0, 0.0, 1.0};
+    EXPECT_DOUBLE_EQ(q.w, 1.0);
+}
+
+// ── Equality ─────────────────────────────────────────────────────────────────
+
+TEST(Quaternion, EqualityOperators)
+{
+    constexpr Quaternion<float> a{1.f, 2.f, 3.f, 4.f};
+    constexpr Quaternion<float> b{1.f, 2.f, 3.f, 4.f};
+    constexpr Quaternion<float> c{1.f, 2.f, 3.f, 5.f};
+
+    EXPECT_TRUE(a == b);
+    EXPECT_FALSE(a == c);
+    EXPECT_FALSE(a != b);
+    EXPECT_TRUE(a != c);
+}
+
+// ── Arithmetic ───────────────────────────────────────────────────────────────
+
+TEST(Quaternion, ScalarMultiply)
+{
+    constexpr Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    constexpr auto r = q * 2.f;
+    EXPECT_FLOAT_EQ(r.x, 2.f);
+    EXPECT_FLOAT_EQ(r.y, 4.f);
+    EXPECT_FLOAT_EQ(r.z, 6.f);
+    EXPECT_FLOAT_EQ(r.w, 8.f);
+}
+
+TEST(Quaternion, ScalarMultiplyAssign)
+{
+    Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    q *= 3.f;
+    EXPECT_FLOAT_EQ(q.x, 3.f);
+    EXPECT_FLOAT_EQ(q.y, 6.f);
+    EXPECT_FLOAT_EQ(q.z, 9.f);
+    EXPECT_FLOAT_EQ(q.w, 12.f);
+}
+
+TEST(Quaternion, Addition)
+{
+    constexpr Quaternion<float> a{1.f, 2.f, 3.f, 4.f};
+    constexpr Quaternion<float> b{4.f, 3.f, 2.f, 1.f};
+    constexpr auto r = a + b;
+    EXPECT_FLOAT_EQ(r.x, 5.f);
+    EXPECT_FLOAT_EQ(r.y, 5.f);
+    EXPECT_FLOAT_EQ(r.z, 5.f);
+    EXPECT_FLOAT_EQ(r.w, 5.f);
+}
+
+TEST(Quaternion, AdditionAssign)
+{
+    Quaternion<float> a{1.f, 0.f, 0.f, 0.f};
+    const Quaternion<float> b{0.f, 1.f, 0.f, 0.f};
+    a += b;
+    EXPECT_FLOAT_EQ(a.x, 1.f);
+    EXPECT_FLOAT_EQ(a.y, 1.f);
+}
+
+TEST(Quaternion, UnaryNegation)
+{
+    constexpr Quaternion<float> q{1.f, -2.f, 3.f, -4.f};
+    constexpr auto r = -q;
+    EXPECT_FLOAT_EQ(r.x, -1.f);
+    EXPECT_FLOAT_EQ(r.y,  2.f);
+    EXPECT_FLOAT_EQ(r.z, -3.f);
+    EXPECT_FLOAT_EQ(r.w,  4.f);
+}
+
+// ── Hamilton product ──────────────────────────────────────────────────────────
+
+TEST(Quaternion, MultiplyByIdentityIsNoop)
+{
+    constexpr Quaternion<float> identity;
+    constexpr Quaternion<float> q{0.5f, 0.5f, 0.5f, 0.5f};
+    expect_quat_near(q * identity, q);
+    expect_quat_near(identity * q, q);
+}
+
+TEST(Quaternion, MultiplyAssign)
+{
+    constexpr Quaternion<float> identity;
+    Quaternion<float> q{0.5f, 0.5f, 0.5f, 0.5f};
+    q *= identity;
+    expect_quat_near(q, {0.5f, 0.5f, 0.5f, 0.5f});
+}
+
+TEST(Quaternion, MultiplyKnownResult)
+{
+    // i * j = k  →  (1,0,0,0) * (0,1,0,0) = (0,0,1,0)
+    constexpr Quaternion<float> i{1.f, 0.f, 0.f, 0.f};
+    constexpr Quaternion<float> j{0.f, 1.f, 0.f, 0.f};
+    constexpr auto k = i * j;
+    EXPECT_FLOAT_EQ(k.x, 0.f);
+    EXPECT_FLOAT_EQ(k.y, 0.f);
+    EXPECT_FLOAT_EQ(k.z, 1.f);
+    EXPECT_FLOAT_EQ(k.w, 0.f);
+}
+
+TEST(Quaternion, MultiplyByInverseGivesIdentity)
+{
+    const Quaternion<float> q = Quaternion<float>::from_axis_angle({0.f, 0.f, 1.f},
+                                                                     std::numbers::pi_v<float> / 3.f);
+    const auto result = q * q.inverse();
+    expect_quat_near(result, Quaternion<float>{});
+}
+
+// ── Conjugate ────────────────────────────────────────────────────────────────
+
+TEST(Quaternion, Conjugate)
+{
+    constexpr Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    constexpr auto c = q.conjugate();
+    EXPECT_FLOAT_EQ(c.x, -1.f);
+    EXPECT_FLOAT_EQ(c.y, -2.f);
+    EXPECT_FLOAT_EQ(c.z, -3.f);
+    EXPECT_FLOAT_EQ(c.w,  4.f);
+}
+
+TEST(Quaternion, ConjugateOfIdentityIsIdentity)
+{
+    constexpr Quaternion<float> id;
+    constexpr auto c = id.conjugate();
+    EXPECT_FLOAT_EQ(c.x, 0.f);
+    EXPECT_FLOAT_EQ(c.y, 0.f);
+    EXPECT_FLOAT_EQ(c.z, 0.f);
+    EXPECT_FLOAT_EQ(c.w, 1.f);
+}
+
+// ── Dot / length ─────────────────────────────────────────────────────────────
+
+TEST(Quaternion, Dot)
+{
+    constexpr Quaternion<float> a{1.f, 0.f, 0.f, 0.f};
+    constexpr Quaternion<float> b{0.f, 1.f, 0.f, 0.f};
+    EXPECT_FLOAT_EQ(a.dot(b), 0.f);
+    EXPECT_FLOAT_EQ(a.dot(a), 1.f);
+}
+
+TEST(Quaternion, LengthSqrIdentity)
+{
+    constexpr Quaternion<float> id;
+    EXPECT_FLOAT_EQ(id.length_sqr(), 1.f);
+}
+
+TEST(Quaternion, LengthSqrGeneral)
+{
+    constexpr Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    EXPECT_FLOAT_EQ(q.length_sqr(), 30.f);
+}
+
+TEST(Quaternion, LengthIdentity)
+{
+    const Quaternion<float> id;
+    EXPECT_NEAR(id.length(), 1.f, kEps);
+}
+
+TEST(Quaternion, Normalized)
+{
+    const Quaternion<float> q{1.f, 1.f, 1.f, 1.f};
+    const auto n = q.normalized();
+    EXPECT_NEAR(n.length(), 1.f, kEps);
+    EXPECT_NEAR(n.x, 0.5f, kEps);
+    EXPECT_NEAR(n.y, 0.5f, kEps);
+    EXPECT_NEAR(n.z, 0.5f, kEps);
+    EXPECT_NEAR(n.w, 0.5f, kEps);
+}
+
+TEST(Quaternion, NormalizedOfZeroLengthReturnsSelf)
+{
+    // length_sqr = 0 would be UB, but zero-vector part + zero w is degenerate;
+    // we just verify the guard branch (divides by zero) doesn't crash by
+    // keeping length > 0 via the default constructor path.
+    const Quaternion<float> unit;
+    const auto n = unit.normalized();
+    expect_quat_near(n, unit);
+}
+
+// ── Inverse ───────────────────────────────────────────────────────────────────
+
+TEST(Quaternion, InverseOfUnitIsConjugate)
+{
+    const Quaternion<float> q = Quaternion<float>::from_axis_angle({1.f, 0.f, 0.f},
+                                                                     std::numbers::pi_v<float> / 4.f);
+    const auto inv = q.inverse();
+    const auto conj = q.conjugate();
+    expect_quat_near(inv, conj);
+}
+
+// ── from_axis_angle ──────────────────────────────────────────────────────────
+
+TEST(Quaternion, FromAxisAngleZeroAngleIsIdentity)
+{
+    const auto q = Quaternion<float>::from_axis_angle({1.f, 0.f, 0.f}, 0.f);
+    EXPECT_NEAR(q.x, 0.f, kEps);
+    EXPECT_NEAR(q.y, 0.f, kEps);
+    EXPECT_NEAR(q.z, 0.f, kEps);
+    EXPECT_NEAR(q.w, 1.f, kEps);
+}
+
+TEST(Quaternion, FromAxisAngle90DegZ)
+{
+    const float half_pi = std::numbers::pi_v<float> / 2.f;
+    const auto q = Quaternion<float>::from_axis_angle({0.f, 0.f, 1.f}, half_pi);
+    const float s = std::sin(half_pi / 2.f);
+    const float c = std::cos(half_pi / 2.f);
+    EXPECT_NEAR(q.x, 0.f, kEps);
+    EXPECT_NEAR(q.y, 0.f, kEps);
+    EXPECT_NEAR(q.z, s,   kEps);
+    EXPECT_NEAR(q.w, c,   kEps);
+}
+
+// ── rotate ───────────────────────────────────────────────────────────────────
+
+TEST(Quaternion, RotateByIdentityIsNoop)
+{
+    constexpr Quaternion<float> id;
+    constexpr Vector3<float> v{1.f, 2.f, 3.f};
+    const auto r = id.rotate(v);
+    expect_vec3_near(r, v);
+}
+
+TEST(Quaternion, Rotate90DegAroundZ)
+{
+    // Rotating (1,0,0) by 90° around Z should give (0,1,0)
+    const auto q = Quaternion<float>::from_axis_angle({0.f, 0.f, 1.f}, std::numbers::pi_v<float> / 2.f);
+    const auto r = q.rotate({1.f, 0.f, 0.f});
+    expect_vec3_near(r, {0.f, 1.f, 0.f});
+}
+
+TEST(Quaternion, Rotate180DegAroundY)
+{
+    // Rotating (1,0,0) by 180° around Y should give (-1,0,0)
+    const auto q = Quaternion<float>::from_axis_angle({0.f, 1.f, 0.f}, std::numbers::pi_v<float>);
+    const auto r = q.rotate({1.f, 0.f, 0.f});
+    expect_vec3_near(r, {-1.f, 0.f, 0.f});
+}
+
+TEST(Quaternion, Rotate90DegAroundX)
+{
+    // Rotating (0,1,0) by 90° around X should give (0,0,1)
+    const auto q = Quaternion<float>::from_axis_angle({1.f, 0.f, 0.f}, std::numbers::pi_v<float> / 2.f);
+    const auto r = q.rotate({0.f, 1.f, 0.f});
+    expect_vec3_near(r, {0.f, 0.f, 1.f});
+}
+
+// ── to_rotation_matrix3 ───────────────────────────────────────────────────────
+
+TEST(Quaternion, RotationMatrix3FromIdentityIsIdentityMatrix)
+{
+    constexpr Quaternion<float> id;
+    constexpr auto m = id.to_rotation_matrix3();
+    for (size_t i = 0; i < 3; ++i)
+        for (size_t j = 0; j < 3; ++j)
+            EXPECT_NEAR(m.at(i, j), i == j ? 1.f : 0.f, kEps);
+}
+
+TEST(Quaternion, RotationMatrix3From90DegZ)
+{
+    //  Expected:  | 0 -1  0 |
+    //             | 1  0  0 |
+    //             | 0  0  1 |
+    const auto q = Quaternion<float>::from_axis_angle({0.f, 0.f, 1.f}, std::numbers::pi_v<float> / 2.f);
+    const auto m = q.to_rotation_matrix3();
+    EXPECT_NEAR(m.at(0, 0),  0.f, kEps);
+    EXPECT_NEAR(m.at(0, 1), -1.f, kEps);
+    EXPECT_NEAR(m.at(0, 2),  0.f, kEps);
+    EXPECT_NEAR(m.at(1, 0),  1.f, kEps);
+    EXPECT_NEAR(m.at(1, 1),  0.f, kEps);
+    EXPECT_NEAR(m.at(1, 2),  0.f, kEps);
+    EXPECT_NEAR(m.at(2, 0),  0.f, kEps);
+    EXPECT_NEAR(m.at(2, 1),  0.f, kEps);
+    EXPECT_NEAR(m.at(2, 2),  1.f, kEps);
+}
+
+TEST(Quaternion, RotationMatrix3ConsistentWithRotate)
+{
+    // Matrix-vector multiply must agree with the rotate() method
+    const auto q = Quaternion<float>::from_axis_angle({1.f, 1.f, 0.f}, std::numbers::pi_v<float> / 3.f);
+    const Vector3<float> v{2.f, -1.f, 0.5f};
+
+    const auto rotated = q.rotate(v);
+    const auto m = q.to_rotation_matrix3();
+
+    // manual mat-vec multiply (row-major)
+    const float rx = m.at(0, 0) * v.x + m.at(0, 1) * v.y + m.at(0, 2) * v.z;
+    const float ry = m.at(1, 0) * v.x + m.at(1, 1) * v.y + m.at(1, 2) * v.z;
+    const float rz = m.at(2, 0) * v.x + m.at(2, 1) * v.y + m.at(2, 2) * v.z;
+
+    EXPECT_NEAR(rotated.x, rx, kEps);
+    EXPECT_NEAR(rotated.y, ry, kEps);
+    EXPECT_NEAR(rotated.z, rz, kEps);
+}
+
+// ── to_rotation_matrix4 ───────────────────────────────────────────────────────
+
+TEST(Quaternion, RotationMatrix4FromIdentityIsIdentityMatrix)
+{
+    constexpr Quaternion<float> id;
+    constexpr auto m = id.to_rotation_matrix4();
+    for (size_t i = 0; i < 4; ++i)
+        for (size_t j = 0; j < 4; ++j)
+            EXPECT_NEAR(m.at(i, j), i == j ? 1.f : 0.f, kEps);
+}
+
+TEST(Quaternion, RotationMatrix4HomogeneousRowAndColumn)
+{
+    const auto q = Quaternion<float>::from_axis_angle({1.f, 0.f, 0.f}, std::numbers::pi_v<float> / 5.f);
+    const auto m = q.to_rotation_matrix4();
+
+    // Last row and last column must be (0,0,0,1)
+    for (size_t i = 0; i < 3; ++i)
+    {
+        EXPECT_NEAR(m.at(3, i), 0.f, kEps);
+        EXPECT_NEAR(m.at(i, 3), 0.f, kEps);
+    }
+    EXPECT_NEAR(m.at(3, 3), 1.f, kEps);
+}
+
+TEST(Quaternion, RotationMatrix4Upper3x3MatchesMatrix3)
+{
+    const auto q = Quaternion<float>::from_axis_angle({0.f, 1.f, 0.f}, std::numbers::pi_v<float> / 7.f);
+    const auto m3 = q.to_rotation_matrix3();
+    const auto m4 = q.to_rotation_matrix4();
+
+    for (size_t i = 0; i < 3; ++i)
+        for (size_t j = 0; j < 3; ++j)
+            EXPECT_NEAR(m4.at(i, j), m3.at(i, j), kEps);
+}
+
+// ── as_array ──────────────────────────────────────────────────────────────────
+
+TEST(Quaternion, AsArray)
+{
+    constexpr Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    constexpr auto arr = q.as_array();
+    EXPECT_FLOAT_EQ(arr[0], 1.f);
+    EXPECT_FLOAT_EQ(arr[1], 2.f);
+    EXPECT_FLOAT_EQ(arr[2], 3.f);
+    EXPECT_FLOAT_EQ(arr[3], 4.f);
+}
+
+// ── std::formatter ────────────────────────────────────────────────────────────
+
+TEST(Quaternion, Formatter)
+{
+    const Quaternion<float> q{1.f, 2.f, 3.f, 4.f};
+    const auto s = std::format("{}", q);
+    EXPECT_EQ(s, "[1, 2, 3, 4]");
+}

vcpkg-configuration.json

@@ -1,7 +1,7 @@
 {
   "default-registry": {
     "kind": "git",
-    "baseline": "05442024c3fda64320bd25d2251cc9807b84fb6f",
+    "baseline": "b1b19307e2d2ec1eefbdb7ea069de7d4bcd31f01",
     "repository": "https://github.com/microsoft/vcpkg"
   },
   "registries": [

vcpkg.json

@@ -17,16 +17,9 @@
   ],
   "features": {
     "avx2": {
-      "description": "omath will use AVX2 to boost performance",
+      "description": "Omath will use AVX2 to boost performance",
       "supports": "!arm"
     },
-    "vmprotect": {
-      "description": "omath will use vmprotect sdk to protect sensitive parts of code from reverse engineering",
-      "supports": "windows | linux | osx | android",
-      "dependencies": [
-        "orange-vmprotect-sdk"
-      ]
-    },
     "benchmark": {
       "description": "Build benchmarks",
       "dependencies": [
@@ -42,7 +35,7 @@
       ]
     },
     "imgui": {
-      "description": "omath will define method to convert omath types to imgui types",
+      "description": "Omath will define method to convert omath types to imgui types",
       "dependencies": [
         "imgui"
       ]