diff --git a/Distance_2/include/CGAL/Distance_2/Line_2_Line_2.h b/Distance_2/include/CGAL/Distance_2/Line_2_Line_2.h
index b87124a8ee9..53fef5cbd53 100644
--- a/Distance_2/include/CGAL/Distance_2/Line_2_Line_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Line_2_Line_2.h
@@ -20,7 +20,6 @@
 #define CGAL_DISTANCE_2_LINE_2_LINE_2_H
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
-#include <CGAL/Distance_2/Point_2_Line_2.h>
 
 #include <CGAL/Kernel/global_functions.h>
 #include <CGAL/Line_2.h>
@@ -35,8 +34,11 @@ squared_distance(const typename K::Line_2& line1,
                  const K& k)
 {
   typedef typename K::FT FT;
+
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
+
   if(internal::parallel(line1, line2, k))
-    return internal::squared_distance(line1.point(), line2, k);
+    return sq_dist(line1.point(), line2);
   else
     return FT(0);
 }
@@ -48,7 +50,7 @@ inline typename K::FT
 squared_distance(const Line_2<K>& line1,
                  const Line_2<K>& line2)
 {
-  return internal::squared_distance(line1, line2, K());
+  return K().compute_squared_distance_2_object()(line1, line2);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Line_2_Triangle_2.h b/Distance_2/include/CGAL/Distance_2/Line_2_Triangle_2.h
index 6ba20ce2eef..3f9fc69ac3d 100644
--- a/Distance_2/include/CGAL/Distance_2/Line_2_Triangle_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Line_2_Triangle_2.h
@@ -20,7 +20,6 @@
 #define CGAL_DISTANCE_2_LINE_2_TRIANGLE_2_H
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
-#include <CGAL/Distance_2/Point_2_Line_2.h>
 
 #include <CGAL/enum.h>
 
@@ -38,17 +37,19 @@ squared_distance(const typename K::Line_2& line,
 {
   typedef typename K::FT FT;
 
-  Oriented_side side0 = line.oriented_side(triangle.vertex(0));
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
+
+  const Oriented_side side0 = line.oriented_side(triangle.vertex(0));
   if(line.oriented_side(triangle.vertex(1)) != side0)
     return FT(0);
 
   if(line.oriented_side(triangle.vertex(2)) != side0)
     return FT(0);
 
-  FT mindist = internal::squared_distance(triangle.vertex(0), line, k);
+  FT mindist = sq_dist(triangle.vertex(0), line);
   for(int i=1; i<3; ++i)
   {
-    FT dist = internal::squared_distance(triangle.vertex(i), line, k);
+    FT dist = sq_dist(triangle.vertex(i), line);
     if(dist < mindist)
       mindist = dist;
   }
@@ -72,7 +73,7 @@ inline typename K::FT
 squared_distance(const Line_2<K>& line,
                  const Triangle_2<K>& triangle)
 {
-  return internal::squared_distance(line, triangle, K());
+  return K().compute_squared_distance_2_object()(line, triangle);
 }
 
 template <class K>
@@ -80,7 +81,7 @@ inline typename K::FT
 squared_distance(const Triangle_2<K>& triangle,
                  const Line_2<K>& line)
 {
-  return internal::squared_distance(line, triangle, K());
+  return K().compute_squared_distance_2_object()(triangle, line);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Point_2_Line_2.h b/Distance_2/include/CGAL/Distance_2/Point_2_Line_2.h
index 43f1ff622ff..61550cb6908 100644
--- a/Distance_2/include/CGAL/Distance_2/Point_2_Line_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Point_2_Line_2.h
@@ -92,7 +92,7 @@ inline typename K::FT
 squared_distance(const Point_2<K>& pt,
                  const Line_2<K>& line)
 {
-  return internal::squared_distance(pt, line, K());
+  return K().compute_squared_distance_2_object()(pt, line);
 }
 
 template <class K>
@@ -100,7 +100,7 @@ inline typename K::FT
 squared_distance(const Line_2<K>& line,
                  const Point_2<K>& pt)
 {
-  return internal::squared_distance(pt, line, K());
+  return K().compute_squared_distance_2_object()(line, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Point_2_Point_2.h b/Distance_2/include/CGAL/Distance_2/Point_2_Point_2.h
index d8ae488e0a7..b7405fddd6b 100644
--- a/Distance_2/include/CGAL/Distance_2/Point_2_Point_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Point_2_Point_2.h
@@ -41,7 +41,7 @@ inline typename K::FT
 squared_distance(const Point_2<K>& pt1,
                  const Point_2<K>& pt2)
 {
-  return internal::squared_distance(pt1, pt2, K());
+  return K().compute_squared_distance_2_object()(pt1, pt2);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Point_2_Segment_2.h b/Distance_2/include/CGAL/Distance_2/Point_2_Segment_2.h
index c8dd3486141..cc591908a2b 100644
--- a/Distance_2/include/CGAL/Distance_2/Point_2_Segment_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Point_2_Segment_2.h
@@ -20,8 +20,6 @@
 #define CGAL_DISTANCE_2_POINT_2_SEGMENT_2_H
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
-#include <CGAL/Distance_2/Point_2_Point_2.h>
-#include <CGAL/Distance_2/Point_2_Line_2.h>
 
 #include <CGAL/Point_2.h>
 #include <CGAL/Segment_2.h>
@@ -57,11 +55,15 @@ squared_distance_indexed(const typename K::Point_2 &pt,
                          int ind,
                          const K& k)
 {
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
+
   if(ind == 0)
-    return internal::squared_distance(pt, seg.source(), k);
+    return sq_dist(pt, seg.source());
+
   if(ind == 1)
-    return internal::squared_distance(pt, seg.target(), k);
-  return internal::squared_distance(pt, seg.supporting_line(), k);
+    return sq_dist(pt, seg.target());
+
+  return sq_dist(pt, seg.supporting_line());
 }
 
 template <class K>
@@ -73,21 +75,23 @@ squared_distance(const typename K::Point_2& pt,
   typedef typename K::Vector_2 Vector_2;
   typedef typename K::RT RT;
 
-  typename K::Construct_vector_2 construct_vector = k.construct_vector_2_object();
+  typename K::Construct_vector_2 vector = k.construct_vector_2_object();
+  typename K::Compute_squared_length_2 sq_length = k.compute_squared_length_2_object();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
 
   // assert that the segment is valid (non zero length).
-  Vector_2 diff = construct_vector(seg.source(), pt);
-  Vector_2 segvec = construct_vector(seg.source(), seg.target());
+  const Vector_2 diff = vector(seg.source(), pt);
+  const Vector_2 segvec = vector(seg.source(), seg.target());
 
-  RT d = wdot(diff, segvec, k);
+  const RT d = wdot(diff, segvec, k);
   if(d <= RT(0))
-    return k.compute_squared_length_2_object()(diff);
+    return sq_length(diff);
 
-  RT e = wdot(segvec,segvec, k);
+  const RT e = wdot(segvec,segvec, k);
   if(wmult((K*)0 ,d, segvec.hw()) > wmult((K*)0, e, diff.hw()))
-    return internal::squared_distance(pt, seg.target(), k);
+    return sq_dist(pt, seg.target());
 
-  return internal::squared_distance(pt, seg.supporting_line(), k);
+  return sq_dist(pt, seg.supporting_line());
 }
 
 template <class K>
@@ -106,7 +110,7 @@ inline typename K::FT
 squared_distance(const Point_2<K>& pt,
                  const Segment_2<K>& seg)
 {
-  return internal::squared_distance(pt, seg, K());
+  return K().compute_squared_distance_2_object()(pt, seg);
 }
 
 template <class K>
@@ -114,7 +118,7 @@ inline typename K::FT
 squared_distance(const Segment_2<K>& seg,
                  const Point_2<K>& pt)
 {
-  return internal::squared_distance(pt, seg, K());
+  return K().compute_squared_distance_2_object()(seg, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Point_2_Triangle_2.h b/Distance_2/include/CGAL/Distance_2/Point_2_Triangle_2.h
index 9cd6c83fe14..2448737e116 100644
--- a/Distance_2/include/CGAL/Distance_2/Point_2_Triangle_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Point_2_Triangle_2.h
@@ -212,13 +212,15 @@ squared_distance_indexed(const typename K::Point_2& pt,
   typedef typename K::FT      FT;
   typedef typename K::Line_2  Line_2;
 
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
+
   if(ind1 == -1)
     return FT(0);
 
   if(ind2 == -1)
-    return internal::squared_distance(pt, triangle.vertex(ind1), k);
+    return sq_dist(pt, triangle.vertex(ind1));
 
-  return internal::squared_distance(pt, Line_2(triangle.vertex(ind1), triangle.vertex(ind2)), k);
+  return sq_dist(pt, Line_2{triangle.vertex(ind1), triangle.vertex(ind2)});
 }
 
 template <class K>
@@ -248,7 +250,7 @@ inline typename K::FT
 squared_distance(const Point_2<K>& pt,
                  const Triangle_2<K>& triangle)
 {
-  return internal::squared_distance(pt, triangle, K());
+  return K().compute_squared_distance_2_object()(pt, triangle);
 }
 
 template <class K>
@@ -256,7 +258,7 @@ inline typename K::FT
 squared_distance(const Triangle_2<K>& triangle,
                  const Point_2<K>& pt)
 {
-  return internal::squared_distance(pt, triangle, K());
+  return K().compute_squared_distance_2_object()(triangle, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Ray_2_Line_2.h b/Distance_2/include/CGAL/Distance_2/Ray_2_Line_2.h
index bef607cf383..ea3c8fb7af0 100644
--- a/Distance_2/include/CGAL/Distance_2/Ray_2_Line_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Ray_2_Line_2.h
@@ -72,7 +72,7 @@ inline typename K::FT
 squared_distance(const Line_2<K>& line,
                  const Ray_2<K>& ray)
 {
-  return internal::squared_distance(line, ray, K());
+  return K().compute_squared_distance_2_object()(line, ray);
 }
 
 template <class K>
@@ -80,7 +80,7 @@ inline typename K::FT
 squared_distance(const Ray_2<K>& ray,
                  const Line_2<K>& line)
 {
-  return internal::squared_distance(line, ray, K());
+  return K().compute_squared_distance_2_object()(ray, line);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Ray_2_Ray_2.h b/Distance_2/include/CGAL/Distance_2/Ray_2_Ray_2.h
index 87fdf1903f6..f8d046f602a 100644
--- a/Distance_2/include/CGAL/Distance_2/Ray_2_Ray_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Ray_2_Ray_2.h
@@ -46,7 +46,7 @@ ray_ray_squared_distance_parallel(const typename K::Vector_2& ray1dir,
   RT wcr = wcross(ray1dir, from1to2, k);
   RT w = from1to2.hw();
 
-  return (square(wcr) / FT(wmult((K*)0, RT(wdot(ray1dir, ray1dir, k)), w, w)));
+  return (square(wcr) / FT(wmult((K*)0, wdot(ray1dir, ray1dir, k), w, w)));
 }
 
 template <class K>
@@ -58,11 +58,12 @@ squared_distance(const typename K::Ray_2& ray1,
   typedef typename K::Vector_2 Vector_2;
   typedef typename K::FT FT;
 
-  typename K::Construct_vector_2 construct_vector = k.construct_vector_2_object();
+  typename K::Construct_vector_2 vector = k.construct_vector_2_object();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
 
-  const Vector_2& ray1dir = ray1.direction().vector();
-  const Vector_2& ray2dir = ray2.direction().vector();
-  Vector_2 diffvec(construct_vector(ray1.source(),ray2.source()));
+  const Vector_2 ray1dir = ray1.direction().vector();
+  const Vector_2 ray2dir = ray2.direction().vector();
+  const Vector_2 diffvec = vector(ray1.source(),ray2.source());
 
   bool crossing1, crossing2;
   switch(orientation(ray1dir, ray2dir, k))
@@ -83,18 +84,18 @@ squared_distance(const typename K::Ray_2& ray1,
   {
     if(crossing2)
       return FT(0);
-    return internal::squared_distance(ray2.source(), ray1, k);
+    return sq_dist(ray2.source(), ray1);
   }
   else
   {
     if(crossing2)
     {
-      return internal::squared_distance(ray1.source(), ray2, k);
+      return sq_dist(ray1.source(), ray2);
     }
     else
     {
-      FT min1 = internal::squared_distance(ray1.source(), ray2, k);
-      FT min2 = internal::squared_distance(ray2.source(), ray1, k);
+      FT min1 = sq_dist(ray1.source(), ray2);
+      FT min2 = sq_dist(ray2.source(), ray1);
       return (min1 < min2) ? min1 : min2;
     }
   }
@@ -106,7 +107,7 @@ template <class K>
 inline typename K::FT
 squared_distance(const Ray_2<K> &ray1, const Ray_2<K> &ray2)
 {
-  return internal::squared_distance(ray1, ray2, K());
+  return K().compute_squared_distance_2_object()(ray1, ray2);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Ray_2_Triangle_2.h b/Distance_2/include/CGAL/Distance_2/Ray_2_Triangle_2.h
index f549ce3d880..8c7bbd0be3f 100644
--- a/Distance_2/include/CGAL/Distance_2/Ray_2_Triangle_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Ray_2_Triangle_2.h
@@ -20,6 +20,7 @@
 #define CGAL_DISTANCE_2_RAY_2_TRIANGLE_2_H
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
+#include <CGAL/Distance_2/Point_2_Triangle_2.h>
 
 #include <CGAL/Ray_2.h>
 #include <CGAL/Triangle_2.h>
@@ -64,11 +65,10 @@ squared_distance(const typename K::Ray_2& ray,
   {
     // Check if all the segment vertices lie at the same side of
     // the triangle segment.
-    const Point_2 &vt1 = triangle.vertex(ind_tr1);
-    const Point_2 &vt2 = triangle.vertex(ind_tr2);
+    const Point_2& vt1 = triangle.vertex(ind_tr1);
+    const Point_2& vt2 = triangle.vertex(ind_tr2);
     if(clockwise(ray.direction().vector(), vt2-vt1, k))
       mindist = FT(0);
-
   }
   else
   {
@@ -105,7 +105,7 @@ inline typename K::FT
 squared_distance(const Ray_2<K>& ray,
                  const Triangle_2<K>& triangle)
 {
-  return internal::squared_distance(ray, triangle, K());
+  return K().compute_squared_distance_2_object()(ray, triangle);
 }
 
 template <class K>
@@ -113,7 +113,7 @@ inline typename K::FT
 squared_distance(const Triangle_2<K>& triangle,
                  const Ray_2<K>& ray)
 {
-  return internal::squared_distance(ray, triangle, K());
+  return K().compute_squared_distance_2_object()(triangle, ray);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Segment_2_Line_2.h b/Distance_2/include/CGAL/Distance_2/Segment_2_Line_2.h
index b5ba324207f..a4d8537b08a 100644
--- a/Distance_2/include/CGAL/Distance_2/Segment_2_Line_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Segment_2_Line_2.h
@@ -21,6 +21,7 @@
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
 #include <CGAL/Distance_2/Point_2_Line_2.h>
+#include <CGAL/Distance_2/Segment_2_Ray_2.h>
 
 #include <CGAL/number_utils.h>
 #include <CGAL/Rational_traits.h>
@@ -56,20 +57,20 @@ squared_distance(const typename K::Segment_2& seg,
   typedef typename K::Vector_2 Vector_2;
   typedef typename K::Point_2  Point_2;
 
-  typename K::Construct_vector_2 construct_vector = k.construct_vector_2_object();
+  typename K::Construct_vector_2 vector = k.construct_vector_2_object();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
 
-  const Vector_2 &linedir = line.direction().vector();
-  const Point_2 &linepoint = line.point();
-  Vector_2 startvec(construct_vector(linepoint, seg.source()));
-  Vector_2 endvec(construct_vector(linepoint, seg.target()));
+  const Vector_2 linedir = line.direction().vector();
+  const Point_2& linepoint = line.point();
+  const Vector_2 startvec = vector(linepoint, seg.source());
+  const Vector_2 endvec = vector(linepoint, seg.target());
+
+  if(seg.source() == seg.target())
+    return sq_dist(seg.source(), line);
 
   bool crossing1;
-  RT c1s, c1e;
-  if(seg.source() == seg.target())
-    return internal::squared_distance(seg.source(), line, k);
-
-  c1s = wcross(linedir, startvec, k);
-  c1e = wcross(linedir, endvec, k);
+  const RT c1s = wcross(linedir, startvec, k);
+  const RT c1e = wcross(linedir, endvec, k);
   if(c1s < RT(0))
   {
     crossing1 = (c1e >= RT(0));
@@ -88,8 +89,7 @@ squared_distance(const typename K::Segment_2& seg,
   }
   else
   {
-    RT dm;
-    dm = _distance_measure_sub<K>(c1s, c1e, startvec, endvec);
+    const RT dm = _distance_measure_sub<K>(c1s, c1e, startvec, endvec);
     if(dm <= RT(0))
       return _sqd_to_line<K>(startvec, c1s, linedir);
      else
@@ -113,7 +113,7 @@ inline typename K::FT
 squared_distance(const Segment_2<K>& seg,
                  const Line_2<K>& line)
 {
-  return internal::squared_distance(seg, line, K());
+  return K().compute_squared_distance_2_object()(seg, line);
 }
 
 template <class K>
@@ -121,7 +121,7 @@ inline typename K::FT
 squared_distance(const Line_2<K>& line,
                  const Segment_2<K>& seg)
 {
-  return internal::squared_distance(seg, line, K());
+  return K().compute_squared_distance_2_object()(line, seg);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Segment_2_Ray_2.h b/Distance_2/include/CGAL/Distance_2/Segment_2_Ray_2.h
index 1c5c300648d..0c17e595c01 100644
--- a/Distance_2/include/CGAL/Distance_2/Segment_2_Ray_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Segment_2_Ray_2.h
@@ -50,48 +50,49 @@ squared_distance_parallel(const typename K::Segment_2& seg,
 {
   typedef typename K::Vector_2 Vector_2;
 
-  const Vector_2& dir1 = seg.direction().vector();
-  const Vector_2& dir2 = ray.direction().vector();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
+
+  const Vector_2 dir1 = seg.direction().vector();
+  const Vector_2 dir2 = ray.direction().vector();
 
   if(same_direction(dir1, dir2, k))
   {
     if(!is_acute_angle(seg.source(), seg.target(), ray.source(), k))
-      return internal::squared_distance(seg.target(), ray.source(), k);
+      return sq_dist(seg.target(), ray.source());
   }
   else
   {
     if(!is_acute_angle(seg.target(), seg.source(), ray.source(), k))
-      return internal::squared_distance(seg.source(), ray.source(), k);
+      return sq_dist(seg.source(), ray.source());
   }
 
-  return internal::squared_distance(ray.source(), seg.supporting_line(), k);
+  return sq_dist(ray.source(), seg.supporting_line());
 }
 
 template <class K>
 typename K::FT
-squared_distance(const typename K::Segment_2 &seg,
-                 const typename K::Ray_2 &ray,
+squared_distance(const typename K::Segment_2& seg,
+                 const typename K::Ray_2& ray,
                  const K& k)
 {
   typedef typename K::RT RT;
   typedef typename K::FT FT;
   typedef typename K::Vector_2 Vector_2;
 
-  typename K::Construct_vector_2 construct_vector = k.construct_vector_2_object();
+  typename K::Construct_vector_2 vector = k.construct_vector_2_object();
   typename K::Orientation_2 orientation = k.orientation_2_object();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
 
   const Vector_2& raydir = ray.direction().vector();
-  Vector_2 startvec(construct_vector(ray.source(), seg.source()));
-  Vector_2 endvec(construct_vector(ray.source(), seg.target()));
-
-  bool crossing1, crossing2;
-  RT c1s, c1e;
+  const Vector_2 startvec = vector(ray.source(), seg.source());
+  const Vector_2 endvec = vector(ray.source(), seg.target());
 
   if(seg.source() == seg.target())
     return internal::squared_distance(seg.source(), ray, k);
 
-  c1s = wcross(raydir, startvec, k);
-  c1e = wcross(raydir, endvec, k);
+  bool crossing1, crossing2;
+  const RT c1s = wcross(raydir, startvec, k);
+  const RT c1e = wcross(raydir, endvec, k);
   if(c1s < RT(0))
   {
     crossing1 = (c1e >= RT(0));
@@ -114,10 +115,10 @@ squared_distance(const typename K::Segment_2 &seg,
   switch (orientation(seg.source(), seg.target(), ray.source()))
   {
     case LEFT_TURN:
-      crossing2 = right_turn(construct_vector(seg.source(), seg.target()), raydir, k);
+      crossing2 = right_turn(vector(seg.source(), seg.target()), raydir, k);
       break;
     case RIGHT_TURN:
-      crossing2 = left_turn(construct_vector(seg.source(), seg.target()), raydir, k);
+      crossing2 = left_turn(vector(seg.source(), seg.target()), raydir, k);
       break;
     default:
       crossing2 = true;
@@ -129,36 +130,34 @@ squared_distance(const typename K::Segment_2 &seg,
     if(crossing2)
       return FT(0);
 
-    return internal::squared_distance(ray.source(), seg, k);
+    return sq_dist(ray.source(), seg);
   }
   else
   {
     if(crossing2)
     {
-      RT dm;
-      dm = _distance_measure_sub<K>(c1s, c1e, startvec, endvec);
+      const RT dm = _distance_measure_sub<K>(c1s, c1e, startvec, endvec);
       if(dm < RT(0))
       {
-        return internal::squared_distance(seg.source(), ray, k);
+        return sq_dist(seg.source(), ray);
       }
       else
       {
         if(dm > RT(0))
-          return internal::squared_distance(seg.target(), ray, k);
+          return sq_dist(seg.target(), ray);
         else // parallel, should not happen (no crossing)
           return internal::squared_distance_parallel(seg, ray, k);
       }
     }
     else
     {
-      FT min1, min2;
-      RT dm = _distance_measure_sub<K>(c1s, c1e, startvec, endvec);
+      const RT dm = _distance_measure_sub<K>(c1s, c1e, startvec, endvec);
       if(dm == RT(0))
         return internal::squared_distance_parallel(seg, ray, k);
 
-      min1 = (dm < RT(0)) ? internal::squared_distance(seg.source(), ray, k)
-                          : internal::squared_distance(seg.target(), ray, k);
-      min2 = internal::squared_distance(ray.source(), seg, k);
+      const FT min1 = (dm < RT(0)) ? sq_dist(seg.source(), ray)
+                                   : sq_dist(seg.target(), ray);
+      const FT min2 = sq_dist(ray.source(), seg);
 
       return (min1 < min2) ? min1 : min2;
     }
@@ -181,7 +180,7 @@ inline typename K::FT
 squared_distance(const Segment_2<K>& seg,
                  const Ray_2<K>& ray)
 {
-  return internal::squared_distance(seg, ray, K());
+  return K().compute_squared_distance_2_object()(seg, ray);
 }
 
 template <class K>
@@ -189,7 +188,7 @@ inline typename K::FT
 squared_distance(const Ray_2<K>& ray,
                  const Segment_2<K>& seg)
 {
-  return internal::squared_distance(seg, ray, K());
+  return K().compute_squared_distance_2_object()(ray, seg);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Segment_2_Segment_2.h b/Distance_2/include/CGAL/Distance_2/Segment_2_Segment_2.h
index 5757e39e637..23e2f07c827 100644
--- a/Distance_2/include/CGAL/Distance_2/Segment_2_Segment_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Segment_2_Segment_2.h
@@ -40,27 +40,29 @@ squared_distance_parallel(const typename K::Segment_2& seg1,
 {
   typedef typename K::Vector_2 Vector_2;
 
-  const Vector_2 &dir1 = seg1.direction().vector();
-  const Vector_2 &dir2 = seg2.direction().vector();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
+
+  const Vector_2 dir1 = seg1.direction().vector();
+  const Vector_2 dir2 = seg2.direction().vector();
 
   if(same_direction(dir1, dir2, k))
   {
     if(!is_acute_angle(seg1.source(), seg1.target(), seg2.source(), k))
-      return internal::squared_distance(seg1.target(), seg2.source(), k);
+      return sq_dist(seg1.target(), seg2.source());
 
     if(!is_acute_angle(seg1.target(), seg1.source(), seg2.target(), k))
-      return internal::squared_distance(seg1.source(), seg2.target(), k);
+      return sq_dist(seg1.source(), seg2.target());
   }
   else
   {
     if(!is_acute_angle(seg1.source(), seg1.target(), seg2.target(), k))
-      return internal::squared_distance(seg1.target(), seg2.target(), k);
+      return sq_dist(seg1.target(), seg2.target());
 
     if(!is_acute_angle(seg1.target(), seg1.source(), seg2.source(), k))
-      return internal::squared_distance(seg1.source(), seg2.source(), k);
+      return sq_dist(seg1.source(), seg2.source());
   }
 
-  return internal::squared_distance(seg2.source(), seg1.supporting_line(), k);
+  return sq_dist(seg2.source(), seg1.supporting_line());
 }
 
 template <class K>
@@ -84,17 +86,19 @@ squared_distance(const typename K::Segment_2& seg1,
   typedef typename K::RT RT;
   typedef typename K::FT FT;
 
-  bool crossing1, crossing2;
-  RT c1s, c1e, c2s, c2e;
+  typename K::Orientation_2 orientation = k.orientation_2_object();
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
 
   if(seg1.source() == seg1.target())
-    return internal::squared_distance(seg1.source(), seg2, k);
+    return sq_dist(seg1.source(), seg2);
 
   if(seg2.source() == seg2.target())
-    return internal::squared_distance(seg2.source(), seg1, k);
+    return sq_dist(seg2.source(), seg1);
 
-  Orientation o1s = orientation(seg2.source(), seg2.target(), seg1.source());
-  Orientation o1e = orientation(seg2.source(), seg2.target(), seg1.target());
+  const Orientation o1s = orientation(seg2.source(), seg2.target(), seg1.source());
+  const Orientation o1e = orientation(seg2.source(), seg2.target(), seg1.target());
+
+  bool crossing1, crossing2;
   if(o1s == RIGHT_TURN)
   {
     crossing1 = (o1e != RIGHT_TURN);
@@ -114,8 +118,8 @@ squared_distance(const typename K::Segment_2& seg1,
     }
   }
 
-  Orientation o2s = orientation(seg1.source(), seg1.target(), seg2.source());
-  Orientation o2e = orientation(seg1.source(), seg1.target(), seg2.target());
+  const Orientation o2s = orientation(seg1.source(), seg1.target(), seg2.source());
+  const Orientation o2e = orientation(seg1.source(), seg1.target(), seg2.target());
   if(o2s == RIGHT_TURN)
   {
     crossing2 = (o2e != RIGHT_TURN);
@@ -140,19 +144,19 @@ squared_distance(const typename K::Segment_2& seg1,
     if(crossing2)
       return (FT)0;
 
-    c2s = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.source(), k));
-    c2e = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.target(), k));
-    Comparison_result dm = CGAL_NTS compare(c2s,c2e);
+    const RT c2s = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.source(), k));
+    const RT c2e = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.target(), k));
+    Comparison_result dm = CGAL_NTS compare(c2s, c2e);
 
     if(dm == SMALLER)
     {
-      return internal::squared_distance(seg2.source(), seg1, k);
+      return sq_dist(seg2.source(), seg1);
     }
     else
     {
       if(dm == LARGER)
       {
-        return internal::squared_distance(seg2.target(), seg1, k);
+        return sq_dist(seg2.target(), seg1);
       }
       else
       {
@@ -163,42 +167,40 @@ squared_distance(const typename K::Segment_2& seg1,
   }
   else
   {
-    c1s = CGAL_NTS abs(wcross(seg2.source(), seg2.target(), seg1.source(), k));
-    c1e = CGAL_NTS abs(wcross(seg2.source(), seg2.target(), seg1.target(), k));
-    Comparison_result dm = CGAL_NTS compare(c1s,c1e);
+    const RT c1s = CGAL_NTS abs(wcross(seg2.source(), seg2.target(), seg1.source(), k));
+    const RT c1e = CGAL_NTS abs(wcross(seg2.source(), seg2.target(), seg1.target(), k));
+    Comparison_result dm = CGAL_NTS compare(c1s, c1e);
     if(crossing2)
     {
       if(dm == SMALLER)
       {
-        return internal::squared_distance(seg1.source(), seg2, k);
+        return sq_dist(seg1.source(), seg2);
       }
       else
       {
         if(dm == LARGER)
-          return internal::squared_distance(seg1.target(), seg2, k);
+          return sq_dist(seg1.target(), seg2);
         else  // parallel, should not happen (no crossing)
           return internal::squared_distance_parallel(seg1, seg2, k);
       }
     }
     else
     {
-      FT min1, min2;
-
       if(dm == EQUAL)
         return internal::squared_distance_parallel(seg1, seg2, k);
 
-      min1 = (dm == SMALLER) ? internal::squared_distance(seg1.source(), seg2, k):
-                               internal::squared_distance(seg1.target(), seg2, k);
+      FT min1 = (dm == SMALLER) ? sq_dist(seg1.source(), seg2):
+                                  sq_dist(seg1.target(), seg2);
+
+      const RT c2s = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.source(), k));
+      const RT c2e = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.target(), k));
 
-      c2s = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.source(), k));
-      c2e = CGAL_NTS abs(wcross(seg1.source(), seg1.target(), seg2.target(), k));
       dm = CGAL_NTS compare(c2s,c2e);
-
       if(dm == EQUAL)  // should not happen.
         return internal::squared_distance_parallel(seg1, seg2, k);
 
-      min2 = (dm == SMALLER) ? internal::squared_distance(seg2.source(), seg1, k):
-                               internal::squared_distance(seg2.target(), seg1, k);
+      FT min2 = (dm == SMALLER) ? sq_dist(seg2.source(), seg1):
+                                  sq_dist(seg2.target(), seg1);
 
       return (min1 < min2) ? min1 : min2;
     }
@@ -207,27 +209,28 @@ squared_distance(const typename K::Segment_2& seg1,
 
 template <class K>
 typename K::FT
-squared_distance(const typename K::Segment_2 &seg1,
-                 const typename K::Segment_2 &seg2,
+squared_distance(const typename K::Segment_2& seg1,
+                 const typename K::Segment_2& seg2,
                  const K& k,
                  const Homogeneous_tag&)
 {
   typedef typename K::RT RT;
   typedef typename K::FT FT;
 
-  bool crossing1, crossing2;
-  RT c1s, c1e, c2s, c2e;
+  typename K::Compute_squared_distance_2 sq_dist = k.compute_squared_distance_2_object();
 
   if(seg1.source() == seg1.target())
-    return internal::squared_distance(seg1.source(), seg2, k);
+    return sq_dist(seg1.source(), seg2);
 
   if(seg2.source() == seg2.target())
-    return internal::squared_distance(seg2.source(), seg1, k);
+    return sq_dist(seg2.source(), seg1);
 
-  c1s = wcross(seg2.source(), seg2.target(), seg1.source(), k);
-  c1e = wcross(seg2.source(), seg2.target(), seg1.target(), k);
-  c2s = wcross(seg1.source(), seg1.target(), seg2.source(), k);
-  c2e = wcross(seg1.source(), seg1.target(), seg2.target(), k);
+  const RT c1s = wcross(seg2.source(), seg2.target(), seg1.source(), k);
+  const RT c1e = wcross(seg2.source(), seg2.target(), seg1.target(), k);
+  const RT c2s = wcross(seg1.source(), seg1.target(), seg2.source(), k);
+  const RT c2e = wcross(seg1.source(), seg1.target(), seg2.target(), k);
+
+  bool crossing1, crossing2;
 
   if(c1s < RT(0))
   {
@@ -272,16 +275,15 @@ squared_distance(const typename K::Segment_2 &seg1,
     if(crossing2)
       return (FT)0;
 
-    RT dm;
-    dm = _distance_measure_sub<K>(c2s,c2e, seg2.source(), seg2.target());
+    const RT dm = _distance_measure_sub<K>(c2s,c2e, seg2.source(), seg2.target());
     if(dm < RT(0))
     {
-      return internal::squared_distance(seg2.source(), seg1, k);
+      return sq_dist(seg2.source(), seg1);
     }
     else
     {
       if(dm > RT(0))
-        return internal::squared_distance(seg2.target(), seg1, k);
+        return sq_dist(seg2.target(), seg1);
       else // parallel, should not happen (no crossing)
         return internal::squared_distance_parallel(seg1, seg2, k);
     }
@@ -290,41 +292,50 @@ squared_distance(const typename K::Segment_2 &seg1,
   {
     if(crossing2)
     {
-      RT dm = _distance_measure_sub<K>(c1s, c1e,seg1.source(),seg1.target());
+      const RT dm = _distance_measure_sub<K>(c1s, c1e,seg1.source(),seg1.target());
       if(dm < RT(0))
       {
-        return internal::squared_distance(seg1.source(), seg2, k);
+        return sq_dist(seg1.source(), seg2);
       }
       else
       {
         if(dm > RT(0))
-          return internal::squared_distance(seg1.target(), seg2, k);
+          return sq_dist(seg1.target(), seg2);
         else // parallel, should not happen (no crossing)
           return internal::squared_distance_parallel(seg1, seg2, k);
       }
     }
     else
     {
-      FT min1, min2;
       RT dm = _distance_measure_sub<K>(c1s, c1e, seg1.source(), seg1.target());
       if(dm == RT(0))
         return internal::squared_distance_parallel(seg1, seg2, k);
 
-      min1 = (dm < RT(0)) ? internal::squared_distance(seg1.source(), seg2, k):
-                            internal::squared_distance(seg1.target(), seg2, k);
-      dm = _distance_measure_sub<K>(c2s, c2e, seg2.source(), seg2.target());
+      FT min1 = (dm < RT(0)) ? sq_dist(seg1.source(), seg2):
+                               sq_dist(seg1.target(), seg2);
 
+      dm = _distance_measure_sub<K>(c2s, c2e, seg2.source(), seg2.target());
       if(dm == RT(0))  // should not happen.
         return internal::squared_distance_parallel(seg1, seg2, k);
 
-      min2 = (dm < RT(0)) ? internal::squared_distance(seg2.source(), seg1, k):
-                            internal::squared_distance(seg2.target(), seg1, k);
+      FT min2 = (dm < RT(0)) ? sq_dist(seg2.source(), seg1):
+                               sq_dist(seg2.target(), seg1);
 
       return (min1 < min2) ? min1 : min2;
     }
   }
 }
 
+template <class K>
+inline typename K::FT
+squared_distance(const Segment_2<K>& seg1,
+                 const Segment_2<K>& seg2,
+                 const K& k)
+{
+  typedef typename K::Kernel_tag Tag;
+  return squared_distance(seg1, seg2, k, Tag());
+}
+
 } // namespace internal
 
 template <class K>
@@ -332,8 +343,7 @@ inline typename K::FT
 squared_distance(const Segment_2<K>& seg1,
                  const Segment_2<K>& seg2)
 {
-  typedef typename K::Kernel_tag Tag;
-  return internal::squared_distance(seg1, seg2, K(), Tag());
+  return K().compute_squared_distance_2_object()(seg1, seg2);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Segment_2_Triangle_2.h b/Distance_2/include/CGAL/Distance_2/Segment_2_Triangle_2.h
index c8234370526..13424998bf8 100644
--- a/Distance_2/include/CGAL/Distance_2/Segment_2_Triangle_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Segment_2_Triangle_2.h
@@ -20,6 +20,7 @@
 #define CGAL_DISTANCE_2_SEGMENT_2_TRIANGLE_2_H
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
+#include <CGAL/Distance_2/Point_2_Triangle_2.h>
 
 #include <CGAL/Segment_2.h>
 #include <CGAL/Triangle_2.h>
@@ -79,7 +80,7 @@ squared_distance(const typename K::Segment_2& seg,
     // the triangle segment.
     const Point_2 &vt1 = triangle.vertex(ind_tr1);
     const Point_2 &vt2 = triangle.vertex(ind_tr2);
-    Orientation or_s = orientation(vt1, vt2, seg.source());
+    const Orientation or_s = orientation(vt1, vt2, seg.source());
     if(orientation(vt1, vt2, seg.target()) != or_s)
       mindist = FT(0);
   }
@@ -87,9 +88,9 @@ squared_distance(const typename K::Segment_2& seg,
   {
     // Check if all the triangle vertices lie
     // at the same side of the segment.
-    const Point_2 &vt1 = seg.source();
-    const Point_2 &vt2 = seg.target();
-    Orientation or_s = orientation(vt1, vt2, triangle.vertex(0));
+    const Point_2& vt1 = seg.source();
+    const Point_2& vt2 = seg.target();
+    const Orientation or_s = orientation(vt1, vt2, triangle.vertex(0));
     for(i=1; i<3; ++i)
     {
       if(orientation(vt1, vt2, triangle.vertex(i)) != or_s)
@@ -119,7 +120,7 @@ inline typename K::FT
 squared_distance(const Segment_2<K>& seg,
                  const Triangle_2<K>& triangle)
 {
-  return internal::squared_distance(seg, triangle, K());
+  return K().compute_squared_distance_2_object()(seg, triangle);
 }
 
 template <class K>
@@ -127,7 +128,7 @@ inline typename K::FT
 squared_distance(const Triangle_2<K>& triangle,
                  const Segment_2<K>& seg)
 {
-  return internal::squared_distance(seg, triangle, K());
+  return K().compute_squared_distance_2_object()(triangle, seg);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/include/CGAL/Distance_2/Triangle_2_Triangle_2.h b/Distance_2/include/CGAL/Distance_2/Triangle_2_Triangle_2.h
index 67cb787b215..b308ccff35f 100644
--- a/Distance_2/include/CGAL/Distance_2/Triangle_2_Triangle_2.h
+++ b/Distance_2/include/CGAL/Distance_2/Triangle_2_Triangle_2.h
@@ -21,6 +21,7 @@
 
 #include <CGAL/Distance_2/internal/squared_distance_utils_2.h>
 #include <CGAL/Distance_2/Point_2_Point_2.h>
+#include <CGAL/Distance_2/Point_2_Triangle_2.h>
 
 #include <CGAL/Triangle_2.h>
 
@@ -68,8 +69,7 @@ squared_distance(const typename K::Triangle_2& triangle1,
     }
   }
 
-  // now check if all vertices are on the right side of the
-  // separating line.
+  // now check if all vertices are on the right side of the separating line.
   if(ind1_2 == -1 && ind2_2 == -1)
     return mindist;
 
@@ -80,7 +80,7 @@ squared_distance(const typename K::Triangle_2& triangle1,
   {
     const Point_2& vt1 = triangle1.vertex(ind1_1);
     const Point_2& vt2 = triangle1.vertex(ind1_2);
-    Orientation or_s = orientation(vt1, vt2, triangle2.vertex(0));
+    const Orientation or_s = orientation(vt1, vt2, triangle2.vertex(0));
     for(int i=1; i<3; ++i)
     {
       if(orientation(vt1, vt2, triangle2.vertex(i)) != or_s)
@@ -94,7 +94,7 @@ squared_distance(const typename K::Triangle_2& triangle1,
   {
     const Point_2& vt1 = triangle2.vertex(ind2_1);
     const Point_2& vt2 = triangle2.vertex(ind2_2);
-    Orientation or_s = orientation(vt1, vt2, triangle1.vertex(0));
+    const Orientation or_s = orientation(vt1, vt2, triangle1.vertex(0));
     for(int i=1; i<3; ++i)
     {
       if(orientation(vt1, vt2, triangle1.vertex(i)) != or_s)
@@ -115,7 +115,7 @@ inline typename K::FT
 squared_distance(const Triangle_2<K>& triangle1,
                  const Triangle_2<K>& triangle2)
 {
-  return internal::squared_distance(triangle1, triangle2, K());
+  return K().compute_squared_distance_2_object()(triangle1, triangle2);
 }
 
 } // namespace CGAL
diff --git a/Distance_2/test/Distance_2/test_distance_2.cpp b/Distance_2/test/Distance_2/test_distance_2.cpp
index 371086c6185..2736f8fe551 100644
--- a/Distance_2/test/Distance_2/test_distance_2.cpp
+++ b/Distance_2/test/Distance_2/test_distance_2.cpp
@@ -1,12 +1,12 @@
-// 2D distance tests.
-
 #ifdef NDEBUG
 #undef NDEBUG //this testsuite requires NDEBUG to be not defined
 #endif
 
-
 #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/Homogeneous.h>
 
 #include <vector>
 #include <iostream>
@@ -268,7 +268,8 @@ struct Test {
 
   void run()
   {
-    std::cout << "2D Distance tests\n";
+    std::cout << "-- Kernel: " << typeid(K).name() << std::endl;
+
     P_P();
     P_S();
     P_R();
@@ -293,8 +294,12 @@ struct Test {
 
 int main()
 {
-  Test< CGAL::Simple_cartesian<double>   >().run();
-  Test< CGAL::Simple_homogeneous<double> >().run();
+  std::cout << "2D Distance tests\n";
 
-  // TODO : test more kernels.
+  Test<CGAL::Simple_cartesian<double> >().run();
+  Test<CGAL::Simple_homogeneous<double> >().run();
+
+  Test<CGAL::Homogeneous<CGAL::Epeck_ft> >().run();
+  Test<CGAL::Exact_predicates_inexact_constructions_kernel>().run();
+  Test<CGAL::Exact_predicates_exact_constructions_kernel>().run();
 }
diff --git a/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h
index 88a69f727e3..03495be04fa 100644
--- a/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Line_3_Line_3.h
@@ -34,11 +34,11 @@ squared_distance(const typename K::Line_3& line1,
 
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
 
-  Vector_3 dir1, dir2, normal, diff;
-  dir1 = line1.direction().vector();
-  dir2 = line2.direction().vector();
-  normal = wcross(dir1, dir2, k);
-  diff = vector(line1.point(), line2.point());
+  const Vector_3 dir1 = line1.direction().vector();
+  const Vector_3 dir2 = line2.direction().vector();
+  const Vector_3 normal = wcross(dir1, dir2, k);
+  const Vector_3 diff = vector(line1.point(), line2.point());
+
   if (is_null(normal, k))
     return squared_distance_to_line(dir2, diff, k);
 
@@ -53,7 +53,7 @@ typename K::FT
 squared_distance(const Line_3<K>& line1,
                  const Line_3<K>& line2)
 {
-  return internal::squared_distance(line1, line2, K());
+  return K().compute_squared_distance_3_object()(line1, line2);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h
index 69bea4dcbe5..6bb4d125c59 100644
--- a/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Line_3_Plane_3.h
@@ -18,7 +18,6 @@
 #define CGAL_DISTANCE_3_LINE_3_PLANE_3_H
 
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
-#include <CGAL/Distance_3/Point_3_Point_3.h>
 
 #include <CGAL/Line_3.h>
 #include <CGAL/Plane_3.h>
@@ -45,8 +44,10 @@ squared_distance(const typename K::Line_3& l,
 {
   typedef typename K::FT FT;
 
-  if (contains_vector(pl, l.direction().vector(), k))
-    return squared_distance(pl, l.point(), k);
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
+
+  if(contains_vector(pl, l.direction().vector(), k))
+    return sq_dist(pl, l.point());
 
   return FT(0);
 }
@@ -68,7 +69,7 @@ typename K::FT
 squared_distance(const Line_3<K>& line,
                  const Plane_3<K>& plane)
 {
-  return internal::squared_distance(line, plane, K());
+  return K().compute_squared_distance_3_object()(line, plane);
 }
 
 template <class K>
@@ -77,7 +78,7 @@ typename K::FT
 squared_distance(const Plane_3<K>& plane,
                  const Line_3<K>& line)
 {
-  return internal::squared_distance(line, plane, K());
+  return K().compute_squared_distance_3_object()(plane, line);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h
index 8feada6266d..d18fb366b7f 100644
--- a/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Plane_3_Plane_3.h
@@ -18,11 +18,28 @@
 #define CGAL_DISTANCE_3_PLANE_3_PLANE_3_H
 
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
-#include <CGAL/Distance_3/Point_3_Plane_3.h>
 
 #include <CGAL/Plane_3.h>
 
 namespace CGAL {
+namespace internal {
+
+template <class K>
+inline typename K::FT
+squared_distance(const typename K::Plane_3& plane1,
+                 const typename K::Plane_3& plane2,
+                 const K& k)
+{
+  typename K::Construct_orthogonal_vector_3 ortho_vec = k.construct_orthogonal_vector_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
+
+  if(!is_null(wcross(ortho_vec(plane1), ortho_vec(plane2), k), k))
+    return typename K::FT(0);
+  else
+    return sq_dist(plane1.point(), plane2);
+}
+
+} // namespace internal
 
 template <class K>
 inline
@@ -30,13 +47,7 @@ typename K::FT
 squared_distance(const Plane_3<K>& plane1,
                  const Plane_3<K>& plane2)
 {
-  K k;
-  typename K::Construct_orthogonal_vector_3 ortho_vec = k.construct_orthogonal_vector_3_object();
-
-  if (!internal::is_null(internal::wcross(ortho_vec(plane1), ortho_vec(plane2), k), k))
-    return typename K::FT(0);
-  else
-    return internal::squared_distance(plane1.point(), plane2, k);
+  return K().compute_squared_distance_3_object()(plane1, plane2);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h
index 6ac84fc2223..24f7593d51d 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Line_3.h
@@ -33,9 +33,10 @@ squared_distance(const typename K::Point_3& pt,
 {
   typedef typename K::Vector_3 Vector_3;
 
-  typename K::Construct_vector_3 construct_vector;
-  Vector_3 dir(line.direction().vector());
-  Vector_3 diff = construct_vector(line.point(), pt);
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+
+  const Vector_3 dir(line.direction().vector());
+  const Vector_3 diff = vector(line.point(), pt);
 
   return internal::squared_distance_to_line(dir, diff, k);
 }
@@ -58,7 +59,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Line_3<K>& line)
 {
-  return internal::squared_distance(pt, line, K());
+  return K().compute_squared_distance_3_object()(pt, line);
 }
 
 template <class K>
@@ -67,7 +68,7 @@ typename K::FT
 squared_distance(const Line_3<K>& line,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(pt, line, K());
+  return K().compute_squared_distance_3_object()(line, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h
index 6d0da93dfe0..7544a72d963 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Plane_3.h
@@ -56,7 +56,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Plane_3<K>& plane)
 {
-  return internal::squared_distance(pt, plane, K());
+  return K().compute_squared_distance_3_object()(pt, plane);
 }
 
 template <class K>
@@ -65,7 +65,7 @@ typename K::FT
 squared_distance(const Plane_3<K>& plane,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(pt, plane, K());
+  return K().compute_squared_distance_3_object()(plane, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h
index a3c4ab12dc9..91c8965d26f 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Ray_3.h
@@ -33,13 +33,13 @@ squared_distance(const typename K::Point_3& pt,
 {
   typedef typename K::Vector_3 Vector_3;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
 
-  Vector_3 diff = construct_vector(ray.source(), pt);
-  const Vector_3 &dir = ray.direction().vector();
+  const Vector_3 diff = vector(ray.source(), pt);
+  const Vector_3 dir = ray.direction().vector();
 
   if(!is_acute_angle(dir, diff, k))
-    return (typename K::FT)(diff*diff);
+    return diff*diff;
 
   return squared_distance_to_line(dir, diff, k);
 }
@@ -61,7 +61,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Ray_3<K>& ray)
 {
-  return internal::squared_distance(pt, ray, K());
+  return K().compute_squared_distance_3_object()(pt, ray);
 }
 
 template <class K>
@@ -70,7 +70,7 @@ typename K::FT
 squared_distance(const Ray_3<K>& ray,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(pt, ray, K());
+  return K().compute_squared_distance_3_object()(ray, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h
index a28dd10a33c..34ee169c8e1 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Segment_3.h
@@ -37,21 +37,22 @@ squared_distance(const typename K::Point_3& pt,
   typedef typename K::RT RT;
   typedef typename K::FT FT;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   // assert that the segment is valid (non zero length).
-  Vector_3 diff = construct_vector(seg.source(), pt);
-  Vector_3 segvec = construct_vector(seg.source(), seg.target());
+  Vector_3 diff = vector(seg.source(), pt);
+  Vector_3 segvec = vector(seg.source(), seg.target());
 
   RT d = wdot(diff,segvec, k);
   if(d <= RT(0))
-    return FT(diff*diff);
+    return diff*diff;
   RT e = wdot(segvec,segvec, k);
   if((d * segvec.hw()) > (e * diff.hw()))
-    return squared_distance(pt, seg.target(), k);
+    return sq_dist(pt, seg.target());
 
   Vector_3 wcr = wcross(segvec, diff, k);
-  return FT(wcr*wcr)/FT(e * diff.hw() * diff.hw());
+  return FT(wcr*wcr) / FT(e * diff.hw() * diff.hw());
 }
 
 template <class K>
@@ -65,19 +66,20 @@ squared_distance(const typename K::Point_3& pt,
   typedef typename K::RT RT;
   typedef typename K::FT FT;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+  typename K::Compute_squared_distance_3 sqd = k.compute_squared_distance_3_object();
 
   // assert that the segment is valid (non zero length).
-  Vector_3 diff = construct_vector(seg.source(), pt);
-  Vector_3 segvec = construct_vector(seg.source(), seg.target());
+  Vector_3 diff = vector(seg.source(), pt);
+  Vector_3 segvec = vector(seg.source(), seg.target());
 
   RT d = wdot(diff,segvec, k);
   if(d <= RT(0))
-    return (FT(diff*diff));
+    return diff*diff;
 
   RT e = wdot(segvec,segvec, k);
   if(d > e)
-    return squared_distance(pt, seg.target(), k);
+    return sqd(pt, seg.target());
 
   Vector_3 wcr = wcross(segvec, diff, k);
   return FT(wcr*wcr)/e;
@@ -113,7 +115,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Segment_3<K>& seg)
 {
-  return internal::squared_distance(pt, seg, K());
+  return K().compute_squared_distance_3_object()(pt, seg);
 }
 
 template <class K>
@@ -122,7 +124,7 @@ typename K::FT
 squared_distance(const Segment_3<K>& seg,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(pt, seg, K());
+  return K().compute_squared_distance_3_object()(seg, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h
index c2bedd0aba1..04d148818a3 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Tetrahedron_3.h
@@ -132,7 +132,7 @@ typename K::FT
 squared_distance(const Tetrahedron_3<K>& tet,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(pt, tet, K());
+  return K().compute_squared_distance_3_object()(tet, pt);
 }
 
 template <class K>
@@ -140,7 +140,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Tetrahedron_3<K>& tet)
 {
-  return internal::squared_distance(pt, tet, K());
+  return K().compute_squared_distance_3_object()(pt, tet);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h
index 4231977b22d..8859831c6f0 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Triangle_3.h
@@ -26,6 +26,7 @@
 namespace CGAL {
 namespace internal {
 
+// returns true iff pt is on the negative side of the plane defined by (ep0, ep1) and normal
 template <class K>
 inline bool
 on_left_of_triangle_edge(const typename K::Point_3& pt,
@@ -34,16 +35,15 @@ on_left_of_triangle_edge(const typename K::Point_3& pt,
                          const typename K::Point_3& ep1,
                          const K& k)
 {
-  // returns true iff pt is on the negative side of the plane defined by (ep0, ep1) and normal
+  typedef typename K::RT RT;
+  typedef typename K::Vector_3 Vector_3;
 
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
-  typename K::Vector_3 edge = vector(ep0, ep1);
-  typename K::Vector_3 diff = vector(ep0, pt);
 
-  typedef typename K::RT RT;
+  const Vector_3 edge = vector(ep0, ep1);
+  const Vector_3 diff = vector(ep0, pt);
 
-  const bool result = RT(wdot(wcross(edge, normal, k), diff, k)) <= RT(0);
-  return result;
+  return (wdot(wcross(edge, normal, k), diff, k) <= RT(0));
 }
 
 template <class K>
@@ -59,6 +59,7 @@ squared_distance_to_triangle(const typename K::Point_3& pt,
 
   typename K::Construct_segment_3 segment = k.construct_segment_3_object();
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   const Vector_3 e1 = vector(t0, t1);
   const Vector_3 oe3 = vector(t0, t2);
@@ -74,24 +75,24 @@ squared_distance_to_triangle(const typename K::Point_3& pt,
     // and only two distances could be used, but leaving 3 for the case of
     // inexact constructions as it might improve the accuracy.
 
-    typename K::FT d1 = internal::squared_distance(pt, segment(t2, t0), k);
-    typename K::FT d2 = internal::squared_distance(pt, segment(t1, t2), k);
-    typename K::FT d3 = internal::squared_distance(pt, segment(t0, t1), k);
+    typename K::FT d1 = sq_dist(pt, segment(t2, t0));
+    typename K::FT d2 = sq_dist(pt, segment(t1, t2));
+    typename K::FT d3 = sq_dist(pt, segment(t0, t1));
 
     return (std::min)( (std::min)(d1, d2), d3);
   }
 
   const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k);
   if(!b01)
-    return internal::squared_distance(pt, segment(t0, t1), k);
+    return sq_dist(pt, segment(t0, t1));
 
   const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k);
   if(!b12)
-    return internal::squared_distance(pt, segment(t1, t2), k);
+    return sq_dist(pt, segment(t1, t2));
 
   const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k);
   if(!b20)
-    return internal::squared_distance(pt, segment(t2, t0), k);
+    return sq_dist(pt, segment(t2, t0));
 
   // The projection of pt is inside the triangle
   inside = true;
@@ -115,6 +116,15 @@ squared_distance(const typename K::Point_3& pt,
                                       unused_inside);
 }
 
+template <class K>
+inline typename K::FT
+squared_distance(const typename K::Triangle_3& t,
+                 const typename K::Point_3& pt,
+                 const K& k)
+{
+  return squared_distance(pt, t, k);
+}
+
 } // namespace internal
 
 template <class K>
@@ -123,7 +133,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Triangle_3<K>& t)
 {
-  return internal::squared_distance(pt, t, K());
+  return K().compute_squared_distance_3_object()(pt, t);
 }
 
 template <class K>
@@ -132,7 +142,7 @@ typename K::FT
 squared_distance(const Triangle_3<K>& t,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(pt, t, K());
+  return K().compute_squared_distance_3_object()(t, pt);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Point_3_Weighted_point_3.h b/Distance_3/include/CGAL/Distance_3/Point_3_Weighted_point_3.h
index 5f19203a64b..f5db1c3c46f 100644
--- a/Distance_3/include/CGAL/Distance_3/Point_3_Weighted_point_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Point_3_Weighted_point_3.h
@@ -30,7 +30,7 @@ typename K::FT
 squared_distance(const Weighted_point_3<K>& wpt,
                  const Point_3<K>& pt)
 {
-  return internal::squared_distance(wpt.point(), pt, K());
+  return K().compute_squared_distance_3_object()(wpt.point(), pt);
 }
 
 template <class K>
@@ -39,7 +39,7 @@ typename K::FT
 squared_distance(const Point_3<K>& pt,
                  const Weighted_point_3<K>& wpt)
 {
-  return internal::squared_distance(pt, wpt.point(), K());
+  return K().compute_squared_distance_3_object()(pt, wpt.point());
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Ray_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Ray_3_Line_3.h
index 63b98004ac0..d4b51b59cbe 100644
--- a/Distance_3/include/CGAL/Distance_3/Ray_3_Line_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Ray_3_Line_3.h
@@ -27,60 +27,65 @@ namespace internal {
 
 template <class K>
 typename K::FT
-squared_distance(const typename K::Line_3& line,
-                 const typename K::Ray_3& ray,
+squared_distance(const typename K::Ray_3& ray,
+                 const typename K::Line_3& line,
                  const K& k)
 {
   typedef typename K::Vector_3 Vector_3;
   typedef typename K::Point_3 Point_3;
   typedef typename K::RT RT;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
 
   const Point_3& rs = ray.source();
 
-  Vector_3 raydir, linedir, normal;
-  linedir = line.direction().vector();
-  raydir = ray.direction().vector();
-  normal = wcross(raydir, linedir, k);
+  const Vector_3 linedir = line.direction().vector();
+  const Vector_3 raydir = ray.direction().vector();
+  const Vector_3 normal = wcross(raydir, linedir, k);
 
-  Vector_3 rs_min_lp = construct_vector(line.point(), rs);
-  if (is_null(normal, k))
+  const Vector_3 rs_min_lp = vector(line.point(), rs);
+  if(is_null(normal, k))
     return squared_distance_to_line(linedir, rs_min_lp, k);
 
   bool crossing;
-  RT sdm_sr_l;
-  Vector_3 perpend2l;
-  perpend2l = wcross(linedir, normal, k);
-
-  sdm_sr_l = wdot(perpend2l, rs_min_lp, k);
-  if (sdm_sr_l < RT(0))
+  const Vector_3 perpend2l = wcross(linedir, normal, k);
+  const RT sdm_sr_l = wdot(perpend2l, rs_min_lp, k);
+  if(sdm_sr_l < RT(0))
   {
-    crossing = (RT(wdot(perpend2l, raydir, k)) >= RT(0));
+    crossing = (wdot(perpend2l, raydir, k) >= RT(0));
   }
   else
   {
-    if (RT(wdot(perpend2l, raydir, k)) <= RT(0))
+    if(wdot(perpend2l, raydir, k) <= RT(0))
       crossing = true;
     else
       crossing = (sdm_sr_l == RT(0));
   }
 
-  if (crossing)
+  if(crossing)
     return squared_distance_to_plane(normal, rs_min_lp, k);
   else
     return squared_distance_to_line(linedir, rs_min_lp, k);
 }
 
+template <class K>
+typename K::FT
+squared_distance(const typename K::Line_3& line,
+                 const typename K::Ray_3& ray,
+                 const K& k)
+{
+  return squared_distance(ray, line, k);
+}
+
 } // namespace internal
 
 template <class K>
 inline
 typename K::FT
-squared_distance(const Line_3<K> &line,
-                 const Ray_3<K> &ray)
+squared_distance(const Line_3<K>& line,
+                 const Ray_3<K>& ray)
 {
-  return internal::squared_distance(line, ray, K());
+  return K().compute_squared_distance_3_object()(line, ray);
 }
 
 template <class K>
@@ -89,7 +94,7 @@ typename K::FT
 squared_distance(const Ray_3<K>& ray,
                  const Line_3<K>& line)
 {
-  return internal::squared_distance(line, ray, K());
+  return K().compute_squared_distance_3_object()(ray, line);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Ray_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Ray_3_Plane_3.h
index c72ee2975ee..bdc9a956146 100644
--- a/Distance_3/include/CGAL/Distance_3/Ray_3_Plane_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Ray_3_Plane_3.h
@@ -18,10 +18,10 @@
 #define CGAL_DISTANCE_3_RAY_3_PLANE_3_H
 
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
+#include <CGAL/number_utils.h>
 
 #include <CGAL/Plane_3.h>
 #include <CGAL/Ray_3.h>
-#include <CGAL/number_utils.h>
 
 namespace CGAL {
 namespace internal {
@@ -32,32 +32,32 @@ squared_distance(const typename K::Ray_3 &ray,
                  const typename K::Plane_3 &plane,
                  const K& k)
 {
-  typedef typename K::Point_3 Point_3;
-  typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object();
 
-  const Point_3 &start = ray.start();
-  const Point_3 &planepoint = plane.point();
-  Vector_3 start_min_pp = construct_vector(planepoint, start);
-  Vector_3 end_min_pp = ray.direction().vector();
-  const Vector_3 &normal = plane.orthogonal_vector();
-  RT sdm_rs2pp = wdot(normal, start_min_pp, k);
-  RT sdm_re2pp = wdot(normal, end_min_pp, k);
+  const Point_3& start = ray.start();
+  const Point_3& planepoint = plane.point();
+  const Vector_3 start_min_pp = construct_vector(planepoint, start);
+  const Vector_3 end_min_pp = ray.direction().vector();
+  const Vector_3 normal = plane.orthogonal_vector();
+  const RT sdm_rs2pp = wdot(normal, start_min_pp, k);
+  const RT sdm_re2pp = wdot(normal, end_min_pp, k);
 
   switch (CGAL_NTS sign(sdm_rs2pp))
   {
     case -1:
-      if (sdm_re2pp > RT(0))
+      if(sdm_re2pp > RT(0))
         return FT(0);
       return squared_distance_to_plane(normal, start_min_pp, k);
     case 0:
     default:
       return FT(0);
     case 1:
-      if (sdm_re2pp < RT(0))
+      if(sdm_re2pp < RT(0))
         return FT(0);
       return squared_distance_to_plane(normal, start_min_pp, k);
   }
@@ -80,7 +80,7 @@ typename K::FT
 squared_distance(const Ray_3<K>& ray,
                  const Plane_3<K>& plane)
 {
-  return internal::squared_distance(ray, plane, K());
+  return K().compute_squared_distance_3_object()(ray, plane);
 }
 
 template <class K>
@@ -89,7 +89,7 @@ typename K::FT
 squared_distance(const Plane_3<K>& plane,
                  const Ray_3<K>& ray)
 {
-  return internal::squared_distance(ray, plane, K());
+  return K().compute_squared_distance_3_object()(plane, ray);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Ray_3_Ray_3.h b/Distance_3/include/CGAL/Distance_3/Ray_3_Ray_3.h
index 015e9a734c1..8ca1eca308c 100644
--- a/Distance_3/include/CGAL/Distance_3/Ray_3_Ray_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Ray_3_Ray_3.h
@@ -19,7 +19,7 @@
 
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
 
-#include <CGAL/Plane_3.h>
+#include <CGAL/Ray_3.h>
 #include <CGAL/Vector_3.h>
 
 namespace CGAL {
@@ -32,8 +32,8 @@ ray_ray_squared_distance_parallel(const typename K::Vector_3& ray1dir,
                                   const typename K::Vector_3& s1_min_s2,
                                   const K& k)
 {
-  if (!is_acute_angle(ray2dir, s1_min_s2, k))
-    if (!same_direction(ray1dir, ray2dir, k))
+  if(!is_acute_angle(ray2dir, s1_min_s2, k))
+    if(!same_direction(ray1dir, ray2dir, k))
       return typename K::FT(s1_min_s2*s1_min_s2);
 
   return squared_distance_to_line(ray1dir, s1_min_s2, k);
@@ -45,62 +45,73 @@ squared_distance(const typename K::Ray_3& ray1,
                  const typename K::Ray_3& ray2,
                  const K& k)
 {
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::Point_3 Point_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   const Point_3& s1 = ray1.source();
   const Point_3& s2 = ray2.source();
-  Vector_3 dir1, dir2, normal;
-  dir1 = ray1.direction().vector();
-  dir2 = ray2.direction().vector();
-  normal = wcross(dir1, dir2, k);
-  Vector_3 s1_min_s2 = construct_vector(s2, s1);
-  if (is_null(normal, k))
+  const Vector_3 dir1 = ray1.direction().vector();
+  const Vector_3 dir2 = ray2.direction().vector();
+  const Vector_3 normal = wcross(dir1, dir2, k);
+  const Vector_3 s1_min_s2 = construct_vector(s2, s1);
+
+  if(is_null(normal, k))
     return ray_ray_squared_distance_parallel(dir1, dir2, s1_min_s2, k);
 
   bool crossing1, crossing2;
-  RT sdm_s1_2, sdm_s2_1;
-  Vector_3 perpend1, perpend2;
-  perpend1 = wcross(dir1, normal, k);
-  perpend2 = wcross(dir2, normal, k);
 
-  sdm_s1_2 = wdot(perpend2, s1_min_s2, k);
-  if (sdm_s1_2 < RT(0)) {
-    crossing1 = (RT(wdot(perpend2, dir1, k)) >= RT(0));
-  } else {
-    if (RT(wdot(perpend2, dir1, k)) <= RT(0)) {
+  const Vector_3 perpend1 = wcross(dir1, normal, k);
+  const Vector_3 perpend2 = wcross(dir2, normal, k);
+
+  const RT sdm_s1_2 = wdot(perpend2, s1_min_s2, k);
+  if(sdm_s1_2 < RT(0))
+  {
+    crossing1 = (wdot(perpend2, dir1, k) >= RT(0));
+  }
+  else
+  {
+    if(RT(wdot(perpend2, dir1, k)) <= RT(0))
       crossing1 = true;
-    } else {
+    else
       crossing1 = (sdm_s1_2 == RT(0));
-    }
   }
 
-  sdm_s2_1 = -RT(wdot(perpend1, s1_min_s2, k));
-  if (sdm_s2_1 < RT(0)) {
-    crossing2 = (RT(wdot(perpend1, dir2, k)) >= RT(0));
-  } else {
-    if (RT(wdot(perpend1, dir2, k)) <= RT(0)) {
+  const RT sdm_s2_1 = - wdot(perpend1, s1_min_s2, k);
+  if(sdm_s2_1 < RT(0))
+  {
+    crossing2 = (wdot(perpend1, dir2, k) >= RT(0));
+  }
+  else
+  {
+    if(wdot(perpend1, dir2, k) <= RT(0))
       crossing2 = true;
-    } else {
+    else
       crossing2 = (sdm_s2_1 == RT(0));
-    }
   }
 
-  if (crossing1) {
-    if (crossing2)
+  if(crossing1)
+  {
+    if(crossing2)
       return squared_distance_to_plane(normal, s1_min_s2, k);
-    return squared_distance(ray2.source(), ray1, k);
-  } else {
-    if (crossing2) {
-      return squared_distance(ray1.source(), ray2, k);
-    } else {
+
+    return sq_dist(s2, ray1);
+  }
+  else
+  {
+    if(crossing2)
+    {
+      return sq_dist(s1, ray2);
+    }
+    else
+    {
       FT min1, min2;
-      min1 = squared_distance(ray1.source(), ray2, k);
-      min2 = squared_distance(ray2.source(), ray1, k);
+      min1 = sq_dist(s1, ray2);
+      min2 = sq_dist(s2, ray1);
       return (min1 < min2) ? min1 : min2;
     }
   }
@@ -124,7 +135,7 @@ typename K::FT
 squared_distance(const Ray_3<K>& ray1,
                  const Ray_3<K>& ray2)
 {
-  return internal::squared_distance(ray1, ray2, K());
+  return K().compute_squared_distance_3_object()(ray1, ray2);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h b/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h
index a6717e4e853..5781c5ceb74 100644
--- a/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Segment_3_Line_3.h
@@ -31,51 +31,50 @@ squared_distance(const typename K::Segment_3& seg,
                  const typename K::Line_3& line,
                  const K& k)
 {
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::Point_3 Point_3;
   typedef typename K::RT RT;
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
 
   typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   const Point_3& linepoint = line.point();
   const Point_3& start = seg.source();
   const Point_3& end = seg.target();
 
-  if (start == end)
-    return squared_distance(start, line, k);
+  if(start == end)
+    return sq_dist(start, line);
 
-  Vector_3 linedir = line.direction().vector();
-  Vector_3 segdir = seg.direction().vector();
-  Vector_3 normal = wcross(segdir, linedir, k);
+  const Vector_3 linedir = line.direction().vector();
+  const Vector_3 segdir = seg.direction().vector();
+  const Vector_3 normal = wcross(segdir, linedir, k);
 
-  if (is_null(normal, k))
+  if(is_null(normal, k))
     return squared_distance_to_line(linedir, vector(linepoint,start), k);
 
   bool crossing;
-  RT sdm_ss2l, sdm_se2l;
-  Vector_3 perpend2line, start_min_lp, end_min_lp;
-  perpend2line = wcross(linedir, normal, k);
-  start_min_lp = vector(linepoint, start);
-  end_min_lp = vector(linepoint, end);
-  sdm_ss2l = wdot(perpend2line, start_min_lp, k);
-  sdm_se2l = wdot(perpend2line, end_min_lp, k);
 
-  if (sdm_ss2l < RT(0)) {
+  const Vector_3 perpend2line = wcross(linedir, normal, k);
+  const Vector_3 start_min_lp = vector(linepoint, start);
+  const Vector_3 end_min_lp = vector(linepoint, end);
+  const RT sdm_ss2l = wdot(perpend2line, start_min_lp, k);
+  const RT sdm_se2l = wdot(perpend2line, end_min_lp, k);
+
+  if(sdm_ss2l < RT(0)) {
     crossing = (sdm_se2l >= RT(0));
   } else {
-    if (sdm_se2l <= RT(0)) {
+    if(sdm_se2l <= RT(0)) {
       crossing = true;
     } else {
       crossing = (sdm_ss2l == RT(0));
     }
   }
 
-  if (crossing) {
+  if(crossing) {
     return squared_distance_to_plane(normal, start_min_lp, k);
   } else {
-    RT dm;
-    dm = distance_measure_sub(sdm_ss2l, sdm_se2l, start_min_lp, end_min_lp, k);
-    if (dm <= RT(0)) {
+    const RT dm = distance_measure_sub(sdm_ss2l, sdm_se2l, start_min_lp, end_min_lp, k);
+    if(dm <= RT(0)) {
       return squared_distance_to_line(linedir, start_min_lp, k);
     } else {
       return squared_distance_to_line(linedir, end_min_lp, k);
@@ -100,7 +99,7 @@ typename K::FT
 squared_distance(const Segment_3<K>& seg,
                  const Line_3<K>& line)
 {
-  return internal::squared_distance(seg, line, K());
+  return K().compute_squared_distance_3_object()(seg, line);
 }
 
 template <class K>
@@ -109,7 +108,7 @@ typename K::FT
 squared_distance(const Line_3<K>& line,
                  const Segment_3<K>& seg)
 {
-  return internal::squared_distance(seg, line, K());
+  return K().compute_squared_distance_3_object()(line, seg);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Segment_3_Plane_3.h b/Distance_3/include/CGAL/Distance_3/Segment_3_Plane_3.h
index d67c6c3e543..1c8d6799556 100644
--- a/Distance_3/include/CGAL/Distance_3/Segment_3_Plane_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Segment_3_Plane_3.h
@@ -31,24 +31,27 @@ squared_distance(const typename K::Segment_3 &seg,
                  const typename K::Plane_3 &plane,
                  const K& k)
 {
-  typedef typename K::Point_3 Point_3;
-  typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::Point_3 Point_3;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+
+  const Point_3& start = seg.start();
+  const Point_3& end = seg.end();
 
-  const Point_3 &start = seg.start();
-  const Point_3 &end = seg.end();
   if (start == end)
     return squared_distance(start, plane, k);
 
-  const Point_3 &planepoint = plane.point();
-  Vector_3 start_min_pp = construct_vector(planepoint, start);
-  Vector_3 end_min_pp = construct_vector(planepoint, end);
-  const Vector_3 &normal = plane.orthogonal_vector();
-  RT sdm_ss2pp = wdot(normal, start_min_pp, k);
-  RT sdm_se2pp = wdot(normal, end_min_pp, k);
+  const Point_3& planepoint = plane.point();
+  const Vector_3 start_min_pp = vector(planepoint, start);
+  const Vector_3 end_min_pp = vector(planepoint, end);
+  const Vector_3& normal = plane.orthogonal_vector();
+
+  const RT sdm_ss2pp = wdot(normal, start_min_pp, k);
+  const RT sdm_se2pp = wdot(normal, end_min_pp, k);
+
   switch (CGAL_NTS sign(sdm_ss2pp))
   {
     case -1:
@@ -64,7 +67,7 @@ squared_distance(const typename K::Segment_3 &seg,
     case 1:
       if (sdm_se2pp <= RT(0))
         return FT(0);
-      if (sdm_ss2pp  * end_min_pp.hw() <= sdm_se2pp * start_min_pp.hw())
+      if (sdm_ss2pp * end_min_pp.hw() <= sdm_se2pp * start_min_pp.hw())
         return squared_distance_to_plane(normal, start_min_pp, k);
       else
         return squared_distance_to_plane(normal, end_min_pp, k);
@@ -88,7 +91,7 @@ typename K::FT
 squared_distance(const Segment_3<K>& seg,
                  const Plane_3<K>& plane)
 {
-  return internal::squared_distance(seg, plane, K());
+  return K().compute_squared_distance_3_object()(seg, plane);
 }
 
 template <class K>
@@ -97,7 +100,7 @@ typename K::FT
 squared_distance(const Plane_3<K>& plane,
                  const Segment_3<K>& seg)
 {
-  return internal::squared_distance(seg, plane, K());
+  return K().compute_squared_distance_3_object()(plane, seg);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Segment_3_Ray_3.h b/Distance_3/include/CGAL/Distance_3/Segment_3_Ray_3.h
index aa3fb6c1ee7..5943630d76f 100644
--- a/Distance_3/include/CGAL/Distance_3/Segment_3_Ray_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Segment_3_Ray_3.h
@@ -18,10 +18,6 @@
 #define CGAL_DISTANCE_3_SEGMENT_3_RAY_3_H
 
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
-#include <CGAL/Distance_3/Point_3_Line_3.h>
-#include <CGAL/Distance_3/Point_3_Point_3.h>
-#include <CGAL/Distance_3/Point_3_Ray_3.h>
-#include <CGAL/Distance_3/Point_3_Segment_3.h>
 
 #include <CGAL/Segment_3.h>
 #include <CGAL/Ray_3.h>
@@ -37,21 +33,23 @@ squared_distance_parallel(const typename K::Segment_3& seg,
 {
   typedef typename K::Vector_3 Vector_3;
 
+  const Vector_3 dir1 = seg.direction().vector();
+  const Vector_3 dir2 = ray.direction().vector();
+
   bool same_direction;
-  const Vector_3 &dir1 = seg.direction().vector();
-  const Vector_3 &dir2 = ray.direction().vector();
-
-  if (CGAL_NTS abs(dir1.hx()) > CGAL_NTS abs(dir1.hy())) {
+  if(CGAL_NTS abs(dir1.hx()) > CGAL_NTS abs(dir1.hy()))
     same_direction = (CGAL_NTS sign(dir1.hx()) == CGAL_NTS sign(dir2.hx()));
-  } else {
+  else
     same_direction = (CGAL_NTS sign(dir1.hy()) == CGAL_NTS sign(dir2.hy()));
-  }
 
-  if (same_direction) {
-    if (!is_acute_angle(seg.source(), seg.target(), ray.source(), k))
+  if(same_direction)
+  {
+    if(!is_acute_angle(seg.source(), seg.target(), ray.source(), k))
       return squared_distance(seg.target(), ray.source(), k);
-  } else {
-    if (!is_acute_angle(seg.target(), seg.source(), ray.source(), k))
+  }
+  else
+  {
+    if(!is_acute_angle(seg.target(), seg.source(), ray.source(), k))
       return squared_distance(seg.source(), ray.source(), k);
   }
 
@@ -64,88 +62,96 @@ squared_distance(const typename K::Segment_3& seg,
                  const typename K::Ray_3& ray,
                  const K& k)
 {
-  typedef typename K::Point_3 Point_3;
-  typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
+  typedef typename K::Point_3 Point_3;
+  typedef typename K::Vector_3 Vector_3;
 
-  typename K::Construct_vector_3 construct_vector;
+  typename K::Construct_vector_3 vector = k.construct_vector_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   const Point_3& ss = seg.source();
   const Point_3& se = seg.target();
 
-  if (ss == se)
-    return squared_distance(ss, ray, k);
+  if(ss == se)
+    return sq_dist(ss, ray);
 
-  Vector_3 raydir, segdir, normal;
-  raydir = ray.direction().vector();
-  segdir = seg.direction().vector();
-  normal = wcross(segdir, raydir, k);
+  const Vector_3 raydir = ray.direction().vector();
+  const Vector_3 segdir = seg.direction().vector();
+  const Vector_3 normal = wcross(segdir, raydir, k);
 
-  if (is_null(normal, k))
+  if(is_null(normal, k))
     return squared_distance_parallel(seg, ray, k);
 
   bool crossing1, crossing2;
-  RT sdm_ss2r, sdm_se2r, sdm_rs2s, sdm_re2s;
-  Vector_3 perpend2seg, perpend2ray, ss_min_rs, se_min_rs;
-  perpend2seg = wcross(segdir, normal, k);
-  perpend2ray = wcross(raydir, normal, k);
-  ss_min_rs = construct_vector(ray.source(), ss);
-  se_min_rs = construct_vector(ray.source(), se);
-  sdm_ss2r = wdot(perpend2ray, ss_min_rs, k);
-  sdm_se2r = wdot(perpend2ray, se_min_rs, k);
 
-  if (sdm_ss2r < RT(0)) {
+  const Vector_3 perpend2seg = wcross(segdir, normal, k);
+  const Vector_3 perpend2ray = wcross(raydir, normal, k);
+  const Vector_3 ss_min_rs = vector(ray.source(), ss);
+  const Vector_3 se_min_rs = vector(ray.source(), se);
+  const RT sdm_ss2r = wdot(perpend2ray, ss_min_rs, k);
+  const RT sdm_se2r = wdot(perpend2ray, se_min_rs, k);
+
+  if(sdm_ss2r < RT(0))
+  {
     crossing1 = (sdm_se2r >= RT(0));
-  } else {
-    if (sdm_se2r <= RT(0)) {
+  }
+  else
+  {
+    if(sdm_se2r <= RT(0))
       crossing1 = true;
-    } else {
+    else
       crossing1 = (sdm_ss2r == RT(0));
-    }
   }
 
-  sdm_rs2s = -RT(wdot(perpend2seg, ss_min_rs, k));
-  sdm_re2s = wdot(perpend2seg, raydir, k);
-  if (sdm_rs2s < RT(0)) {
+  const RT sdm_rs2s = - wdot(perpend2seg, ss_min_rs, k);
+  const RT sdm_re2s = wdot(perpend2seg, raydir, k);
+  if(sdm_rs2s < RT(0))
+  {
     crossing2 = (sdm_re2s >= RT(0));
-  } else {
-    if (sdm_re2s <= RT(0)) {
+  } else
+  {
+    if(sdm_re2s <= RT(0))
       crossing2 = true;
-    } else {
+    else
       crossing2 = (sdm_rs2s == RT(0));
-    }
   }
 
-  if (crossing1) {
-    if (crossing2) {
+  if(crossing1)
+  {
+    if(crossing2)
       return squared_distance_to_plane(normal, ss_min_rs, k);
-    }
-    return squared_distance(ray.source(), seg, k);
-  } else {
-    if (crossing2) {
-      RT dm;
-      dm = distance_measure_sub(sdm_ss2r, sdm_se2r, ss_min_rs, se_min_rs, k);
-      if (dm < RT(0)) {
-        return squared_distance(ss, ray, k);
-      } else {
-        if (dm > RT(0)) {
-          return squared_distance(se, ray, k);
-        } else {
+
+    return sq_dist(ray.source(), seg);
+  }
+  else
+  {
+    if(crossing2)
+    {
+      const RT dm = distance_measure_sub(sdm_ss2r, sdm_se2r, ss_min_rs, se_min_rs, k);
+      if(dm < RT(0))
+      {
+        return sq_dist(ss, ray);
+      }
+      else
+      {
+        if(dm > RT(0))
+          return sq_dist(se, ray);
+        else
           // parallel, should not happen (no crossing)
           return squared_distance_parallel(seg, ray, k);
-        }
       }
-    } else {
-      FT min1, min2;
-      RT dm;
-      dm = distance_measure_sub(sdm_ss2r, sdm_se2r, ss_min_rs, se_min_rs, k);
-      if (dm == RT(0))
+    }
+    else
+    {
+      const RT dm = distance_measure_sub(sdm_ss2r, sdm_se2r, ss_min_rs, se_min_rs, k);
+      if(dm == RT(0))
         return squared_distance_parallel(seg, ray, k);
-      min1 = (dm < RT(0))
-             ? squared_distance(ss, ray, k)
-             : squared_distance(se, ray, k);
-      min2 = squared_distance(ray.source(), seg, k);
+
+      const FT min1 = (dm < RT(0)) ? sq_dist(ss, ray)
+                                   : sq_dist(se, ray);
+      const FT min2 = sq_dist(ray.source(), seg);
+
       return (min1 < min2) ? min1 : min2;
     }
   }
@@ -168,7 +174,7 @@ typename K::FT
 squared_distance_parallel(const Segment_3<K>& seg,
                           const Ray_3<K>& ray)
 {
-  return internal::squared_distance_parallel(ray,seg, K());
+  return internal::squared_distance_parallel(ray, seg, K());
 }
 
 template <class K>
@@ -177,7 +183,7 @@ typename K::FT
 squared_distance(const Segment_3<K>& seg,
                  const Ray_3<K>& ray)
 {
-  return internal::squared_distance(seg, ray, K());
+  return K().compute_squared_distance_3_object()(seg, ray);
 }
 
 template <class K>
@@ -186,7 +192,7 @@ typename K::FT
 squared_distance(const Ray_3<K>& ray,
                  const Segment_3<K>& seg)
 {
-  return internal::squared_distance(seg, ray, K());
+  return K().compute_squared_distance_3_object()(ray, seg);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h b/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h
index 8232a7647b7..ad933787204 100644
--- a/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Segment_3_Segment_3.h
@@ -18,7 +18,6 @@
 #define CGAL_DISTANCE_3_SEGMENT_3_SEGMENT_3_H
 
 #include <CGAL/Distance_3/internal/squared_distance_utils_3.h>
-#include <CGAL/Distance_3/Point_3_Point_3.h>
 
 #include <CGAL/Segment_3.h>
 
@@ -49,6 +48,7 @@ squared_distance(const typename K::Segment_3& s1,
   typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object();
   typename K::Construct_vector_3 cv = k.construct_vector_3_object();
   typename K::Compute_scalar_product_3 sp = k.compute_scalar_product_3_object();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
 
   Segment_3_Segment_3_Result<K> res;
 
@@ -73,7 +73,7 @@ squared_distance(const typename K::Segment_3& s1,
     {
       res.x = 0;
       res.y = 0;
-      res.squared_distance = CGAL::squared_distance(p1, p2);
+      res.squared_distance = sq_dist(p1, p2);
       return res;
     }
 
@@ -81,7 +81,7 @@ squared_distance(const typename K::Segment_3& s1,
 
     res.x = 0;
     res.y = boost::algorithm::clamp<FT>(f/d, 0, 1); // (f - x*c) / d
-    res.squared_distance = CGAL::squared_distance(p1, p2 + res.y*v2);
+    res.squared_distance = sq_dist(p1, p2 + res.y*v2);
 
     return res;
   }
@@ -91,7 +91,7 @@ squared_distance(const typename K::Segment_3& s1,
 
     res.y = 0;
     res.x = boost::algorithm::clamp<FT>(e/a, 0, 1); // (e + y*c) / a
-    res.squared_distance = CGAL::squared_distance(p1 + res.x*v1, p2);
+    res.squared_distance = sq_dist(p1 + res.x*v1, p2);
 
     return res;
   }
@@ -128,7 +128,7 @@ squared_distance(const typename K::Segment_3& s1,
   CGAL_postcondition(FT(0) <= res.x && res.x <= FT(1));
   CGAL_postcondition(FT(0) <= res.y && res.y <= FT(1));
 
-  res.squared_distance = CGAL::squared_distance(p1 + res.x*v1, p2 + res.y*v2);
+  res.squared_distance = sq_dist(p1 + res.x*v1, p2 + res.y*v2);
 
   CGAL_postcondition(res.squared_distance >= FT(0));
 
@@ -148,24 +148,27 @@ squared_distance_parallel(const typename K::Segment_3& seg1,
 {
   typedef typename K::Vector_3 Vector_3;
 
-  const Vector_3& dir1 = seg1.direction().vector();
-  const Vector_3& dir2 = seg2.direction().vector();
+  typename K::Compute_squared_distance_3 sq_dist = k.compute_squared_distance_3_object();
+
+  const Vector_3 dir1 = seg1.direction().vector();
+  const Vector_3 dir2 = seg2.direction().vector();
 
   if(same_direction(dir1, dir2, k))
   {
     if(!is_acute_angle(seg1.source(), seg1.target(), seg2.source(), k))
-      return squared_distance(seg1.target(), seg2.source(), k);
+      return sq_dist(seg1.target(), seg2.source(), k);
     if(!is_acute_angle(seg1.target(), seg1.source(), seg2.target(), k))
-      return squared_distance(seg1.source(), seg2.target(), k);
+      return sq_dist(seg1.source(), seg2.target(), k);
   }
   else
   {
     if(!is_acute_angle(seg1.source(), seg1.target(), seg2.target(), k))
-      return squared_distance(seg1.target(), seg2.target(), k);
+      return sq_dist(seg1.target(), seg2.target(), k);
     if(!is_acute_angle(seg1.target(), seg1.source(), seg2.source(), k))
-      return squared_distance(seg1.source(), seg2.source(), k);
+      return sq_dist(seg1.source(), seg2.source(), k);
   }
-  return squared_distance(seg2.source(), seg1.supporting_line(), k);
+
+  return sq_dist(seg2.source(), seg1.supporting_line(), k);
 }
 
 template <typename K>
@@ -197,7 +200,7 @@ typename K::FT
 squared_distance(const Segment_3<K>& seg1,
                  const Segment_3<K>& seg2)
 {
-  return internal::squared_distance(seg1, seg2, K());
+  return K().compute_squared_distance_3_object()(seg1, seg2);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h b/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h
index 5454d73d67f..73f4b0eb69c 100644
--- a/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Triangle_3_Triangle_3.h
@@ -161,7 +161,7 @@ squared_distance(const typename K::Triangle_3& tr1,
   typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object();
 
   // ideally just limits<FT>::infinity|max(), but it is not available for exact NTs...
-  FT global_min_sqd = CGAL::squared_distance(vertex(tr1, 0), vertex(tr2, 0));
+  FT global_min_sqd = squared_distance(vertex(tr1, 0), vertex(tr2, 0));
 
   bool are_triangles_known_to_be_disjoint = false;
   std::pair<Distance_3::internal::Segment_3_Segment_3_Result<K>, bool> ss_res;
@@ -221,7 +221,7 @@ typename K::FT
 squared_distance(const Triangle_3<K>& tr1,
                  const Triangle_3<K>& tr2)
 {
-  return internal::squared_distance(tr1, tr2, K());
+  return K().compute_squared_distance_3_object()(tr1, tr2);
 }
 
 } // namespace CGAL
diff --git a/Distance_3/include/CGAL/Distance_3/Weighted_point_3_Weighted_point_3.h b/Distance_3/include/CGAL/Distance_3/Weighted_point_3_Weighted_point_3.h
index bd880062faf..54a41981103 100644
--- a/Distance_3/include/CGAL/Distance_3/Weighted_point_3_Weighted_point_3.h
+++ b/Distance_3/include/CGAL/Distance_3/Weighted_point_3_Weighted_point_3.h
@@ -17,9 +17,7 @@
 #ifndef CGAL_DISTANCE_3_WEIGHTED_POINT_3_WEIGHTED_POINT_3_H
 #define CGAL_DISTANCE_3_WEIGHTED_POINT_3_WEIGHTED_POINT_3_H
 
-#include <CGAL/Distance_3/Point_3_Point_3.h>
-
-#include <CGAL/Point_3.h>
+#include <CGAL/Weighted_point_3.h>
 
 namespace CGAL {
 
@@ -29,7 +27,7 @@ typename K::FT
 squared_distance(const Weighted_point_3<K>& wpt1,
                  const Weighted_point_3<K>& wpt2)
 {
-  return internal::squared_distance(wpt1.point(), wpt2.point(), K());
+  return K().compute_squared_distance_3_object()(wpt1.point(), wpt2.point());
 }
 
 } // namespace CGAL
diff --git a/Distance_3/test/Distance_3/test_distance_3.cpp b/Distance_3/test/Distance_3/test_distance_3.cpp
index 4354d158286..ff375560b55 100644
--- a/Distance_3/test/Distance_3/test_distance_3.cpp
+++ b/Distance_3/test/Distance_3/test_distance_3.cpp
@@ -2,6 +2,7 @@
 #include <CGAL/Simple_homogeneous.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Homogeneous.h>
 
 #include <CGAL/squared_distance_3.h>
 
@@ -252,9 +253,9 @@ private:
       P p2 = random_point();
       P p3 = random_point();
       p0 = CGAL::midpoint(p0, p1);
-      p1 = p0 + 0.1 * (p1 - p0);
-      p2 = p2 + (p2 - CGAL::ORIGIN) / CGAL::approximate_sqrt(CGAL::square(p2.x()) + CGAL::square(p2.y()) + CGAL::square(p2.z()) + 3);
-      p3 = p3 + (p3 - CGAL::ORIGIN) * std::cos(1.3);
+      p1 = p0 + FT(0.1) * V{p1 - p0};
+      p2 = p2 + V{p2 - CGAL::ORIGIN} / CGAL::approximate_sqrt(CGAL::square(p2.x()) + CGAL::square(p2.y()) + CGAL::square(p2.z()) + 3);
+      p3 = p3 + V{p3 - CGAL::ORIGIN} * FT(std::cos(1.3));
 
       // degenerate inputs
       check_squared_distance(S{p0, p0}, S{p0, p0}, 0); // both degen
@@ -273,9 +274,9 @@ private:
 
       // collinear segments
       const double lambda_4 = r.get_double(0, 1);
-      const P p4 = p2 + lambda_4 * (p3 - p2);
+      const P p4 = p2 + FT(lambda_4) * V{p3 - p2};
       const double lambda_5 = r.get_double(0, 1);
-      const P p5 = p2 + lambda_5 * (p3 - p2);
+      const P p5 = p2 + FT(lambda_5) * V{p3 - p2};
 
       // [p2;p3) fully contains [p4;p5]
       check_squared_distance(S{p2, p3}, S{p4, p5}, 0);
@@ -284,8 +285,7 @@ private:
       check_squared_distance(S{p3, p2}, S{p4, p5}, 0);
 
       const double lambda_6 = r.get_double(0, 1);
-      const P p6 = p3 + lambda_6 * (p3 - p2);
-
+      const P p6 = p3 + FT(lambda_6) * V{p3 - p2};
       // [p2;p3] overlaps [p5;p6]
       check_squared_distance(S{p2, p3}, S{p6, p5}, 0);
       check_squared_distance(S{p2, p3}, S{p5, p6}, 0);
@@ -293,7 +293,7 @@ private:
       check_squared_distance(S{p3, p2}, S{p6, p5}, 0);
 
       const double lambda_7 = r.get_double(1, 2);
-      const P p7 = p3 + lambda_7 * (p3 - p2);
+      const P p7 = p3 + FT(lambda_7) * V{p3 - p2};
 
       // [p2;p3] disjoint && left of [p6;p7]
       check_squared_distance(S{p2, p3}, S{p6, p7}, CGAL::squared_distance(p3, p6));
@@ -303,9 +303,9 @@ private:
 
       // Generic collinear
       const double lambda_8 = r.get_double(-M, M);
-      const P p8 = p2 + lambda_8 * (p3 - p2);
+      const P p8 = p2 + FT(lambda_8) * V{p3 - p2};
       const double lambda_9 = r.get_double(-M, M);
-      const P p9 = p2 + lambda_9 * (p3 - p2);
+      const P p9 = p2 + FT(lambda_9) * V{p3 - p2};
 
       S s89(p8, p9);
       S s32(p3, p2);
@@ -364,7 +364,7 @@ private:
 
       FT upper_bound = CGAL::squared_distance(p0, p2);
 
-      V p01 = V{p0, p1} / N;
+      V p01 = V{p0, p1} / FT(N);
       for(int l=0; l<N; ++l)
       {
         P tp = p0 + FT(l) * p01;
@@ -505,13 +505,13 @@ private:
 
       // these are still exact with EPECK
       p0 = CGAL::midpoint(p0, p1);
-      p1 = p0 + 0.1 * (p1 - p0);
-      p2 = p2 + (p2 - p0) / CGAL_PI;
+      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 + (p4 - CGAL::ORIGIN) / CGAL::approximate_sqrt(CGAL::square(p4.x()) + CGAL::square(p4.y()) + CGAL::square(p4.z()) + 3);
+      p4 = p4 + V{p4 - CGAL::ORIGIN} / CGAL::approximate_sqrt(CGAL::square(p4.x()) + CGAL::square(p4.y()) + CGAL::square(p4.z()) + 3);
 
-      p5 = p5 + (p5 - CGAL::ORIGIN) * std::cos(1.3);
+      p5 = p5 + V{p5 - CGAL::ORIGIN} * FT(std::cos(1.3));
 
       // degenerate inputs
       check_squared_distance(T{p3, p3, p3}, T{p3, p3, p3}, 0); // both degen
@@ -670,6 +670,8 @@ int main()
 //  Test<CGAL::Simple_homogeneous<double> >(r, 1e-5).run();
 //  Test<CGAL::Simple_cartesian<CGAL::Interval_nt<true> > >(r, 1e-5).run();
 
+  Test<CGAL::Homogeneous<CGAL::Epeck_ft> >(r, 0).run();
+
   const double epick_eps = 10 * std::numeric_limits<double>::epsilon();
   Test<CGAL::Exact_predicates_inexact_constructions_kernel>(r, epick_eps).run();