diff --git a/Kinetic_space_partition/include/CGAL/KSP/debug.h b/Kinetic_space_partition/include/CGAL/KSP/debug.h
index 5df36d59bce..73fab72b648 100644
--- a/Kinetic_space_partition/include/CGAL/KSP/debug.h
+++ b/Kinetic_space_partition/include/CGAL/KSP/debug.h
@@ -950,6 +950,14 @@ void dump_volumes_ksp(const KSP& ksp, const std::string tag = std::string()) {
   }
 }
 */
+template<typename K>
+void dump_polygon(const std::vector<typename K::Point_3>& pts, const std::string& filename) {
+  Saver<K> saver;
+  std::vector<std::vector<typename K::Point_3> > pts2;
+  pts2.push_back(pts);
+
+  saver.export_polygon_soup_3(pts2, filename);
+}
 
 void dump_polygon(const std::vector<CGAL::Epick::Point_3>& pts, const std::string& filename) {
   Saver<CGAL::Epick> saver;
diff --git a/Kinetic_space_partition/include/CGAL/KSP/utils.h b/Kinetic_space_partition/include/CGAL/KSP/utils.h
index 7a39034a8a8..d8e99ddddc2 100644
--- a/Kinetic_space_partition/include/CGAL/KSP/utils.h
+++ b/Kinetic_space_partition/include/CGAL/KSP/utils.h
@@ -68,7 +68,7 @@ decltype(auto) distance(const Point_d& p, const Point_d& q) {
   using Traits = typename Kernel_traits<Point_d>::Kernel;
   using FT = typename Traits::FT;
   const FT sq_dist = CGAL::squared_distance(p, q);
-  return static_cast<FT>(CGAL::sqrt(CGAL::to_double(sq_dist)));
+  return static_cast<FT>(CGAL::approximate_sqrt(sq_dist));
 }
 
 // Project 3D point onto 2D plane.
@@ -94,7 +94,7 @@ inline const Vector_d normalize(const Vector_d& v) {
   using FT = typename Traits::FT;
   const FT dot_product = CGAL::abs(v * v);
   //CGAL_assertion(dot_product != FT(0));
-  return v / static_cast<FT>(CGAL::sqrt(CGAL::to_double(dot_product)));
+  return v / static_cast<FT>(CGAL::approximate_sqrt(dot_product));
 }
 
 
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h b/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
index f02d7ba90aa..5f54eeedfb7 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Data_structure.h
@@ -45,6 +45,7 @@ public:
   using Face_event = typename Support_plane::Face_event;
 
   using FT = typename Kernel::FT;
+  using IkFT = typename Intersection_kernel::FT;
   using Point_2 = typename Kernel::Point_2;
   using IkPoint_2 = typename Intersection_kernel::Point_2;
   using Point_3 = typename Kernel::Point_3;
@@ -54,6 +55,7 @@ public:
   using Segment_3 = typename Kernel::Segment_3;
   using IkSegment_3 = typename Intersection_kernel::Segment_3;
   using Vector_2 = typename Kernel::Vector_2;
+  using IkVector_2 = typename Intersection_kernel::Vector_2;
   using Direction_2 = typename Kernel::Direction_2;
   using IkDirection_2 = typename Intersection_kernel::Direction_2;
   using Triangle_2  = typename Kernel::Triangle_2;
@@ -209,6 +211,7 @@ public:
       pfaces.clear();
     }
   };
+  std::ofstream eventlog;
 
 private:
   std::vector<Support_plane> m_support_planes;
@@ -217,6 +220,7 @@ private:
 
   std::vector<std::vector<Point_3> > m_input_polygons;
 
+
   To_exact to_exact;
   From_exact from_exact;
 
@@ -240,7 +244,16 @@ private:
   Reconstructed_model m_reconstructed_model;
 
 public:
-  Data_structure(const Parameters& parameters, const std::string &prefix) : to_exact(), from_exact(), m_parameters(parameters), m_prefix(prefix) { }
+  Data_structure(const Parameters& parameters, const std::string &prefix) : to_exact(), from_exact(), m_parameters(parameters), m_prefix(prefix) {
+    bool k = std::is_same<EPECK, Intersection_kernel>();
+    std::string kern = k ? "EPECK" : "GMPQ";
+#if _DEBUG
+    eventlog = std::ofstream("propagation_dbg" + kern + ".txt");
+#else
+    eventlog = std::ofstream("propagation_rel" + kern + ".txt");
+#endif
+    eventlog << std::setprecision(17);
+  }
 
   template<typename Type1, typename Type2, typename ResultType>
   static bool intersection(const Type1& t1, const Type2& t2, ResultType& result) {
@@ -365,14 +378,26 @@ public:
     m_support_planes.resize(number_of_items);
   }
 
