added std namspace to int64_t type

include/omath/collision/epa_algorithm.hpp

@@ -8,6 +8,7 @@
 #include <memory>
 #include <memory_resource>
 #include <queue>
+#include <unordered_map>
 #include <utility>
 #include <vector>
 
@@ -50,7 +51,6 @@ namespace omath::collision
             int max_iterations{64};
             FloatingType tolerance{1e-4}; // absolute tolerance on distance growth
         };
-
         // Precondition: simplex.size()==4 and contains the origin.
         [[nodiscard]]
         static std::optional<Result> solve(const ColliderInterfaceType& a, const ColliderInterfaceType& b,
@@ -65,20 +65,14 @@ namespace omath::collision
 
             Result out{};
 
-            // Hoisted outside the loop to reuse the allocation across iterations.
-            std::pmr::vector<Edge> boundary{&mem_resource};
-            boundary.reserve(16);
+            // 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)
             {
                 // Lazily discard stale (deleted or index-mismatched) heap entries.
-                while (!heap.empty())
-                {
-                    const auto& top = heap.top();
-                    if (!faces[top.idx].deleted && faces[top.idx].d == top.d)
-                        break;
-                    heap.pop();
-                }
+                discard_stale_heap_entries(faces, heap);
 
                 if (heap.empty())
                     break;
@@ -105,21 +99,11 @@ namespace omath::collision
 
                 // Tombstone visible faces and collect the horizon boundary.
                 // This avoids copying the faces array (O(n)) each iteration.
-                boundary.clear();
-                for (auto& f : faces)
-                {
-                    if (!f.deleted && visible_from(f, p))
-                    {
-                        f.deleted = true;
-                        add_edge_boundary(boundary, f.i0, f.i1);
-                        add_edge_boundary(boundary, f.i1, f.i2);
-                        add_edge_boundary(boundary, f.i2, f.i0);
-                    }
-                }
+                tombstone_visible_faces(faces, boundary, p);
 
                 // Stitch new faces around the horizon and push them directly onto the
                 // heap — no full O(n log n) rebuild needed.
-                for (const auto& e : boundary)
+                for (const auto& [key, e] : boundary)
                 {
                     const int fi = static_cast<int>(faces.size());
                     faces.emplace_back(make_face(vertexes, e.a, e.b, new_idx));
@@ -133,10 +117,7 @@ namespace omath::collision
             }
 
             // Find the best surviving (non-deleted) face.
-            const Face* best = nullptr;
-            for (const auto& f : faces)
-                if (!f.deleted && (best == nullptr || f.d < best->d))
-                    best = &f;
+            const Face* best = find_best_surviving_face(faces);
 
             if (!best)
                 return std::nullopt;
@@ -153,8 +134,8 @@ namespace omath::collision
         struct Face final
         {
             int i0, i1, i2;
-            VectorType n;       // unit outward normal
-            FloatingType d;     // n · v0  (>= 0 ideally because origin is inside)
+            VectorType n; // unit outward normal
+            FloatingType d; // n · v0  (>= 0 ideally because origin is inside)
             bool deleted{false}; // tombstone flag — avoids O(n) compaction per iteration
         };
 
@@ -180,6 +161,16 @@ 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)
         {
@@ -198,14 +189,16 @@ namespace omath::collision
             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; cancel if opposite already present.
-            auto itb = std::ranges::find_if(boundary, [&](const Edge& e) { return e.a == b && e.b == a; });
-            if (itb != boundary.end())
-                boundary.erase(itb);
+            // O(1) cancel: if the opposite edge (b→a) is already in the map it is an
+            // internal edge shared by two visible faces and must be removed.
+            // Otherwise this is a horizon edge and we insert it.
+            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);
+                boundary.emplace(pack_edge(a, b), Edge{a, b});
         }
 
         [[nodiscard]]
@@ -277,5 +270,41 @@ namespace omath::collision
                 vertexes.emplace_back(simplex[i]);
             return vertexes;
         }
+
+        static const Face* find_best_surviving_face(const std::pmr::vector<Face>& faces)
+        {
+            const Face* best = nullptr;
+            for (const auto& f : faces)
+                if (!f.deleted && (best == nullptr || f.d < best->d))
+                    best = &f;
+            return best;
+        }
+        static void tombstone_visible_faces(std::pmr::vector<Face>& faces, BoundaryMap& boundary,
+                                            const VectorType& p)
+        {
+            boundary.clear();
+            for (auto& f : faces)
+            {
+                if (!f.deleted && visible_from(f, p))
+                {
+                    f.deleted = true;
+                    add_edge_boundary(boundary, f.i0, f.i1);
+                    add_edge_boundary(boundary, f.i1, f.i2);
+                    add_edge_boundary(boundary, f.i2, f.i0);
+                }
+            }
+        }
+
+        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/mesh_collider.hpp

@@ -42,6 +42,16 @@ 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

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);
+    }
+}