diff --git a/Distance_3/benchmark/Distance_3/CMakeLists.txt b/Distance_3/benchmark/Distance_3/CMakeLists.txt
new file mode 100644
index 00000000000..3cc6ca18f7d
--- /dev/null
+++ b/Distance_3/benchmark/Distance_3/CMakeLists.txt
@@ -0,0 +1,23 @@
+# Created by the script cgal_create_cmake_script
+# This is the CMake script for compiling a CGAL application.
+
+cmake_minimum_required(VERSION 3.12...3.31)
+project(Distance_3_Benchs)
+
+find_package(CGAL REQUIRED)
+
+if(MSVC)
+  # Turn off a VC++ warning on a potential division by zero
+  # in Cartesian_kernel/include/CGAL/constructions/kernel_ftC3.h
+  # where CGAL_assume() does not help
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /wd4723")
+endif()
+
+# create a target per cppfile
+file(
+  GLOB cppfiles
+  RELATIVE ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp)
+foreach(cppfile ${cppfiles})
+  create_single_source_cgal_program("${cppfile}")
+endforeach()
diff --git a/Distance_3/benchmark/Distance_3/bench_compare_distance_3.cpp b/Distance_3/benchmark/Distance_3/bench_compare_distance_3.cpp
new file mode 100644
index 00000000000..d7005a594ab
--- /dev/null
+++ b/Distance_3/benchmark/Distance_3/bench_compare_distance_3.cpp
@@ -0,0 +1,297 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Simple_homogeneous.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_integer.h>
+#include <CGAL/Homogeneous.h>
+
+#include <CGAL/squared_distance_3.h>
+
+#include <CGAL/Cartesian_converter.h>
+
+#include <CGAL/Random.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
+
+// #define CGAL_USE_GTE_AS_SANITY_CHECK
+#ifdef CGAL_USE_GTE_AS_SANITY_CHECK
+#include <Mathematics/DistTriangle3Triangle3.h>
+#include <Mathematics/DistSegmentSegment.h>
+#endif
+
+#include <cassert>
+#include <iostream>
+
+struct randomint
+{
+  randomint() ;
+  int get() const { return sequence[cur]; }
+  int next() {
+    cur = (cur + 1) % 11;
+    return get();
+  }
+
+private:
+  int sequence[11];
+  int cur;
+};
+
+inline randomint::randomint()
+{
+  cur = 0;
+  sequence[0] = 19;
+  sequence[1] = 5;
+  sequence[2] = 17;
+  sequence[3] = 13;
+  sequence[4] = 29;
+  sequence[5] = 2;
+  sequence[6] = 23;
+  sequence[7] = 31;
+  sequence[8] = 3;
+  sequence[9] = 37;
+  sequence[10] = 11;
+}
+
+randomint ri;
+
+template <typename K>
+struct Test
+{
+  typedef typename K::RT              RT;
+  typedef typename K::FT              FT;
+  typedef typename K::Comparison_result Comparison_result;
+  typedef typename K::Point_3         P;
+  typedef typename K::Segment_3       S;
+  typedef typename K::Vector_3        V;
+  typedef typename K::Ray_3           R;
+  typedef typename K::Line_3          L;
+  typedef typename K::Triangle_3      T;
+  typedef typename K::Plane_3         Pl;
+  typedef typename K::Tetrahedron_3   Tet;
+  typedef typename K::Iso_cuboid_3    Cub;
+
+private:
+  CGAL::Random& r;
+  double m = 0, M = 1;
+
+public:
+   Test(CGAL::Random& r) : r(r) { }
+
+private:
+  inline RT to_nt(int d) const { return RT(d); }
+
+  P p(int x, int y, int z)
+  {
+    int w = ri.next();
+    return P(to_nt(x*w), to_nt(y*w), to_nt(z*w), to_nt(w));
+  }
+
+  P random_point() const
+  {
+    return P(FT(r.get_double(m, M)), FT(r.get_double(m, M)), FT(r.get_double(m, M)));
+  }
+
+  Pl pl(int a, int b, int c, int d)
+  {
+    int w = ri.next();
+    return Pl(to_nt(a*w), to_nt(b*w), to_nt(c*w), to_nt(d*w));
+  }
+
+
+private:
+  void P_P(int N, FT d2)
+  {
+    std::cout << "Point - Point" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      K().compare_squared_distance_3_object()(p0, p1, d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void P_S(int N, FT d2)
+  {
+    std::cout << "Point - Segment" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      K().compare_squared_distance_3_object()(p0, S(p1,p2), d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void P_T(int N, FT d2)
+  {
+    std::cout << "Point - Triangle" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P q = random_point();
+
+      K().compare_squared_distance_3_object()(q, T(p0, p1, p2), d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void P_Tet(int N, FT d2)
+  {
+    std::cout << "Point - Tetrahedron\n";
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P p3 = random_point();
+      P q = random_point();
+
+      K().compare_squared_distance_3_object()(q, Tet(p0, p1, p2, p3), d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void S_S(int N, FT d2)
+  {
+    std::cout << "Segment - Segment" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P q0 = random_point();
+      P q1 = random_point();
+
+      K().compare_squared_distance_3_object()(S(p0, p1), S(q0, q1), d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+
+  }
+
+  void S_L(int N, FT d2)
+  {
+    std::cout << "Segment - Line" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P q0 = random_point();
+      P q1 = random_point();
+
+      K().compare_squared_distance_3_object()(S(p0, p1), L(q0, q1), d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+
+  }
+
+  void T_T(int N, FT d2)
+  {
+    std::cout << "Triangle - Triangle with distance " << d2 << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P p3 = random_point();
+      P p4 = random_point();
+      P p5 = random_point();
+
+      // // these are still exact with EPECK
+      // p0 = CGAL::midpoint(p0, p1);
+      // p1 = p0 + FT(FT(0.1)) * V{p1 - p0};
+      // p2 = p2 + V{p2 - p0} / FT(CGAL_PI);
+
+      // // this is still exact with EPECK_with_sqrt
+      // p4 = p4 + V{p4 - CGAL::ORIGIN} / CGAL::approximate_sqrt(CGAL::square(p4.x()) + CGAL::square(p4.y()) + CGAL::square(p4.z()) + 3);
+
+      // p5 = p5 + V{p5 - CGAL::ORIGIN} * FT(std::cos(1.3));
+      // generic triangles
+      T tr1{p0, p1, p2}, tr2{p3, p4, p5};
+      K().compare_squared_distance_3_object()(tr1, tr2, d2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+public:
+  void run()
+  {
+    std::cout << "Kernel: " << typeid(K).name() << std::endl;
+    auto test_3=[]<typename F>(F& f){
+
+    };
+    P_P(10000000, FT(0.1));
+    P_S(1000000, FT(0.1));
+
+    std::cout << std::endl;
+    P_T(500000, FT(10));
+    P_T(500000, FT(0.1));
+    P_T(500000, FT(0.001));
+    std::cout << std::endl;
+
+    P_Tet(200000, FT(0.1));
+
+    std::cout << std::endl;
+    S_S(500000, FT(10));
+    S_S(500000, FT(0.1));
+    S_S(500000, FT(0.001));
+    std::cout << std::endl;
+    S_L(500000, FT(0.1));
+
+    std::cout << std::endl;
+    T_T(500000, FT(10));
+    T_T(500000, FT(0.1));
+    T_T(500000, FT(0.001));
+    std::cout << std::endl;
+    std::cout << std::endl;
+  }
+};
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  std::cout << "3D Distance tests" << std::endl;
+
+  CGAL::Random rp;
+  CGAL::Random r(argc==1?rp.get_seed():std::stoi(argv[1]));
+  std::cout << "random seed = " << r.get_seed() << std::endl;
+
+ Test<CGAL::Simple_cartesian<double> >(r).run();
+ Test<CGAL::Simple_homogeneous<double> >(r).run();
+//  Test<CGAL::Simple_cartesian<CGAL::Interval_nt<true> > >(r).run();
+
+  // Test<CGAL::Homogeneous<CGAL::Exact_integer> >(r).run();
+
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>(r).run();
+
+  Test<CGAL::Exact_predicates_exact_constructions_kernel>(r).run();
+
+  std::cout << "Done!" << std::endl;
+
+  return EXIT_SUCCESS;
+}
diff --git a/Distance_3/benchmark/Distance_3/bench_distance_3.cpp b/Distance_3/benchmark/Distance_3/bench_distance_3.cpp
new file mode 100644
index 00000000000..641bd01f0ae
--- /dev/null
+++ b/Distance_3/benchmark/Distance_3/bench_distance_3.cpp
@@ -0,0 +1,280 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Simple_homogeneous.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_integer.h>
+#include <CGAL/Homogeneous.h>
+
+#include <CGAL/squared_distance_3.h>
+
+#include <CGAL/Cartesian_converter.h>
+
+#include <CGAL/Random.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Real_timer.h>
+
+// #define CGAL_USE_GTE_AS_SANITY_CHECK
+#ifdef CGAL_USE_GTE_AS_SANITY_CHECK
+#include <Mathematics/DistTriangle3Triangle3.h>
+#include <Mathematics/DistSegmentSegment.h>
+#endif
+
+#include <cassert>
+#include <iostream>
+
+struct randomint
+{
+  randomint() ;
+  int get() const { return sequence[cur]; }
+  int next() {
+    cur = (cur + 1) % 11;
+    return get();
+  }
+
+private:
+  int sequence[11];
+  int cur;
+};
+
+inline randomint::randomint()
+{
+  cur = 0;
+  sequence[0] = 19;
+  sequence[1] = 5;
+  sequence[2] = 17;
+  sequence[3] = 13;
+  sequence[4] = 29;
+  sequence[5] = 2;
+  sequence[6] = 23;
+  sequence[7] = 31;
+  sequence[8] = 3;
+  sequence[9] = 37;
+  sequence[10] = 11;
+}
+
+randomint ri;
+
+template <typename K>
+struct Test
+{
+  typedef typename K::RT              RT;
+  typedef typename K::FT              FT;
+  typedef typename K::Comparison_result Comparison_result;
+  typedef typename K::Point_3         P;
+  typedef typename K::Segment_3       S;
+  typedef typename K::Vector_3        V;
+  typedef typename K::Ray_3           R;
+  typedef typename K::Line_3          L;
+  typedef typename K::Triangle_3      T;
+  typedef typename K::Plane_3         Pl;
+  typedef typename K::Tetrahedron_3   Tet;
+  typedef typename K::Iso_cuboid_3    Cub;
+
+private:
+  CGAL::Random& r;
+  double m = 0, M = 1;
+
+public:
+   Test(CGAL::Random& r) : r(r) { }
+
+private:
+  inline RT to_nt(int d) const { return RT(d); }
+
+  P p(int x, int y, int z)
+  {
+    int w = ri.next();
+    return P(to_nt(x*w), to_nt(y*w), to_nt(z*w), to_nt(w));
+  }
+
+  P random_point() const
+  {
+    return P(FT(r.get_double(m, M)), FT(r.get_double(m, M)), FT(r.get_double(m, M)));
+  }
+
+  Pl pl(int a, int b, int c, int d)
+  {
+    int w = ri.next();
+    return Pl(to_nt(a*w), to_nt(b*w), to_nt(c*w), to_nt(d*w));
+  }
+
+
+private:
+  void P_P(int N)
+  {
+    std::cout << "Point - Point" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      CGAL::squared_distance(p0, p1);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void P_S(int N)
+  {
+    std::cout << "Point - Segment" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      CGAL::squared_distance(p0, S(p1,p2));
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void P_T(int N)
+  {
+    std::cout << "Point - Triangle" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P q = random_point();
+
+      CGAL::squared_distance(q, T(p0, p1, p2));
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void P_Tet(int N)
+  {
+    std::cout << "Point - Tetrahedron\n";
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P p3 = random_point();
+      P q = random_point();
+
+      CGAL::squared_distance(q, Tet(p0, p1, p2, p3));
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+  void S_S(int N)
+  {
+    std::cout << "Segment - Segment" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P q0 = random_point();
+      P q1 = random_point();
+
+      CGAL::squared_distance(S(p0, p1), S(q0, q1));
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+
+  }
+
+  void S_L(int N)
+  {
+    std::cout << "Segment - Line" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P q0 = random_point();
+      P q1 = random_point();
+
+      CGAL::squared_distance(S(p0, p1), L(q0, q1));
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+
+  }
+
+  void T_T(int N)
+  {
+    std::cout << "Triangle - Triangle" << std::endl;
+    CGAL::Real_timer t;
+    t.start();
+    for(int i=0; i<N; ++i)
+    {
+      P p0 = random_point();
+      P p1 = random_point();
+      P p2 = random_point();
+      P p3 = random_point();
+      P p4 = random_point();
+      P p5 = random_point();
+
+      // // these are still exact with EPECK
+      // p0 = CGAL::midpoint(p0, p1);
+      // p1 = p0 + FT(0.1) * V{p1 - p0};
+      // p2 = p2 + V{p2 - p0} / FT(CGAL_PI);
+
+      // // this is still exact with EPECK_with_sqrt
+      // p4 = p4 + V{p4 - CGAL::ORIGIN} / CGAL::approximate_sqrt(CGAL::square(p4.x()) + CGAL::square(p4.y()) + CGAL::square(p4.z()) + 3);
+
+      // p5 = p5 + V{p5 - CGAL::ORIGIN} * FT(std::cos(1.3));
+      // generic triangles
+      T tr1{p0, p1, p2}, tr2{p3, p4, p5};
+      CGAL::squared_distance(tr1, tr2);
+    }
+    t.stop();
+    std::cout << t.time() << " sec" << std::endl;
+  }
+
+public:
+  void run()
+  {
+    std::cout << "Kernel: " << typeid(K).name() << std::endl;
+
+    P_P(10000000);
+    P_S(1000000);
+    P_T(500000);
+    P_Tet(200000);
+
+    S_S(500000);
+    S_L(500000);
+    T_T(100000);
+    std::cout << std::endl;
+  }
+};
+
+int main(int argc, char** argv)
+{
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  std::cout << "3D Distance tests" << std::endl;
+
+  CGAL::Random rp;
+  CGAL::Random r(argc==1?rp.get_seed():std::stoi(argv[1]));
+  std::cout << "random seed = " << r.get_seed() << std::endl;
+
+ Test<CGAL::Simple_cartesian<double> >(r).run();
+ Test<CGAL::Simple_homogeneous<double> >(r).run();
+//  Test<CGAL::Simple_cartesian<CGAL::Interval_nt<true> > >(r).run();
+
+  // Test<CGAL::Homogeneous<CGAL::Exact_integer> >(r).run();
+
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>(r).run();
+
+  Test<CGAL::Exact_predicates_exact_constructions_kernel>(r).run();
+
+  std::cout << "Done!" << std::endl;
+
+  return EXIT_SUCCESS;
+}
diff --git a/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp b/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp
new file mode 100644
index 00000000000..70fe503d513
--- /dev/null
+++ b/Distance_3/benchmark/Distance_3/bench_thingi10k_3.cpp
@@ -0,0 +1,140 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Simple_homogeneous.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_integer.h>
+#include <CGAL/Homogeneous.h>
+
+#include <CGAL/squared_distance_3.h>
+#include <CGAL/Real_timer.h>
+
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/IO/polygon_soup_io.h>
+
+#include <boost/container/small_vector.hpp>
+
+#include <cassert>
+#include <iostream>
+
+
+template <typename K>
+struct Test
+{
+  typedef typename K::RT              RT;
+  typedef typename K::FT              FT;
+  typedef typename K::Comparison_result Comparison_result;
+  typedef typename K::Point_3         P;
+  typedef typename K::Segment_3       S;
+  typedef typename K::Vector_3        V;
+  typedef typename K::Ray_3           R;
+  typedef typename K::Line_3          L;
+  typedef typename K::Triangle_3      T;
+  typedef typename K::Plane_3         Pl;
+  typedef typename K::Tetrahedron_3   Tet;
+  typedef typename K::Iso_cuboid_3    Cub;
+
+  typedef std::vector<boost::container::small_vector<std::size_t, 3> >::iterator Iterator;
+  typedef CGAL::Box_intersection_d::Box_with_handle_d<double,3,Iterator> CBox;
+
+  size_t nb_closed_pairs;
+
+public:
+   Test() : nb_closed_pairs(0){ }
+
+  void close_triangles(std::vector<P> &points, std::vector<boost::container::small_vector<std::size_t, 3> >& triangles, FT d2){
+    std::vector< CBox > boxes;
+    auto extend_bbox3=[&](const Iterator it, FT& d2){
+      CGAL::Bbox_3 bb = points[(*it)[0]].bbox()+points[(*it)[1]].bbox()+points[(*it)[2]].bbox();
+      return CGAL::Bbox_3(bb.xmin(),bb.ymin(),bb.zmin(),bb.xmax()+CGAL::to_double(d2),bb.ymax()+CGAL::to_double(d2),bb.zmax()+CGAL::to_double(d2));
+    };
+
+    auto callback=[&](const CBox &ba, const CBox &bb){
+      boost::container::small_vector<std::size_t, 3> &a = *(ba.handle());
+      boost::container::small_vector<std::size_t, 3> &b = *(bb.handle());
+
+      std::sort(a.begin(), a.end());
+      std::sort(b.begin(), b.end());
+      std::vector<int> v;
+      std::set_intersection(a.begin(), a.end(), b.begin(), b.end(),std::back_inserter(v));
+
+      if(v.size()!=0) //they have common vertices
+        return; 
+
+      bool comp = K().compare_squared_distance_3_object()(T(points[a[0]], points[a[1]], points[a[2]]),
+                                                 T(points[b[0]], points[b[1]], points[b[2]]),
+                                                 d2)!=CGAL::LARGER;
+      if(comp)
+      {
+        nb_closed_pairs++;
+      }
+    };
+    
+    for(Iterator it=triangles.begin(); it!=triangles.end(); ++it)
+      boxes.emplace_back(extend_bbox3(it, d2), it);
+
+    CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), callback);
+  }
+
+  void run(std::string filename, FT d2)
+  {
+    nb_closed_pairs=0;
+    std::vector<P> input_points;
+    std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+
+    if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+    {
+      std::cerr << "Cannot read " << filename << "\n";
+      return;
+    }
+    CGAL::Real_timer t;
+    t.start();
+    close_triangles(input_points, input_triangles, d2);
+    t.stop();
+    std::cout << "#points = " << input_points.size() << " and #triangles = " << input_triangles.size() << " has " << nb_closed_pairs << " pairs at squared distance " << d2 << " in " << t.time() << " sec." << std::endl;
+  }
+};
+
+int main(int argc, char** argv)
+{
+  const std::string filename = argc == 1 ? CGAL::data_file_path("meshes/elephant.off")
+                                         : std::string(argv[1]);
+
+  // const std::string out_file = argc <= 2 ? "rounded_soup.off"
+  //                                        : std::string(argv[2]);
+
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
+  std::cout << "3D Distance bench" << std::endl;
+
+  std::vector<CGAL::Exact_predicates_inexact_constructions_kernel::Point_3> input_points;
+  std::vector<boost::container::small_vector<std::size_t, 3>> input_triangles;
+
+  if (!CGAL::IO::read_polygon_soup(filename, input_points, input_triangles))
+  {
+    std::cerr << "Cannot read " << filename << "\n";
+    return 1;
+  }
+  CGAL::Bbox_3 bb = CGAL::bbox_3(input_points.begin(), input_points.end());
+  double max= (std::max)((std::max)(bb.xmax()-bb.xmin(),bb.ymax()-bb.ymin()),bb.zmax()-bb.zmin());
+
+//  Test<CGAL::Simple_cartesian<double> >().run(filename);
+//  Test<CGAL::Simple_homogeneous<double> >().run(filename);
+//  Test<CGAL::Simple_cartesian<CGAL::Interval_nt<true> > >(r).run();
+
+  // Test<CGAL::Homogeneous<CGAL::Exact_integer> >(r).run();
+
+  // Test<CGAL::Exact_predicates_inexact_constructions_kernel>().run(filename, max*max);
+  // Test<CGAL::Exact_predicates_inexact_constructions_kernel>().run(filename, 10*max*max/input_points.size());
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>().run(filename, 100*max/input_points.size());
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>().run(filename, max/input_points.size());
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>().run(filename, max/(100*input_points.size()));
+
+  Test<CGAL::Exact_predicates_exact_constructions_kernel>().run(filename, 100*max/input_points.size());
+  Test<CGAL::Exact_predicates_exact_constructions_kernel>().run(filename, max/input_points.size());
+  Test<CGAL::Exact_predicates_exact_constructions_kernel>().run(filename, max/(100*input_points.size()));
+
+  std::cout << "Done!" << std::endl;
+
+  return EXIT_SUCCESS;
+}