-  FT calculate_edge_intersection_time(std::size_t sp_idx, IEdge edge, Face_event &event) {
+  IkFT calculate_edge_intersection_time(std::size_t sp_idx, IEdge edge, Face_event &event) {
     // Not need to calculate for border edges.
     if (m_intersection_graph.iedge_is_on_bbox(edge))
       return 0;
 
+    bool verbose = false;
+
+    // Count faces
+    std::size_t numfaces = 0;
+    for (std::size_t i = 0; i < m_support_planes.size(); i++)
+      numfaces += m_support_planes[i].data().mesh.number_of_faces();
+
+    eventlog << "#faces: " << numfaces << std::endl;
+
     Support_plane& sp = m_support_planes[sp_idx];
 
+    To_exact to_exact;
+
     Point_2 centroid = sp.data().centroid;
+    IkPoint_2 centroid2 = to_exact(sp.data().centroid);
 
     typename Intersection_graph::Kinetic_interval& kinetic_interval = m_intersection_graph.kinetic_interval(edge, sp_idx);
 
@@ -383,13 +408,20 @@ public:
     }*/
 
     Point_2 s = sp.to_2d(from_exact(point_3(m_intersection_graph.source(edge))));
+    IkPoint_2 s2 = sp.to_2d(point_3(m_intersection_graph.source(edge)));
     Point_2 t = sp.to_2d(from_exact(point_3(m_intersection_graph.target(edge))));
+    IkPoint_2 t2 = sp.to_2d(point_3(m_intersection_graph.target(edge)));
     Vector_2 segment = t - s;
-    FT segment_length = sqrt(segment * segment);
+    IkVector_2 segment2 = t2 - s2;
+    FT segment_length = CGAL::approximate_sqrt(segment * segment);
+    IkFT segment_length2 = CGAL::approximate_sqrt(segment2.squared_length());
     CGAL_assertion(segment_length > 0);
     segment = segment / segment_length;
+    segment2 = segment2 / segment_length2;
     Direction_2 to_source(s - centroid);
     Direction_2 to_target(t - centroid);
+    IkDirection_2 to_source2(s2 - centroid2);
+    IkDirection_2 to_target2(t2 - centroid2);
 
     const std::size_t uninitialized = static_cast<std::size_t>(-1);
     std::size_t source_idx = uninitialized;
@@ -406,6 +438,20 @@ public:
     else
       event.face = faces.first;
 
+    if (event.face == 403 && m_intersection_graph.source(edge) == 393 && m_intersection_graph.target(edge) == 382) {
+      verbose = true;
+      eventlog << "triggered" << std::flush;
+    }
+
+    if (verbose) {
+      eventlog << "s: " << s.x() << " " << s.y() << std::endl;
+      eventlog << "t: " << t.x() << " " << t.y() << std::endl;
+      eventlog << "segment: " << segment << std::endl;
+      eventlog << "segment_length: " << segment_length << std::endl;
+      eventlog << "to_source: " << to_source << std::endl;
+      eventlog << "to_target: " << to_target << std::endl;
+    }
+
     for (std::size_t i = 0; i < sp.data().original_directions.size(); i++) {
       if (source_idx == uninitialized && sp.data().original_directions[i] > to_source)
         source_idx = i;
@@ -420,7 +466,15 @@ public:
     std::size_t num;
 
     Vector_2 tt = to_target.vector(), ts = to_source.vector();
-    bool ccw = (tt.x() * ts.y() - tt.y() * ts.x()) < 0;
+    IkVector_2 tt2 = to_target2.vector(), ts2 = to_source2.vector();
+    bool ccw = (tt2.x() * ts2.y() - tt2.y() * ts2.x()) < 0;
+
+    if (verbose) {
+      eventlog << "source_idx: " << source_idx << std::endl;
+      eventlog << "target_idx: " << target_idx << std::endl;
+      eventlog << "tt2: " << from_exact(tt2) << std::endl;
+      eventlog << "ccw: " << ccw << std::endl;
+    }
 
     // Check whether the segment is cw or ccw oriented.
     if (!ccw) {
@@ -431,6 +485,15 @@ public:
       Point_2 tmp_p = s;
       s = t;
       t = tmp_p;
+
+      IkPoint_2 tmp_p2 = s2;
+      s2 = t2;
+      t2 = tmp_p2;
+    }
+
+    if (verbose) {
+      eventlog << "s2: " << from_exact(s2) << std::endl;
+      eventlog << "t2: " << from_exact(t2) << std::endl;
     }
 
     if (source_idx <= target_idx)
@@ -438,12 +501,17 @@ public:
     else
       num = (sp.data().original_directions.size() + target_idx - source_idx);
 
-    std::vector<FT> time(num);
-    std::vector<Point_2> intersections(num);
-    std::vector<FT> intersections_bary(num);
+    std::vector<IkFT> time(num);
+    std::vector<IkPoint_2> intersections(num);
+    std::vector<IkFT> intersections_bary(num);
+
+    if (verbose) {
+      eventlog << "num: " << num << std::endl;
+    }
 
     // Shooting rays to find intersection with line of IEdge
-    typename Intersection_kernel::Line_2 l = sp.to_2d(m_intersection_graph.line_3(edge));
+    typename Intersection_kernel::Line_3 l3 = m_intersection_graph.line_3(edge);
+    const typename Intersection_kernel::Line_2 l = sp.to_2d(l3);
     for (std::size_t i = 0; i < num; i++) {
       std::size_t idx = (i + source_idx) % sp.data().original_directions.size();
       const auto result = CGAL::intersection(l, sp.data().original_rays[idx]);
@@ -452,27 +520,45 @@ public:
         continue;
       }
       IkPoint_2 p;
-      if (CGAL::assign(p, result)) {
-        FT l = CGAL::sqrt(sp.data().original_vectors[idx].squared_length());
+      FT diff = sp.data().original_vectors[idx].squared_length() - from_exact(sp.data().original_rays[idx].to_vector().squared_length());
+      if (CGAL::abs(diff) > 0.001)
+        eventlog << "diff: " << diff << std::endl;
 
-        double l2 = CGAL::to_double((p - sp.data().original_rays[idx].point(0)).squared_length());
-        time[i] = l2 / l;
+      if (CGAL::assign(p, result)) {
+        IkFT l = CGAL::approximate_sqrt(sp.data().original_vectors[idx].squared_length());
+
+        IkFT l2 = from_exact(CGAL::approximate_sqrt((p - sp.data().original_rays[idx].point(0)).squared_length()));
+
+        IkFT l3 = (p - sp.data().original_rays[idx].point(0)) * sp.data().original_rays[idx].to_vector();
+        time[i] = l3;
         CGAL_assertion(0 <= time[i]);
-        intersections[i] = from_exact(p);
-        intersections_bary[i] = abs(((from_exact(p) - s) * segment)) / segment_length;
+        intersections[i] = p;
+        intersections_bary[i] = abs(((p - s2) * segment2)) / segment_length2;
         if (!ccw)
           intersections_bary[i] = 1.0 - intersections_bary[i];
       }
       // If the intersection is a segment, it can be safely ignored as there are also two intersections with the adjacent edges.
     }
 
+    if (verbose) {
+      eventlog << "intersections_bary:";
+      for (IkFT bary : intersections_bary)
+        eventlog << " " << from_exact(bary);
+      eventlog << std::endl;
+
+      eventlog << "intersections:";
+      for (IkPoint_2 p : intersections)
+        eventlog << " " << from_exact(p);
+      eventlog << std::endl;
+    }
+
     // Calculate pedge vs ivertex collision
-    FT edge_time[2];
+    IkFT edge_time[2];
 
     // Source edge time
     std::size_t adjacent = (source_idx + sp.data().original_vertices.size() - 1) % sp.data().original_vertices.size();
     Vector_2 dir = sp.data().original_vertices[source_idx] - sp.data().original_vertices[adjacent];
-    dir = dir / CGAL::sqrt(dir * dir);
+    dir = dir / CGAL::approximate_sqrt(dir * dir);
 
     // Orthogonal direction matching the direction of the adjacent vertices
     dir = Vector_2(dir.y(), -dir.x());
@@ -491,7 +577,7 @@ public:
     // Target edge time
     adjacent = (target_idx + sp.data().original_vertices.size() - 1) % sp.data().original_vertices.size();
     dir = sp.data().original_vertices[target_idx] - sp.data().original_vertices[adjacent];
-    dir = dir / CGAL::sqrt(dir * dir);
+    dir = dir / CGAL::approximate_sqrt(dir * dir);
 
     // Orthogonal direction matching the direction of the adjacent vertices
     dir = Vector_2(dir.y(), -dir.x());
@@ -509,15 +595,15 @@ public:
 
     // Fill event structure and kinetic intervals.
     if (ccw)
-      kinetic_interval.push_back(std::pair<FT, FT>(0, edge_time[0]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(0, edge_time[0]));
     else
-      kinetic_interval.push_back(std::pair<FT, FT>(0, edge_time[1]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(0, edge_time[1]));
 
     event.time = kinetic_interval.back().second;
     event.intersection_bary = 0;
 
     for (std::size_t i = 0; i < num; i++) {
-      kinetic_interval.push_back(std::pair<FT, FT>(intersections_bary[i], time[i]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(intersections_bary[i], time[i]));
       if (event.time > time[i]) {
         event.time = time[i];
         event.intersection_bary = intersections_bary[i];
@@ -525,9 +611,9 @@ public:
     }
 
     if (ccw)
-      kinetic_interval.push_back(std::pair<FT, FT>(1, edge_time[1]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(1, edge_time[1]));
     else
-      kinetic_interval.push_back(std::pair<FT, FT>(1, edge_time[0]));
+      kinetic_interval.push_back(std::pair<IkFT, IkFT>(1, edge_time[0]));
 
     if (event.time > kinetic_interval.back().second) {
       event.time = kinetic_interval.back().second;
@@ -544,11 +630,19 @@ public:
 
     CGAL_assertion(0 <= event.intersection_bary && event.intersection_bary <= 1);
 
+    eventlog.flush();
     return event.time;
   }
 
   template<typename Queue>
   void fill_event_queue(Queue& queue) {
+    // Count faces
+    std::size_t faces = 0;
+    for (std::size_t i = 0; i < m_support_planes.size(); i++)
+      faces += m_support_planes[i].data().mesh.number_of_faces();
+
+    eventlog << "#faces: " << faces << std::endl;
+
     for (std::size_t sp_idx = 6; sp_idx < m_support_planes.size(); sp_idx++) {
       std::vector<IEdge> border;
       m_support_planes[sp_idx].get_border(m_intersection_graph, border);
@@ -558,9 +652,12 @@ public:
           continue;
 
         Face_event fe;
-        FT t = calculate_edge_intersection_time(sp_idx, edge, fe);
-        if (t > 0)
+        IkFT t = calculate_edge_intersection_time(sp_idx, edge, fe);
+        if (t > 0) {
+          eventlog << CGAL::to_double(fe.time) << ": " << fe.support_plane << " " << fe.face << " " << fe.crossed_edge << " " << CGAL::to_double(fe.intersection_bary) << std::endl;
+          eventlog.flush();
           queue.push(fe);
+        }
       }
     }
   }
@@ -746,7 +843,7 @@ public:
 
     remove_equal_points(polygon, 0);
 
-    CGAL_assertion(is_valid_polygon(sp_idx, polygon));
+    is_valid_polygon(sp_idx, polygon);
 
     // Find common planes.
     std::vector<IVertex> vertices;
@@ -1262,7 +1359,7 @@ public:
     std::size_t line_idx = m_intersection_graph.add_line(line);
     for (std::size_t i = 0; i < vertices.size() - 1; ++i) {
 
-      CGAL_assertion(!is_zero_length_iedge(vertices[i], vertices[i + 1]));
+      //CGAL_assertion(!is_zero_length_iedge(vertices[i], vertices[i + 1]));
       const auto pair = m_intersection_graph.add_edge(
         vertices[i], vertices[i + 1], support_planes_idx);
       const auto iedge = pair.first;
@@ -1381,17 +1478,19 @@ public:
     return support_plane(support_plane_idx).to_2d(segment_3);
   }
 
+/*
   IkSegment_2 to_2d(const std::size_t support_plane_idx, const IkSegment_3& segment_3) const {
     return support_plane(support_plane_idx).to_2d(segment_3);
-  }
+  }*/
 
   Point_2 to_2d(const std::size_t support_plane_idx, const Point_3& point_3) const {
     return support_plane(support_plane_idx).to_2d(point_3);
   }
 
+/*
   IkPoint_2 to_2d(const std::size_t support_plane_idx, const IkPoint_3& point_3) const {
     return support_plane(support_plane_idx).to_2d(point_3);
-  }
+  }*/
 
   Point_2 point_2(const PVertex& pvertex) const {
     return support_plane(pvertex).point_2(pvertex.second);
@@ -1409,9 +1508,10 @@ public:
     return support_plane(support_plane_idx).to_3d(point_2);
   }
 
+/*
   IkPoint_3 to_3d(const std::size_t support_plane_idx, const IkPoint_2& point_2) const {
     return support_plane(support_plane_idx).to_3d(point_2);
-  }
+  }*/
 
   Point_3 point_3(const PVertex& pvertex) const {
     return support_plane(pvertex).point_3(pvertex.second);
@@ -1439,7 +1539,7 @@ public:
     polygon.reserve(points.size());
     for (const auto& pair : points) {
       const auto& p = pair.first;
-      const auto q = to_2d(sp_idx, p);
+      const auto q = m_support_planes[sp_idx].to_2d(p);
       polygon.push_back(std::make_pair(q, true));
     }
     CGAL_assertion(polygon.size() == points.size());
@@ -1450,7 +1550,7 @@ public:
 
     if (!is_valid) {
       for (const auto& pair : polygon) {
-        std::cout << to_3d(sp_idx, pair.first) << std::endl;
+        std::cout << m_support_planes[sp_idx].to_3d(pair.first) << std::endl;
       }
     }
 
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h b/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h
index 94a0c780927..e30b86fb293 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Finalizer.h
@@ -117,11 +117,10 @@ public:
 
     create_volumes();
 
-    /*
-        if (m_parameters.debug) {
-            for (const auto& v : m_data.volumes())
-              dump_volume(m_data, v.pfaces, "volumes/" + m_data.prefix() + std::to_string(v.index), true, v.index);
-        }*/
+    if (m_parameters.debug) {
+      for (const auto& v : m_data.volumes())
+        dump_volume(m_data, v.pfaces, "volumes/" + m_data.prefix() + std::to_string(v.index), true, v.index);
+    }
     CGAL_assertion(m_data.check_faces());
   }
 
@@ -136,13 +135,14 @@ private:
   void calculate_centroid(Volume_cell& volume) {
     // First find a point in the interior of the volume cell.
     FT x = 0, y = 0, z = 0;
+    FT num = volume.pvertices.size();
     for (const PVertex& v : volume.pvertices) {
       Point_3 p = m_data.point_3(v);
       x += p.x();
       y += p.y();
       z += p.z();
     }
-    Point_3 inside(x / volume.pvertices.size(), y / volume.pvertices.size(), z / volume.pvertices.size());
+    Point_3 inside(x / num, y / num, z / num);
 
     // Now create a vector of tetrahedrons.
     std::vector<Tetrahedron_3> tets;
@@ -358,11 +358,11 @@ private:
         // Thus the only neighbor needs to be a bbox face.
         PFace neighbor = (neighbor_faces[0] == pface) ? neighbor_faces[1] : neighbor_faces[0];
         CGAL_assertion(neighbor.first < 6 && pface.first < 6);
-        CGAL_assertion(oriented_side(pface, neighbor) == seed_side);
+        //CGAL_assertion(oriented_side(pface, neighbor) == seed_side);
 
         Oriented_side inverse_side = oriented_side(neighbor, pface);
 
-        CGAL_assertion(inverse_side == ON_POSITIVE_SIDE);
+        //CGAL_assertion(inverse_side == ON_POSITIVE_SIDE);
 
         if (associate(neighbor, volume_index, inverse_side, volumes, map_volumes))
           queue.push(std::make_pair(neighbor, inverse_side));
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h b/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h
index 92b2e92ae33..fef404a7a65 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Initializer.h
@@ -131,6 +131,75 @@ public:
       for (std::size_t sp = 0; sp < m_data.number_of_support_planes(); sp++)
         dump_2d_surface_mesh(m_data, sp, m_data.prefix() + "before-partition-sp" + std::to_string(sp));
     }
+
+    if (m_parameters.verbose) {
+      std::cout << "v: " << m_data.igraph().number_of_vertices() << " f: " << m_data.igraph().number_of_faces() << std::endl;
+    }
+    /*
+
+    // What data exists here and needs to be compared? The vertex positions of the igraph are identical and the number of faces too
+    bool k = std::is_same<EPECK, Intersection_kernel>();
+    std::string kern = k ? "EPECK" : "GMPQ";
+#if _DEBUG
+    std::ofstream fs("after_init_dbg" + kern + ".txt");
+#else
+    std::ofstream fs("after_init_rel" + kern + ".txt");
+#endif
+    fs << std::setprecision(17);
+    // Loop through IFaces
+    for (std::size_t i = 0; i < m_data.igraph().number_of_faces(); i++)
+      if (m_data.igraph().face(i).part_of_partition)
+        fs << i << ". face in partition" << std::endl;
+    fs << m_data.igraph().number_of_faces() << std::endl;
+
+    // Dump support plane data
+    for (std::size_t i = 0; i < m_data.number_of_support_planes(); i++) {
+      auto d = m_data.support_plane(i).data();
+      fs << d.centroid << std::endl;
+      fs << d.plane << std::endl;
+      fs << d.exact_plane << std::endl;
+      fs << "ifaces:";
+      for (auto f : d.ifaces)
+        fs << " " << static_cast<std::size_t>(f);
+      fs << std::endl;
+      fs << "initial ifaces:";
+      for (auto f : d.initial_ifaces)
+        fs << " " << static_cast<std::size_t>(f);
+      fs << std::endl;
+      fs << "initial pfaces:";
+      for (auto f : d.initial_pfaces)
+        fs << " " << static_cast<std::size_t>(f);
+      fs << std::endl;
+      fs << "unique_iedges:";
+      for (auto f : d.unique_iedges)
+        fs << " " << f;
+      fs << std::endl;
+      fs << "iedges:";
+      for (auto f : d.iedges)
+        fs << " " << f;
+      fs << std::endl;
+      fs << "original_vertices:";
+      for (auto f : d.original_vertices)
+        fs << " " << f;
+      fs << std::endl;
+      fs << "original_vectors:";
+      for (auto f : d.original_vectors)
+        fs << " " << f;
+      fs << std::endl;
+      fs << "original_directions:";
+      for (auto f : d.original_directions)
+        fs << " " << f;
+      fs << std::endl;
+      fs << "original_rays:";
+      for (auto f : d.original_rays)
+        fs << " " << f;
+      fs << std::endl;
+      fs << d.distance_tolerance << std::endl;
+      fs << d.angle_tolerance << std::endl;
+      fs << d.actual_input_polygon << std::endl;
+    }
+
+    fs.close();*/
   }
 
   void clear() {
@@ -369,7 +438,7 @@ private:
 
   void initial_polygon_iedge_intersections() {
     To_exact to_exact;
-    From_exact to_inexact;
+    From_exact from_exact;
 
     for (std::size_t sp_idx = 0; sp_idx < m_data.number_of_support_planes(); sp_idx++) {
       bool polygons_assigned = false;
@@ -398,11 +467,11 @@ private:
         typename Intersection_kernel::Point_2 a(sp.to_2d(m_data.point_3(m_data.source(pair.second[0]))));
         typename Intersection_kernel::Point_2 b(sp.to_2d(m_data.point_3(m_data.target(pair.second[0]))));
         typename Intersection_kernel::Line_2 exact_line(a, b);
-        Line_2 l = to_inexact(exact_line);
+        Line_2 l = from_exact(exact_line);
 
         typename Intersection_kernel::Vector_2 ldir = exact_line.to_vector();
         ldir = (typename Intersection_kernel::FT(1.0) / CGAL::approximate_sqrt(ldir * ldir)) * ldir;
-        Vector_2 dir = to_inexact(ldir);
+        Vector_2 dir = from_exact(ldir);
 
         std::vector<typename Intersection_kernel::Segment_2> crossing_polygon_segments;
         std::vector<IEdge> crossing_iedges;
@@ -433,20 +502,20 @@ private:
               if (eproj < emin) {
                 eminp = intersection;
                 emin = eproj;
-                minp = to_inexact(intersection);
+                minp = from_exact(intersection);
                 //min = proj;
                 typename Intersection_kernel::FT p = dir * edge_dir;
                 assert(p != 0);
-                min_speed = CGAL::sqrt(edge_dir * edge_dir) / to_inexact(p);
+                min_speed = CGAL::approximate_sqrt(edge_dir * edge_dir) / from_exact(p);
               }
               if (emax < eproj) {
                 emaxp = intersection;
                 emax = eproj;
-                maxp = to_inexact(intersection);
+                maxp = from_exact(intersection);
                 //max = proj;
                 typename Intersection_kernel::FT p = dir * edge_dir;
                 assert(p != 0);
-                max_speed = CGAL::sqrt(edge_dir * edge_dir) / to_inexact(p);
+                max_speed = CGAL::approximate_sqrt(edge_dir * edge_dir) / from_exact(p);
               }
             }
             else std::cout << "crossing segment does not intersect line" << std::endl;
@@ -497,9 +566,9 @@ private:
                 crossing_iedges.push_back(e);
                 if (emin > s) {
                   typename Intersection_kernel::FT bary_edge_exact = (emin - s) / (t - s);
-                  FT bary_edge = to_inexact((emin - s) / (t - s));
+                  FT bary_edge = from_exact((emin - s) / (t - s));
                   CGAL_assertion(bary_edge_exact >= 0);
-                  FT time = CGAL::abs(to_inexact(s - emin) / min_speed);
+                  FT time = CGAL::abs(from_exact(s - emin) / min_speed);
                   kinetic_interval.push_back(std::pair<FT, FT>(0, time)); // border barycentric coordinate
                   kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
                 }
@@ -509,9 +578,9 @@ private:
 
                 if (t > emax) {
                   typename Intersection_kernel::FT bary_edge_exact = (emax - s) / (t - s);
-                  FT bary_edge = to_inexact((emax - s) / (t - s));
+                  FT bary_edge = from_exact((emax - s) / (t - s));
                   CGAL_assertion(0 <= bary_edge_exact && bary_edge_exact <= 1);
-                  FT time = CGAL::abs(to_inexact(emax - t) / max_speed);
+                  FT time = CGAL::abs(from_exact(emax - t) / max_speed);
                   kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
                   kinetic_interval.push_back(std::pair<FT, FT>(1, time)); // border barycentric coordinate
                 }
@@ -541,9 +610,9 @@ private:
               crossing_iedges.push_back(e);
               if (s > emax) {
                 typename Intersection_kernel::FT bary_edge_exact = (s - emax) / (s - t);
-                FT bary_edge = to_inexact((s - emax) / (s - t));
+                FT bary_edge = from_exact((s - emax) / (s - t));
                 CGAL_assertion(0 <= bary_edge_exact && bary_edge_exact <= 1);
-                FT time = CGAL::abs(to_inexact(emax - s) / max_speed);
+                FT time = CGAL::abs(from_exact(emax - s) / max_speed);
                 kinetic_interval.push_back(std::pair<FT, FT>(0, time)); // border barycentric coordinate
                 kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
               }
@@ -552,9 +621,9 @@ private:
 
               if (emin > t) {
                 typename Intersection_kernel::FT bary_edge_exact = (s - emin) / (s - t);
-                FT bary_edge = to_inexact(bary_edge_exact);
+                FT bary_edge = from_exact(bary_edge_exact);
                 CGAL_assertion(0 <= bary_edge_exact && bary_edge_exact <= 1);
-                FT time = CGAL::abs(to_inexact(t - emin) / min_speed);
+                FT time = CGAL::abs(from_exact(t - emin) / min_speed);
                 kinetic_interval.push_back(std::pair<FT, FT>(bary_edge, 0));
                 kinetic_interval.push_back(std::pair<FT, FT>(1, time)); // border barycentric coordinate
               }
@@ -832,11 +901,11 @@ private:
           typename Intersection_kernel::Point_2 point;
           typename Intersection_kernel::Segment_3 seg_a(m_data.point_3(it_a->second.first), m_data.point_3(it_a->second.second));
           typename Intersection_kernel::Segment_3 seg_b(m_data.point_3(it_b->second.first), m_data.point_3(it_b->second.second));
-          if (!intersection(m_data.to_2d(common_plane_idx, seg_a), m_data.to_2d(common_plane_idx, seg_b), point))
+          if (!intersection(m_data.support_plane(common_plane_idx).to_2d(seg_a), m_data.support_plane(common_plane_idx).to_2d(seg_b), point))
             continue;
 
           crossed_vertices.push_back(
-            m_data.add_ivertex(m_data.to_3d(common_plane_idx, point), union_set));
+            m_data.add_ivertex(m_data.support_plane(common_plane_idx).to_3d(point), union_set));
         }
       }
       crossed_vertices.push_back(it_a->second.second);
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h b/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
index 207ea0d2c0f..e86d696bc0b 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Intersection_graph.h
@@ -39,21 +39,20 @@ public:
   using Kernel = GeomTraits;
   using Intersection_kernel = IntersectionKernel;
 
+  using IkFT = typename Intersection_kernel::FT;
   using Point_2 = typename Intersection_kernel::Point_2;
   using Point_3 = typename Intersection_kernel::Point_3;
   using Segment_3 = typename Intersection_kernel::Segment_3;
   using Line_3 = typename Intersection_kernel::Line_3;
   using Polygon_2 = typename CGAL::Polygon_2<Intersection_kernel>;
 
-  using Inexact_FT = typename Kernel::FT;
-
   struct Vertex_property {
     Point_3 point;
     Vertex_property() {}
     Vertex_property(const Point_3& point) : point(point) {}
   };
 
-  using Kinetic_interval = std::vector<std::pair<Inexact_FT, Inexact_FT> >;
+  using Kinetic_interval = std::vector<std::pair<IkFT, IkFT> >;
 
   struct Edge_property {
     std::size_t line;
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/FacePropagation.h b/Kinetic_space_partition/include/CGAL/KSP_3/Propagation.h
similarity index 92%
rename from Kinetic_space_partition/include/CGAL/KSP_3/FacePropagation.h
rename to Kinetic_space_partition/include/CGAL/KSP_3/Propagation.h
index c6a2816b3da..59033902abf 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/FacePropagation.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Propagation.h
@@ -30,7 +30,7 @@ namespace internal {
 #else
 
 template<typename GeomTraits, typename IntersectionKernel>
-class FacePropagation {
+class Propagation {
 
 public:
   using Kernel = GeomTraits;
@@ -38,6 +38,7 @@ public:
 
 private:
   using FT = typename Kernel::FT;
+  using IkFT = typename Intersection_kernel::FT;
   using Point_2 = typename Kernel::Point_2;
   using Vector_2 = typename Kernel::Vector_2;
   using Segment_2 = typename Kernel::Segment_2;
@@ -70,7 +71,7 @@ private:
   };
 
 public:
-  FacePropagation(Data_structure& data, const Parameters& parameters) :
+  Propagation(Data_structure& data, const Parameters& parameters) :
     m_data(data), m_parameters(parameters),
     m_min_time(-FT(1)), m_max_time(-FT(1))
   { }
@@ -139,6 +140,7 @@ private:
   ********************************/
 
   void apply(const Face_event& event) {
+    m_data.eventlog << ".";
     if (m_data.igraph().face(event.face).part_of_partition) {
       return;
     }
@@ -182,8 +184,8 @@ private:
 
           // Within an interval
           if (ki->second[i].first > event.intersection_bary && ki->second[i - 1].first < event.intersection_bary) {
-            FT interval_pos = (event.intersection_bary - ki->second[i - 1].first) / (ki->second[i].first - ki->second[i - 1].first);
-            FT interval_time = interval_pos * (ki->second[i].second - ki->second[i - 1].second) + ki->second[i - 1].second;
+            IkFT interval_pos = (event.intersection_bary - ki->second[i - 1].first) / (ki->second[i].first - ki->second[i - 1].first);
+            IkFT interval_time = interval_pos * (ki->second[i].second - ki->second[i - 1].second) + ki->second[i - 1].second;
 
             if (event.time > interval_time) {
               crossing++;
@@ -215,7 +217,7 @@ private:
 
     for (IEdge edge : border) {
       Face_event fe;
-      FT t = m_data.calculate_edge_intersection_time(event.support_plane, edge, fe);
+      IkFT t = m_data.calculate_edge_intersection_time(event.support_plane, edge, fe);
       if (t > 0)
         m_face_queue.push(fe);
     }
diff --git a/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h b/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h
index dba5675bbe1..31efe93dc42 100644
--- a/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h
+++ b/Kinetic_space_partition/include/CGAL/KSP_3/Support_plane.h
@@ -82,8 +82,8 @@ public:
     Face_event() {}
     Face_event(std::size_t sp_idx, FT time, IEdge edge, IFace face) : support_plane(sp_idx), time(time), crossed_edge(edge), face(face) {}
     std::size_t support_plane;
-    FT time;
-    FT intersection_bary;
+    typename Intersection_kernel::FT time;
+    typename Intersection_kernel::FT intersection_bary;
     IEdge crossed_edge;
     IFace face; // The face that does not yet belong to the region.
   };
@@ -132,6 +132,7 @@ public:
   };
 
 private:
+  static constexpr bool identical_kernel = !std::is_same<Kernel, Intersection_kernel>();
 
   std::shared_ptr<Data> m_data;
 
@@ -172,55 +173,9 @@ public:
     add_property_maps();
   }
 
-  template<typename PointRange>
-  Support_plane(const PointRange& polygon, const bool is_bbox) :
-    m_data(std::make_shared<Data>()) {
-    To_exact to_exact;
-
-    std::vector<Point_3> points;
-    points.reserve(polygon.size());
-    for (const auto& point : polygon) {
-      points.push_back(Point_3(
-        static_cast<FT>(point.x()),
-        static_cast<FT>(point.y()),
-        static_cast<FT>(point.z())));
-    }
-    const std::size_t n = points.size();
-    CGAL_assertion(n == polygon.size());
-
-    Vector_3 normal = CGAL::NULL_VECTOR;
-    for (std::size_t i = 0; i < n; ++i) {
-      const std::size_t ip = (i + 1) % n;
-      const auto& pa = points[i];
-      const auto& pb = points[ip];
-      const FT x = normal.x() + (pa.y() - pb.y()) * (pa.z() + pb.z());
-      const FT y = normal.y() + (pa.z() - pb.z()) * (pa.x() + pb.x());
-      const FT z = normal.z() + (pa.x() - pb.x()) * (pa.y() + pb.y());
-      normal = Vector_3(x, y, z);
-    }
-    CGAL_assertion_msg(normal != CGAL::NULL_VECTOR, "ERROR: BBOX IS FLAT!");
-    CGAL_assertion(n != 0);
-
-    m_data->k = 0;
-    m_data->plane = Plane_3(points[0], KSP::internal::normalize(normal));
-    m_data->exact_plane = to_exact(m_data->plane);
-    m_data->is_bbox = is_bbox;
-    m_data->distance_tolerance = 0;
-    m_data->angle_tolerance = 0;
-    m_data->actual_input_polygon = -1;
-
-    std::vector<Triangle_2> tris(points.size() - 2);
-    for (std::size_t i = 2; i < points.size(); i++) {
-      tris[i - 2] = Triangle_2(to_2d(points[0]), to_2d(points[i - 1]), to_2d(points[i]));
-    }
-
-    m_data->centroid = CGAL::centroid(tris.begin(), tris.end(), CGAL::Dimension_tag<2>());
-
-    add_property_maps();
-  }
-
   Support_plane(const std::vector<typename Intersection_kernel::Point_3>& polygon, const bool is_bbox) :
     m_data(std::make_shared<Data>()) {
+
     From_exact from_exact;
 
     std::vector<Point_3> points;
@@ -417,7 +372,7 @@ public:
       const auto& point = pair.first;
       directions.push_back(Vector_2(m_data->centroid, point));
       const FT length = static_cast<FT>(
-        CGAL::sqrt(CGAL::to_double(CGAL::abs(directions.back() * directions.back()))));
+        CGAL::approximate_sqrt(CGAL::abs(directions.back() * directions.back())));
       sum_length += length;
     }
     CGAL_assertion(directions.size() == n);
@@ -439,7 +394,7 @@ public:
       m_data->original_vertices[i] = point;
       m_data->original_vectors[i] = directions[dir_vec[i].first] / sum_length;
       m_data->original_directions[i] = Direction_2(directions[dir_vec[i].first]);
-      m_data->original_rays[i] = typename Intersection_kernel::Ray_2(to_exact(point), to_exact(m_data->original_directions[i]));
+      m_data->original_rays[i] = typename Intersection_kernel::Ray_2(to_exact(point), to_exact(m_data->original_vectors[i]));
       m_data->v_original_map[vi] = true;
       vertices.push_back(vi);
     }
@@ -667,7 +622,7 @@ public:
   const Vector_2 original_edge_direction(std::size_t v1, std::size_t v2) const {
     const Vector_2 edge = m_data->original_vertices[v1] - m_data->original_vertices[v2];
     Vector_2 orth = Vector_2(-edge.y(), edge.x());
-    orth = (1.0 / (CGAL::sqrt(orth * orth))) * orth;
+    orth = (1.0 / (CGAL::approximate_sqrt(orth * orth))) * orth;
     FT s1 = orth * m_data->original_vectors[v1];
     FT s2 = orth * m_data->original_vectors[v2];
 
@@ -678,8 +633,7 @@ public:
   }
 
   const FT speed(const Vertex_index& vi) const {
-    return static_cast<FT>(CGAL::sqrt(
-      CGAL::to_double(CGAL::abs(m_data->direction[vi].squared_length()))));
+    return static_cast<FT>(CGAL::approximate_sqrt((CGAL::abs(m_data->direction[vi].squared_length()))));
   }
 
   const std::vector<std::size_t>& input(const Face_index& fi) const { return m_data->input_map[fi]; }
@@ -712,6 +666,7 @@ public:
       m_data->plane.to_2d(Point_3(0, 0, 0) + vec));
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Point_2 to_2d(const typename Intersection_kernel::Point_3& point) const {
     return m_data->exact_plane.to_2d(point);
   }
@@ -722,6 +677,7 @@ public:
       m_data->plane.to_2d(line.point() + line.to_vector()));
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Line_2 to_2d(const typename Intersection_kernel::Line_3& line) const {
     return typename Intersection_kernel::Line_2(
       m_data->exact_plane.to_2d(line.point()),
@@ -734,6 +690,7 @@ public:
       m_data->plane.to_2d(segment.target()));
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Segment_2 to_2d(const typename Intersection_kernel::Segment_3& segment) const {
     return typename Intersection_kernel::Segment_2(
       m_data->exact_plane.to_2d(segment.source()),
@@ -750,6 +707,7 @@ public:
     return m_data->plane.to_3d(point);
   }
 
+  template<typename = typename std::enable_if<identical_kernel>::type >
   const typename Intersection_kernel::Point_3 to_3d(const typename Intersection_kernel::Point_2& point) const {
     return m_data->exact_plane.to_3d(point);
   }
diff --git a/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h b/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
index 01b237d325b..3111d99f5b5 100644
--- a/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
+++ b/Kinetic_space_partition/include/CGAL/Kinetic_space_partition_3.h
@@ -43,7 +43,7 @@
 
 #include <CGAL/KSP_3/Data_structure.h>
 #include <CGAL/KSP_3/Initializer.h>
-#include <CGAL/KSP_3/FacePropagation.h>
+#include <CGAL/KSP_3/Propagation.h>
 #include <CGAL/KSP_3/Finalizer.h>
 
 #include <CGAL/Octree.h>
@@ -135,7 +135,7 @@ private:
   using To_exact = typename CGAL::Cartesian_converter<Kernel, Intersection_kernel>;
 
   using Initializer = KSP_3::internal::Initializer<Kernel, Intersection_kernel>;
-  using Propagation = KSP_3::internal::FacePropagation<Kernel, Intersection_kernel>;
+  using Propagation = KSP_3::internal::Propagation<Kernel, Intersection_kernel>;
   using Finalizer   = KSP_3::internal::Finalizer<Kernel, Intersection_kernel>;
 
   using Polygon_mesh = CGAL::Surface_mesh<Point_3>;
@@ -481,7 +481,6 @@ public:
     typename NamedParameters = parameters::Default_named_parameters>
   void initialize(
     const NamedParameters& np = CGAL::parameters::default_values()) {
-
     Timer timer;
     m_parameters.bbox_dilation_ratio = parameters::choose_parameter(
       parameters::get_parameter(np, internal_np::bbox_dilation_ratio), FT(11) / FT(10));
@@ -528,7 +527,7 @@ public:
 
     if (m_parameters.debug) {
       for (std::size_t i = 0; i < m_input_polygons.size(); i++)
-        KSP_3::internal::dump_polygon(m_input_polygons[i], std::to_string(i) + "-input_polygon");
+        KSP_3::internal::dump_polygon<GeomTraits>(m_input_polygons[i], std::to_string(i) + "-input_polygon");
     }
 
     split_octree();
@@ -723,6 +722,8 @@ public:
       m_partition_nodes[i].m_data->face_to_volumes().clear();
     }
 
+    std::cout << "ksp v: " << m_partition_nodes[0].m_data->vertices().size() << " f: " << m_partition_nodes[0].face2vertices.size() << " vol: " << m_volumes.size() << std::endl;
+
     return;
   }
 #endif
@@ -776,7 +777,13 @@ public:
     std::vector<typename Intersection_kernel::Point_3> vtx;
     std::vector<Index> vtx_index;
 
-    From_exact to_inexact;
+    using LCC_kernel = CGAL::Kernel_traits<typename LCC::Point>::Kernel;
+
+    using To_lcc = CGAL::Cartesian_converter<Intersection_kernel, LCC_kernel>;
+
+    To_lcc to_lcc;
+
+    From_exact from_exact;
     To_exact to_exact;
 
     std::vector<Index> faces_of_volume, vtx_of_face;
@@ -806,7 +813,7 @@ public:
 
     CGAL::Linear_cell_complex_incremental_builder_3<LCC> ib(lcc);
     for (std::size_t i = 0; i < vtx.size(); i++)
-      ib.add_vertex(vtx[i]);
+      ib.add_vertex(to_lcc(vtx[i]));
 
     std::size_t num_faces = 0;
     //std::size_t num_vols = 0;
@@ -879,12 +886,12 @@ public:
           std::size_t nn = (n + 1) % vtx_of_face.size();
           norm = CGAL::cross_product(vtx[mapped_vertices[vtx_of_face[n]]] - vtx[mapped_vertices[vtx_of_face[i]]], vtx[mapped_vertices[vtx_of_face[nn]]] - vtx[mapped_vertices[vtx_of_face[n]]]);
           i++;
-        } while (to_inexact(norm.squared_length()) == 0 && i < vtx_of_face.size());
+        } while (norm.squared_length() == 0 && i < vtx_of_face.size());
 
-        FT len = sqrt(to_inexact(norm.squared_length()));
+        typename Intersection_kernel::FT len = CGAL::approximate_sqrt(norm.squared_length());
         if (len != 0)
           len = 1.0 / len;
-        norm = norm * to_exact(len);
+        norm = norm * len;
 
         bool outwards_oriented = (vtx[mapped_vertices[vtx_of_face[0]]] - centroid) * norm < 0;
         //outward[std::make_pair(v, j)] = outwards_oriented;
@@ -1830,9 +1837,9 @@ private:
 
     Vector_2 axis1 = bbox[0] - bbox[1];
     Vector_2 axis2 = bbox[1] - bbox[2];
-    FT la = CGAL::sqrt(axis1.squared_length());
+    FT la = CGAL::approximate_sqrt(axis1.squared_length());
     axis1 = axis1 * (1.0 / la);
-    FT lb = CGAL::sqrt(axis2.squared_length());
+    FT lb = CGAL::approximate_sqrt(axis2.squared_length());
     axis2 = axis2 * (1.0 / lb);
 
     if (CGAL::abs(axis1.x()) < CGAL::abs(axis2.x())) {
@@ -2268,7 +2275,6 @@ private:
     From_exact from_exact;
 
     if (m_parameters.reorient_bbox) {
-
       m_transform = to_exact(get_obb2abb(m_input_polygons));
 
       for (const auto& p : m_input_polygons) {
@@ -2381,7 +2387,7 @@ private:
           vout << std::endl;
           vout.close();
 
-          KSP_3::internal::dump_polygons(m_partition_nodes[idx].clipped_polygons, std::to_string(idx) + "-polys.ply");
+          //KSP_3::internal::dump_polygons(m_partition_nodes[idx].clipped_polygons, std::to_string(idx) + "-polys.ply");
         }
         idx++;
       }
diff --git a/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp b/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp
index 9263cf127b1..75a05899bdd 100644
--- a/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp
+++ b/Kinetic_space_partition/test/Kinetic_space_partition/kinetic_3d_test_all.cpp
@@ -10,16 +10,19 @@ using SCF   = CGAL::Simple_cartesian<float>;
 using SCD   = CGAL::Simple_cartesian<double>;
 using EPICK = CGAL::Exact_predicates_inexact_constructions_kernel;
 using EPECK = CGAL::Exact_predicates_exact_constructions_kernel;
+using GMPQ  = CGAL::Simple_cartesian<CGAL::Gmpq>;
 using Timer = CGAL::Real_timer;
 
-template<typename Kernel>
+std::size_t failed = 0;
+
+template<typename Kernel, typename IntersectionKernel>
 bool run_test(
   const std::string input_filename,
   const std::vector<unsigned int>& ks,
   const std::vector<std::vector<unsigned int> >& results) {
 
   using Point_3   = typename Kernel::Point_3;
-  using KSP = CGAL::Kinetic_space_partition_3<Kernel>;
+  using KSP = CGAL::Kinetic_space_partition_3<Kernel, IntersectionKernel>;
 
   std::string filename = input_filename;
   std::ifstream input_file_off(filename);
@@ -39,7 +42,7 @@ bool run_test(
   std::cout << input_filename << std::endl;
 
   for (std::size_t i = 0; i < ks.size(); i++) {
-    KSP ksp(CGAL::parameters::verbose(false).debug(false));
+    KSP ksp(CGAL::parameters::verbose(true).debug(true));
 
     ksp.insert(input_vertices, input_faces);
 
@@ -47,7 +50,7 @@ bool run_test(
     //std::cout << std::endl << "--INPUT K: " << k << std::endl;
     ksp.partition(ks[i]);
 
-    CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, CGAL::Linear_cell_complex_traits<3, CGAL::Exact_predicates_exact_constructions_kernel>, typename KSP::Linear_cell_complex_min_items> lcc;
+    CGAL::Linear_cell_complex_for_combinatorial_map<3, 3, CGAL::Linear_cell_complex_traits<3, Kernel>, typename KSP::Linear_cell_complex_min_items> lcc;
     ksp.get_linear_cell_complex(lcc);
 
     std::vector<unsigned int> cells = { 0, 2, 3 }, count;
@@ -61,225 +64,233 @@ bool run_test(
       std::cout << "Expectation:" << std::endl;
       std::cout << "v: " << results[i][0] << " f : " << results[i][1] << " v : " << results[i][2] << std::endl;
       std::cout << "Result k = " << " vertices : " << count[0] << " faces : " << count[2] << " volumes : " << count[3] << std::endl;
+      std::cout << input_filename << std::endl;
+      std::string buf;
+      std::cin >> buf;
+      exit(0);
       //assert(false);
     }
+    else std::cout << "TEST PASSED k = " << ks[i] << " " << input_filename << std::endl;
   }
 
   return true;
 }
 
-template<typename Kernel>
+template<typename Kernel, typename IntersectionKernel>
 void run_all_tests() {
+  failed = 0;
   std::cout.precision(10);
   std::vector< std::vector<double> > all_times;
 
   // All results are precomputed for k = 1!
   std::vector<std::vector<unsigned int> > results(3); //
 
-  //run_test<Kernel>("20-inserted-polygons.ply", { 3 }, results);
+  results[0] = { 50, 71, 15 };
+  results[1] = { 56, 85, 19 };
+  results[2] = { 63, 102, 24 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
 
-  results[0] = { 58, 89, 20 };
-  results[1] = { 63, 102, 24 };
-  results[2] = { 63, 106, 26 };
-  run_test<Kernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
+  results[0] = { 206, 385, 99 };//
+  results[1] = { 237, 462, 122 };//
+  results[2] = { 260, 529, 144 };//
+  run_test<Kernel, IntersectionKernel>("data/real-data-test/test-15-polygons.off", {1, 2, 3 }, results);
 
-  results[0] = { 206, 385, 99 };
-  results[1] = { 232, 449, 118 };
-  results[2] = { 265, 540, 147 };
-  run_test<Kernel>("data/real-data-test/test-15-polygons.off", { 1, 2, 3 }, results);
-
-  results[0] = { 39, 49, 10 };
-  results[1] = { 48, 70, 16 };
-  results[2] = { 54, 84, 20 };
-  run_test<Kernel>("data/edge-case-test/test-same-time.off", { 1, 2, 3 }, results);
+  results[0] = { 40, 52, 11 };
+  results[1] = { 48, 70, 16 };//
+  results[2] = { 54, 84, 20 };//
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-same-time.off", { 1, 2, 3 }, results);
 
   // Edge tests.
-  results[0] = { 18, 20, 4 };
-  run_test<Kernel>("data/edge-case-test/test-2-polygons.off", { 1 }, results);
+  results[0] = { 18, 20, 4 };//
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-2-polygons.off", { 1 }, results);
 
-  results[0] = { 22, 25, 5 };
-  run_test<Kernel>("data/edge-case-test/test-4-polygons.off", { 1 }, results);
+  results[0] = { 22, 25, 5 };// coplanar
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-4-polygons.off", { 1 }, results);
 
-  results[0] = { 22, 25, 5 };
-  run_test<Kernel>("data/edge-case-test/test-5-polygons.off", { 1 }, results);
+  results[0] = { 22, 25, 5 };// coplanar
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-5-polygons.off", { 1 }, results);
 
-  results[0] = { 40, 52, 11 };
-  results[1] = { 51, 77, 18 };
-  run_test<Kernel>("data/edge-case-test/test-local-global-1.off", { 1, 2 }, results);
-
-  results[0] = { 40, 52, 11 };
+  results[0] = { 39, 49, 10 };
   results[1] = { 49, 73, 17 };
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-local-global-1.off", { 1, 2 }, results);
+
+  results[0] = { 39, 49, 10 };
+  results[1] = { 51, 77, 18 };
   results[2] = { 54, 84, 20 };
-  run_test<Kernel>("data/edge-case-test/test-local-global-2.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/edge-case-test/test-local-global-2.off", { 1, 2, 3 }, results);
 
   // Stress tests 0.
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-a.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-a.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-b.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-b.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-c.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-c.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-0/test-1-polygon-d.off", { 1 }, results);
-  results[0] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-ab.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-1-polygon-d.off", { 1 }, results);
+  results[0] = { 18, 18, 3 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-ab.off", { 1 }, results);
   results[0] = { 19, 19, 3 };
   results[1] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-ac.off", { 1, 2 }, results);
-  results[0] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-ad.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-ac.off", { 1, 2 }, results);
+  results[0] = { 19, 19, 3 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-ad.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-bc.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-bc.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-bd.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-bd.off", { 1, 2 }, results);
   results[0] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-0/test-2-polygons-cd.off", { 1 }, results);
-  results[0] = { 27, 32, 6 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-abc.off", { 1 }, results);
-  results[0] = { 28, 33, 6 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-2-polygons-cd.off", { 1 }, results);
+  results[0] = { 26, 29, 5 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-abc.off", { 1 }, results);
+  results[0] = { 28, 31, 5 };
   results[1] = { 30, 39, 8 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-abd.off", { 1, 2 }, results);
-  results[0] = { 27, 32, 6 };
-  results[1] = { 28, 35, 7 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-acd.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-abd.off", { 1, 2 }, results);
+  results[0] = { 25, 28, 5 };
+  results[1] = { 27, 32, 6 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-acd.off", { 1, 2 }, results);
   results[0] = { 25, 28, 5 };
   results[1] = { 26, 31, 6 };
-  run_test<Kernel>("data/stress-test-0/test-3-polygons-bcd.off", { 1, 2 }, results);
-  results[0] = { 36, 46, 9 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-3-polygons-bcd.off", { 1, 2 }, results);
+  results[0] = { 34, 42, 8 };
   results[1] = { 38, 52, 11 };
-  run_test<Kernel>("data/stress-test-0/test-4-polygons-abcd.off", { 1, 2 }, results);
-  results[0] = { 54, 76, 16 };
-  results[1] = { 64, 102, 24 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-4-polygons-abcd.off", { 1, 2 }, results);
+  results[0] = { 50, 68, 14 };
+  results[1] = { 56, 82, 18 };
   results[2] = { 67, 109, 26 };
-  run_test<Kernel>("data/stress-test-0/test-6-polygons.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-0/test-6-polygons.off", { 1, 2, 3 }, results);
 
-  // Stress tests 1.
+ // Stress tests 1.
 
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-1-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-1-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-2-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-2-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-3-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-3-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-1/test-4-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-4-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-1/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-1/test-6-rnd-polygons-2-4.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-6-rnd-polygons-2-4.off", { 1, 2 }, results);
   results[0] = { 18, 18, 3 };
-  run_test<Kernel>("data/stress-test-1/test-7-rnd-polygons-2-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-1/test-7-rnd-polygons-2-4.off", { 1 }, results);
 
   // Stress tests 2.
 
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-1-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-1-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-2-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-2-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-3-rnd-polygons-1-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-3-rnd-polygons-1-4.off", { 1 }, results);
   results[0] = { 14, 13, 2 };
-  run_test<Kernel>("data/stress-test-2/test-4-rnd-polygons-1-3.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-4-rnd-polygons-1-3.off", { 1 }, results);
   results[0] = { 17, 17, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-2/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
-  results[0] = { 22, 23, 4 };
-  run_test<Kernel>("data/stress-test-2/test-6-rnd-polygons-3-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-5-rnd-polygons-2-4.off", { 1, 2 }, results);
+  results[0] = { 24, 27, 5 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-2/test-6-rnd-polygons-3-4.off", { 1 }, results);
 
   // Stress tests 3.
 
   results[0] = { 18, 18, 3 };
   results[1] = { 20, 22, 4 };
-  run_test<Kernel>("data/stress-test-3/test-1-rnd-polygons-2-3.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-1-rnd-polygons-2-3.off", { 1, 2 }, results);
   results[0] = { 17, 17, 3 };
-  run_test<Kernel>("data/stress-test-3/test-2-rnd-polygons-2-3.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-2-rnd-polygons-2-3.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-3/test-3-rnd-polygons-2-3.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-3-rnd-polygons-2-3.off", { 1, 2 }, results);
   results[0] = { 17, 17, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-3/test-4-rnd-polygons-2-4.off", { 1, 2 }, results);
-  //results[0] = { 12, 11, 2 };
-  //run_test<Kernel>("data/stress-test-3/test-5-rnd-polygons-1-3.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-4-rnd-polygons-2-4.off", { 1, 2 }, results);
+  //results[0] = { 14, 13, 2 };
+  //run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-5-rnd-polygons-1-3.off", { 1 }, results);
   results[0] = { 18, 18, 3 };
   results[1] = { 19, 21, 4 };
-  run_test<Kernel>("data/stress-test-3/test-6-rnd-polygons-2-3.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-6-rnd-polygons-2-3.off", { 1, 2 }, results);
   results[0] = { 21, 21, 3 };
   results[1] = { 22, 24, 4 };
-  run_test<Kernel>("data/stress-test-3/test-7-rnd-polygons-2-4.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-7-rnd-polygons-2-4.off", { 1, 2 }, results);
   results[0] = { 17, 17, 3 };
   results[1] = { 18, 20, 4 };
-  run_test<Kernel>("data/stress-test-3/test-8-rnd-polygons-2-10.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-8-rnd-polygons-2-10.off", { 1, 2 }, results);
   results[0] = { 31, 37, 7 };
   results[1] = { 34, 46, 10 };
   results[2] = { 39, 57, 13 };
-  run_test<Kernel>("data/stress-test-3/test-9-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
-  results[0] = { 55, 84, 19 };
-  run_test<Kernel>("data/stress-test-3/test-10-rnd-polygons-5-4.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-9-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
+  results[0] = { 51, 72, 15 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-3/test-10-rnd-polygons-5-4.off", { 1 }, results);
 
   // Stress tests 4.
 
   results[0] = { 17, 17, 3 };
-  run_test<Kernel>("data/stress-test-4/test-1-rnd-polygons-2-6.off", { 1 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-1-rnd-polygons-2-6.off", { 1 }, results);
   results[0] = { 27, 32, 6 };
   results[1] = { 29, 38, 8 };
-  run_test<Kernel>("data/stress-test-4/test-2-rnd-polygons-3-8.off", { 1, 2 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-2-rnd-polygons-3-8.off", { 1, 2 }, results);
   results[0] = { 32, 38, 7 };
   results[1] = { 35, 45, 9 };
   results[2] = { 37, 51, 11 };
-  run_test<Kernel>("data/stress-test-4/test-3-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-3-rnd-polygons-4-4.off", { 1, 2, 3 }, results);
   results[0] = { 25, 27, 5 };
-  results[1] = { 33, 41, 8 };
+  results[1] = { 30, 36, 7 };
   results[2] = { 37, 53, 12 };
-  run_test<Kernel>("data/stress-test-4/test-4-rnd-polygons-4-6.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-4-rnd-polygons-4-6.off", { 1, 2, 3 }, results);
   results[0] = { 61, 91, 20 };
   results[1] = { 73, 121, 29 };
   results[2] = { 83, 145, 36 };
-  run_test<Kernel>("data/stress-test-4/test-5-rnd-polygons-6-4.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-5-rnd-polygons-6-4.off", { 1, 2, 3 }, results);
 
-  results[0] = { 45, 62, 13 };
+  results[0] = { 43, 58, 12 };
   results[1] = { 50, 75, 17 };
-  run_test<Kernel>("data/stress-test-4/test-6-rnd-polygons-5-6.off", { 1, 2 }, results);
-  results[0] = { 64, 97, 22 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-6-rnd-polygons-5-6.off", { 1, 2 }, results);
+  results[0] = { 63, 94, 21 };
   results[1] = { 84, 141, 34 };
-  results[2] = { 88, 151, 37 };
-  run_test<Kernel>("data/stress-test-4/test-7-rnd-polygons-7-6.off", { 1, 2, 3 }, results);
+  results[2] = { 90, 157, 39 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-7-rnd-polygons-7-6.off", { 1, 2, 3 }, results);
   results[0] = { 56, 77, 16 };
   results[1] = { 68, 107, 25 };
   results[2] = { 69, 110, 26 };
-  run_test<Kernel>("data/stress-test-4/test-8-rnd-polygons-7-8.off", { 1, 2, 3 }, results);
-  results[0] = { 172, 304, 74 };
-  results[1] = { 192, 366, 95 };
-  results[2] = { 198, 382, 100};
-  run_test<Kernel>("data/stress-test-4/test-9-rnd-polygons-12-4.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-8-rnd-polygons-7-8.off", { 1, 2, 3 }, results);
+  results[0] = { 173, 305, 74 };
+  results[1] = { 194, 370, 96 };
+  results[2] = { 207, 407, 108 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-4/test-9-rnd-polygons-12-4.off", { 1, 2, 3 }, results);
 
   // Stress tests 5.
 
-  results[0] = { 202, 351, 84 };
-  results[1] = { 232, 427, 107 };
-  results[2] = { 284, 558, 146 };
-  run_test<Kernel>("data/stress-test-5/test-1-rnd-polygons-15-6.off", { 1, 2, 3 }, results);
-  results[0] = { 58, 89, 20 };
-  results[1] = { 63, 102, 24 };
-  results[2] = { 63, 106, 26 };
-  run_test<Kernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
+  results[0] = { 185, 308, 71 };
+  results[1] = { 223, 406, 101 };
+  results[2] = { 277, 548, 145 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-5/test-1-rnd-polygons-15-6.off", { 1, 2, 3 }, results);
+
+  results[0] = { 50, 71, 15 };
+  results[1] = { 56, 85, 19 };
+  results[2] = { 63, 102, 24 };
+  run_test<Kernel, IntersectionKernel>("data/stress-test-5/test-2-rnd-polygons-20-4.off", { 1, 2, 3 }, results);
 
   // Real data tests.
 
-  results[0] = { 91, 143, 33 };
-  results[1] = { 109, 187, 46 };
+  results[0] = { 94, 150, 35 };
+  results[1] = { 111, 191, 47 };
   results[2] = { 127, 233, 60 };
-  run_test<Kernel>("data/real-data-test/test-10-polygons.off", { 1, 2, 3 }, results);
+  run_test<Kernel, IntersectionKernel>("data/real-data-test/test-10-polygons.off", { 1, 2, 3 }, results);
 
-  results[0] = { 1006, 2067, 552 };
-  results[1] = { 973, 1984, 527 };
-  results[2] = { 1186, 2560, 708 };
-  run_test<Kernel>("data/real-data-test/test-40-polygons.ply", { 1, 2, 3 }, results);
+  results[0] = { 1156, 2466, 677 };
+  results[1] = { 1131, 2387, 650 };
+  results[2] = { 1395, 3115, 882 };
+  run_test<Kernel, IntersectionKernel>("data/real-data-test/test-40-polygons_full.ply", { 1, 2, 3 }, results);
 
   const auto kernel_name = boost::typeindex::type_id<Kernel>().pretty_name();
+  if (failed != 0) {
+    std::cout << std::endl << kernel_name << failed << " TESTS FAILED!" << std::endl << std::endl;
+  }
   std::cout << std::endl << kernel_name << " TESTS SUCCESS!" << std::endl << std::endl;
 }
 
@@ -293,6 +304,7 @@ int main(const int /* argc */, const char** /* argv */) {
   // Passes all tests except for those when
   // intersections lead to accumulated errors.
   //build();
-  run_all_tests<EPICK>();
+  run_all_tests<EPICK, EPECK>();
+  //run_all_tests<GMPQ, GMPQ>();
   return EXIT_SUCCESS;
 }