From 28a6a91453447d4b6c060b7cc07ff326d5acba87 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Wed, 5 May 2021 16:54:05 +0200
Subject: [PATCH 01/11] Intersections_3: Make do_intersect(sphere/triangle)
 work if FT is division free

---
 .../include/CGAL/squared_distance_3_0.h       |  65 +++++++
 .../include/CGAL/squared_distance_3_1.h       | 164 ++++++++++++++++++
 .../include/CGAL/squared_distance_3_2.h       |  80 +++++++++
 .../Triangle_3_Sphere_3_do_intersect.h        |  74 +++++---
 .../Intersections_3/test_intersections_3.cpp  |  27 ++-
 5 files changed, 384 insertions(+), 26 deletions(-)

diff --git a/Distance_3/include/CGAL/squared_distance_3_0.h b/Distance_3/include/CGAL/squared_distance_3_0.h
index 7b496c5967a..64d3a1af570 100644
--- a/Distance_3/include/CGAL/squared_distance_3_0.h
+++ b/Distance_3/include/CGAL/squared_distance_3_0.h
@@ -202,6 +202,23 @@ squared_distance_to_plane(const typename K::Vector_3 & normal,
   return FT(dot*dot) /
     FT(wmult((K*)0, squared_length, diff.hw(), diff.hw()));
 }
+  
+template <class K>
+void
+squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
+                             const typename K::Vector_3 & diff,
+                             typename K::RT& num,
+                             typename K::RT& den,
+                             const K& k)
+{
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT dot, squared_length;
+  dot = wdot(normal, diff, k);
+  squared_length = wdot(normal, normal, k);
+  num =  square(dot);
+  den = wmult((K*)0, squared_length, diff.hw(), diff.hw());
+}
 
 
 template <class K>
@@ -218,7 +235,55 @@ squared_distance_to_line(const typename K::Vector_3 & dir,
                               (K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw()));
 }
 
+template <class K>
+void
+squared_distance_to_line_RT(const typename K::Vector_3 & dir,
+                            const typename K::Vector_3 & diff,
+                            typename K::RT& num,
+                            typename K::RT& den,
+                            const K& k,
+                            Cartesian_tag)
+{
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  Vector_3 wcr = wcross(dir, diff, k);
+  num = wcr*wcr;
+  den = wmult((K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw());
+}
 
+template <class K>
+void
+squared_distance_to_line_RT(const typename K::Vector_3 & dir,
+                            const typename K::Vector_3 & diff,
+                            typename K::RT& num,
+                            typename K::RT& den,
+                            const K& k,
+                            Homogeneous_tag)
+{
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  Vector_3 wcr = wcross(dir, diff, k);
+  FT ft = FT(wcr*wcr)/FT(wmult(
+                               (K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw()));
+  num = Rational_traits<FT>().numerator(ft);
+  den = Rational_traits<FT>().denominator(ft);
+}
+
+template <class K>
+void
+squared_distance_to_line_RT(const typename K::Vector_3 & dir,
+                            const typename K::Vector_3 & diff,
+                            typename K::RT& num,
+                            typename K::RT& den,
+                            const K& k)
+{
+  typedef typename K::Kernel_tag Tag;
+  Tag tag;
+  squared_distance_to_line_RT(dir,diff,num,den,k,tag);
+}
+                                       
 template <class K>
 inline
 bool
diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h
index e0e63a6fb23..ba2b7d435d4 100644
--- a/Distance_3/include/CGAL/squared_distance_3_1.h
+++ b/Distance_3/include/CGAL/squared_distance_3_1.h
@@ -43,6 +43,22 @@ squared_distance(
   Vector_3 diff = construct_vector(line.point(), pt);
   return internal::squared_distance_to_line(dir, diff, k);
 }
+  
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Line_3 &line,
+    typename K::RT num,
+    typename K::RT den,
+    const K& k)
+{
+  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);
+  return internal::squared_distance_to_line_RT(dir, diff, num, den, k);
+}
 
 template <class K>
 inline
@@ -73,6 +89,70 @@ squared_distance(
     return squared_distance_to_line(dir, diff, k);
 }
 
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Ray_3 &ray,
+    typename K::RT num,
+    typename K::RT den,
+    const K& k,
+    Homogeneous_tag&)
+{
+  typename K::Construct_vector_3 construct_vector;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  
+    Vector_3 diff = construct_vector(ray.source(), pt);
+    const Vector_3 &dir = ray.direction().vector();
+    if (!is_acute_angle(dir,diff, k) ){
+      FT ft = diff*diff;
+      num = Rational_traits<FT>().numerator(ft);
+      den = Rational_traits<FT>().denominator(ft);
+      return;
+    }
+    squared_distance_to_line_RT(dir, diff, num, den, k);
+}
+
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Ray_3 &ray,
+    typename K::RT num,
+    typename K::RT den,
+    const K& k,
+    Cartesian_tag& )
+{
+  typename K::Construct_vector_3 construct_vector;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+  
+    Vector_3 diff = construct_vector(ray.source(), pt);
+    const Vector_3 &dir = ray.direction().vector();
+    if (!is_acute_angle(dir,diff, k) ){
+         num = diff*diff;
+         den =  RT(1);
+         return;
+    }
+    squared_distance_to_line_RT(dir, diff, num, den, k);
+}
+
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Ray_3 &ray,
+    typename K::RT num,
+    typename K::RT den,
+    const K& k)
+{
+  typedef typename K::Kernel_tag Tag;
+  Tag tag;
+  squared_distance_RT(pt,ray,num,den,k,tag);
+}
+  
 
 template <class K>
 inline
@@ -114,6 +194,42 @@ squared_distance(
     return FT(wcr*wcr)/FT(e * diff.hw() * diff.hw());
 }
 
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Segment_3 &seg,
+    typename K::RT& num,
+    typename K::RT& den,
+    const K& k,
+    const Homogeneous_tag)
+{
+    typename K::Construct_vector_3 construct_vector;
+    typedef typename K::Vector_3 Vector_3;
+    typedef typename K::RT RT;
+    typedef typename K::FT FT;
+    // 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());
+    RT d = wdot(diff,segvec, k);
+    if (d <= (RT)0){
+      FT ft = FT(diff*diff);
+      num = Rational_traits<FT>().numerator(ft);
+      den = Rational_traits<FT>().denominator(ft);
+    }
+    RT e = wdot(segvec,segvec, k);
+    if ( (d * segvec.hw()) > (e * diff.hw())){
+      FT sd = squared_distance(pt, seg.target(), k);
+      num = Rational_traits<FT>().numerator(sd);
+      den = Rational_traits<FT>().denominator(sd);
+      return;
+    }
+    Vector_3 wcr = wcross(segvec, diff, k);
+    FT ft((wcr*wcr)/FT(e * diff.hw() * diff.hw()));
+    num = Rational_traits<FT>().numerator(ft);
+    den = Rational_traits<FT>().denominator(ft);
+}
+
 template <class K>
 typename K::FT
 squared_distance(
@@ -140,6 +256,54 @@ squared_distance(
     return FT(wcr*wcr)/e;
 }
 
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Segment_3 &seg,
+    typename K::RT& num,
+    typename K::RT& den,
+    const K& k,
+    const Cartesian_tag&)
+{
+    typename K::Construct_vector_3 construct_vector;
+    typedef typename K::Vector_3 Vector_3;
+    typedef typename K::RT RT;
+    typedef typename K::FT FT;
+    // 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());
+    RT d = wdot(diff,segvec, k);
+    if (d <= (RT)0){
+      num = diff*diff;
+      den = RT(1);
+      return;
+    }
+    RT e = wdot(segvec,segvec, k);
+    if (d > e){
+      num = squared_distance(pt, seg.target(), k);
+      den = RT(1);
+      return;
+    }
+    Vector_3 wcr = wcross(segvec, diff, k);
+    num = wcr*wcr;
+    den = e;
+}
+
+template <class K>
+inline
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Segment_3 &seg,
+    typename K::RT& num,
+    typename K::RT& den,
+    const K& k)
+{
+  typedef typename K::Kernel_tag Tag;
+  Tag tag;
+  return squared_distance_RT(pt, seg, num, den, k, tag);
+}
 
 template <class K>
 inline
diff --git a/Distance_3/include/CGAL/squared_distance_3_2.h b/Distance_3/include/CGAL/squared_distance_3_2.h
index 78af20db08f..785d9365722 100644
--- a/Distance_3/include/CGAL/squared_distance_3_2.h
+++ b/Distance_3/include/CGAL/squared_distance_3_2.h
@@ -261,6 +261,65 @@ squared_distance_to_triangle(
       }
 }
 
+template <class K>
+inline void
+squared_distance_to_triangle_RT(
+    const typename K::Point_3 & pt,
+    const typename K::Point_3 & t0,
+    const typename K::Point_3 & t1,
+    const typename K::Point_3 & t2,
+    bool & inside,
+    typename K::RT& num,
+    typename K::RT& den,
+    const K& k)
+{
+  typename K::Construct_vector_3 vector;
+  typedef typename K::Vector_3 Vector_3;
+  const Vector_3 e1 = vector(t0, t1);
+  const Vector_3 oe3 = vector(t0, t2);
+  const Vector_3 normal = wcross(e1, oe3, k);
+
+  if(normal != NULL_VECTOR
+     && on_left_of_triangle_edge(pt, normal, t0, t1, k)
+     && on_left_of_triangle_edge(pt, normal, t1, t2, k)
+     && on_left_of_triangle_edge(pt, normal, t2, t0, k))
+      {
+        // the projection of pt is inside the triangle
+        inside = true;
+        return squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
+      }
+      else {
+        // The case normal==NULL_VECTOR covers the case when the triangle
+        // is colinear, or even more degenerate. In that case, we can
+        // simply take also the distance to the three segments.
+
+        squared_distance_RT(pt,
+                           typename K::Segment_3(t2, t0),
+                           num,
+                           den,
+                           k);
+        typename K::RT num2, den2;
+        squared_distance_RT(pt,
+                            typename K::Segment_3(t1, t2),
+                            num2,
+                            den2,
+                            k);
+        if(compare_quotients(num2,den2,num,den) == SMALLER){
+          num =  num2;
+          den = den2;
+        }
+        squared_distance_RT(pt,
+                            typename K::Segment_3(t0, t1),
+                            num2,
+                            den2,
+                            k);
+        if(compare_quotients(num2,den2,num,den) == SMALLER){
+          num =  num2;
+          den = den2;
+        }
+      }
+}
+
 template <class K>
 inline typename K::FT
 squared_distance(
@@ -278,6 +337,27 @@ squared_distance(
                                       k);
 }
 
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 & pt,
+    const typename K::Triangle_3 & t,
+    typename K::RT& num,
+    typename K::RT& den,
+    const K& k)
+{
+  typename K::Construct_vertex_3 vertex;
+  bool inside = false;
+  squared_distance_to_triangle_RT(pt,
+                                  vertex(t, 0),
+                                  vertex(t, 1),
+                                  vertex(t, 2),
+                                  inside,
+                                  num,
+                                  den,
+                                  k);
+}
+
 } // namespace internal
 
 
diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
index f8fa015669a..a7aac04faed 100644
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
@@ -15,7 +15,7 @@
 
 #include <CGAL/squared_distance_3_2.h>
 #include <CGAL/Intersection_traits_3.h>
-
+#include <CGAL/Rational_traits.h>
 namespace CGAL {
 
 template <class K>
@@ -36,9 +36,15 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Sphere_3 &sp,
              const typename K::Triangle_3 &tr,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), tr) <= sp.squared_radius();
+  typedef typename K::RT RT;
+  RT num, den;
+  
+  CGAL::internal::squared_distance_RT(sp.center(), tr, num, den, k);
+  return ! (compare_quotients<RT>(num, den,
+                                  Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
+                                  Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
 }
 
 template <class K>
@@ -46,18 +52,9 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Triangle_3 &tr,
              const typename K::Sphere_3 &sp,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), tr) <= sp.squared_radius();
-}
-template <class K>
-inline
-typename K::Boolean
-do_intersect(const typename K::Sphere_3 &sp,
-             const typename K::Line_3 &lin,
-             const K & /* k */)
-{
-  return squared_distance(sp.center(), lin) <= sp.squared_radius();
+  return do_intersect(sp, tr, k); 
 }
 
 
@@ -66,9 +63,26 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Line_3 &lin,
              const typename K::Sphere_3 &sp,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), lin) <= sp.squared_radius();
+  typedef typename K::RT RT;
+  RT num, den;
+  
+  CGAL::internal::squared_distance_RT(sp.center(), lin, num, den, k);
+  return ! (compare_quotients<RT>(num, den,
+                                  Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
+                                  Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
+
+}
+  
+template <class K>
+inline
+typename K::Boolean
+do_intersect(const typename K::Sphere_3 &sp,
+             const typename K::Line_3 &lin,
+             const K & k)
+{
+  return do_intersect(lin,sp,k);
 }
 
 
@@ -78,9 +92,15 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Sphere_3 &sp,
              const typename K::Ray_3 &ray,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), ray) <= sp.squared_radius();
+  typedef typename K::RT RT;
+  RT num, den;
+  
+  CGAL::internal::squared_distance_RT(sp.center(), ray, num, den, k);
+  return ! (compare_quotients<RT>(num, den,
+                                  Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
+                                  Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
 }
 
 
@@ -89,9 +109,9 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Ray_3 &ray,
              const typename K::Sphere_3 &sp,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), ray) <= sp.squared_radius();
+  return do_intersect(sp,ray,k);
 }
 
 template <class K>
@@ -99,9 +119,15 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Sphere_3 &sp,
              const typename K::Segment_3 &seg,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), seg) <= sp.squared_radius();
+  typedef typename K::RT RT;
+  RT num, den;
+  
+  CGAL::internal::squared_distance_RT(sp.center(), seg, num, den, k);
+  return ! (compare_quotients<RT>(num, den,
+                                  Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
+                                  Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
 }
 
 
@@ -110,9 +136,9 @@ inline
 typename K::Boolean
 do_intersect(const typename K::Segment_3 &seg,
              const typename K::Sphere_3 &sp,
-             const K & /* k */)
+             const K & k)
 {
-  return squared_distance(sp.center(), seg) <= sp.squared_radius();
+  return do_intersect(sp,seg,k);
 }
 
 } // namespace internal
diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
index 65ec3326031..d9ce491bc0a 100644
--- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp
+++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
@@ -161,14 +161,36 @@ struct Test {
 
   void P_do_intersect()
   {
-    P p(0,0,0), q(1,0,0), r(2,0,0), s(10,10,10);
+    P p(0,0,0), q(1,0,0), r(2,0,0), s(10,10,10), s2(10,10,0), s3(10,-10,-1), t(3,0,0);
     Sph sph(p,1);
     Cub cub(p,r);
     assert(do_intersect(q,sph));
     assert(do_intersect(sph,q));
     assert(! do_intersect(s,cub));
     assert(! do_intersect(cub,s));
-  }
+
+    assert(do_intersect(sph, Tr(p,q,s)));
+    assert(do_intersect(sph, Tr(r,q,s)));
+    assert(! do_intersect(sph, Tr(r,s,t)));
+
+    assert(! do_intersect(sph, Tr(s,s2,s3)));
+    
+    assert(do_intersect(sph, S(p,r)));
+    assert(do_intersect(sph, S(q,r)));
+    assert(do_intersect(sph, S(s,q)));
+    assert(! do_intersect(sph, S(s,r)));
+    /*
+    assert(do_intersect(sph, L(p,r)));
+    assert(do_intersect(sph, L(q,r)));
+    assert(do_intersect(sph, L(s,q)));
+    assert(! do_intersect(sph, L(s,r)));
+  
+    assert(do_intersect(sph, R(p,r)));
+    assert(do_intersect(sph, R(q,r)));
+    assert(do_intersect(sph, R(s,q)));
+    assert(! do_intersect(sph, R(s,r)));
+    */
+  }  
 
 
   void Cub_Cub()
@@ -1265,6 +1287,7 @@ struct Test {
     Bbox_L();
     Bbox_R();
     Bbox_Tr();
+
   }
 };
 

From 10e1b51ed3ee76652f1e7e5ab5f39dc59aca6f00 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Wed, 5 May 2021 17:24:37 +0200
Subject: [PATCH 02/11] Add &s

---
 Distance_3/include/CGAL/squared_distance_3_1.h   | 16 ++++++++--------
 .../Intersections_3/test_intersections_3.cpp     |  4 ++--
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h
index ba2b7d435d4..39cfe3fd03d 100644
--- a/Distance_3/include/CGAL/squared_distance_3_1.h
+++ b/Distance_3/include/CGAL/squared_distance_3_1.h
@@ -49,8 +49,8 @@ void
 squared_distance_RT(
     const typename K::Point_3 &pt,
     const typename K::Line_3 &line,
-    typename K::RT num,
-    typename K::RT den,
+    typename K::RT& num,
+    typename K::RT& den,
     const K& k)
 {
   typedef typename K::Vector_3 Vector_3;
@@ -94,8 +94,8 @@ void
 squared_distance_RT(
     const typename K::Point_3 &pt,
     const typename K::Ray_3 &ray,
-    typename K::RT num,
-    typename K::RT den,
+    typename K::RT& num,
+    typename K::RT& den,
     const K& k,
     Homogeneous_tag&)
 {
@@ -120,8 +120,8 @@ void
 squared_distance_RT(
     const typename K::Point_3 &pt,
     const typename K::Ray_3 &ray,
-    typename K::RT num,
-    typename K::RT den,
+    typename K::RT& num,
+    typename K::RT& den,
     const K& k,
     Cartesian_tag& )
 {
@@ -144,8 +144,8 @@ void
 squared_distance_RT(
     const typename K::Point_3 &pt,
     const typename K::Ray_3 &ray,
-    typename K::RT num,
-    typename K::RT den,
+    typename K::RT& num,
+    typename K::RT& den,
     const K& k)
 {
   typedef typename K::Kernel_tag Tag;
diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
index d9ce491bc0a..e3c8b68fe28 100644
--- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp
+++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
@@ -179,12 +179,12 @@ struct Test {
     assert(do_intersect(sph, S(q,r)));
     assert(do_intersect(sph, S(s,q)));
     assert(! do_intersect(sph, S(s,r)));
-    /*
+    
     assert(do_intersect(sph, L(p,r)));
     assert(do_intersect(sph, L(q,r)));
     assert(do_intersect(sph, L(s,q)));
     assert(! do_intersect(sph, L(s,r)));
-  
+    /*
     assert(do_intersect(sph, R(p,r)));
     assert(do_intersect(sph, R(q,r)));
     assert(do_intersect(sph, R(s,q)));

From b7b2e1955556defdb60088db16ab7d7404c228c7 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Wed, 5 May 2021 17:39:59 +0200
Subject: [PATCH 03/11] Add the sphere ray tests

---
 Intersections_3/test/Intersections_3/test_intersections_3.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
index e3c8b68fe28..13fb6846007 100644
--- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp
+++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
@@ -184,12 +184,12 @@ struct Test {
     assert(do_intersect(sph, L(q,r)));
     assert(do_intersect(sph, L(s,q)));
     assert(! do_intersect(sph, L(s,r)));
-    /*
+    
     assert(do_intersect(sph, R(p,r)));
     assert(do_intersect(sph, R(q,r)));
     assert(do_intersect(sph, R(s,q)));
     assert(! do_intersect(sph, R(s,r)));
-    */
+    
   }  
 
 

From 020ec5408e2e2e4ca10044791e7a049f445e56da Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Thu, 6 May 2021 07:51:29 +0200
Subject: [PATCH 04/11] Fixes concerning return

---
 Distance_3/include/CGAL/squared_distance_3_1.h | 1 +
 Distance_3/include/CGAL/squared_distance_3_2.h | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h
index 39cfe3fd03d..ff8d87409e8 100644
--- a/Distance_3/include/CGAL/squared_distance_3_1.h
+++ b/Distance_3/include/CGAL/squared_distance_3_1.h
@@ -216,6 +216,7 @@ squared_distance_RT(
       FT ft = FT(diff*diff);
       num = Rational_traits<FT>().numerator(ft);
       den = Rational_traits<FT>().denominator(ft);
+      return;
     }
     RT e = wdot(segvec,segvec, k);
     if ( (d * segvec.hw()) > (e * diff.hw())){
diff --git a/Distance_3/include/CGAL/squared_distance_3_2.h b/Distance_3/include/CGAL/squared_distance_3_2.h
index 785d9365722..6a9acd00153 100644
--- a/Distance_3/include/CGAL/squared_distance_3_2.h
+++ b/Distance_3/include/CGAL/squared_distance_3_2.h
@@ -286,7 +286,7 @@ squared_distance_to_triangle_RT(
       {
         // the projection of pt is inside the triangle
         inside = true;
-        return squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
+        squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
       }
       else {
         // The case normal==NULL_VECTOR covers the case when the triangle

From f45e843cd0c5755c4dd1681b8a3186e1e75ba644 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Thu, 6 May 2021 08:03:53 +0200
Subject: [PATCH 05/11] trailing whitespace

---
 Distance_3/include/CGAL/squared_distance_3_0.h       |  4 ++--
 Distance_3/include/CGAL/squared_distance_3_1.h       |  8 ++++----
 .../internal/Triangle_3_Sphere_3_do_intersect.h      | 12 ++++++------
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/Distance_3/include/CGAL/squared_distance_3_0.h b/Distance_3/include/CGAL/squared_distance_3_0.h
index 64d3a1af570..7ec5e64337e 100644
--- a/Distance_3/include/CGAL/squared_distance_3_0.h
+++ b/Distance_3/include/CGAL/squared_distance_3_0.h
@@ -202,7 +202,7 @@ squared_distance_to_plane(const typename K::Vector_3 & normal,
   return FT(dot*dot) /
     FT(wmult((K*)0, squared_length, diff.hw(), diff.hw()));
 }
-  
+
 template <class K>
 void
 squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
@@ -283,7 +283,7 @@ squared_distance_to_line_RT(const typename K::Vector_3 & dir,
   Tag tag;
   squared_distance_to_line_RT(dir,diff,num,den,k,tag);
 }
-                                       
+
 template <class K>
 inline
 bool
diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h
index ff8d87409e8..1a883a35245 100644
--- a/Distance_3/include/CGAL/squared_distance_3_1.h
+++ b/Distance_3/include/CGAL/squared_distance_3_1.h
@@ -43,7 +43,7 @@ squared_distance(
   Vector_3 diff = construct_vector(line.point(), pt);
   return internal::squared_distance_to_line(dir, diff, k);
 }
-  
+
 template <class K>
 void
 squared_distance_RT(
@@ -103,7 +103,7 @@ squared_distance_RT(
   typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
-  
+
     Vector_3 diff = construct_vector(ray.source(), pt);
     const Vector_3 &dir = ray.direction().vector();
     if (!is_acute_angle(dir,diff, k) ){
@@ -128,7 +128,7 @@ squared_distance_RT(
   typename K::Construct_vector_3 construct_vector;
   typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
-  
+
     Vector_3 diff = construct_vector(ray.source(), pt);
     const Vector_3 &dir = ray.direction().vector();
     if (!is_acute_angle(dir,diff, k) ){
@@ -152,7 +152,7 @@ squared_distance_RT(
   Tag tag;
   squared_distance_RT(pt,ray,num,den,k,tag);
 }
-  
+
 
 template <class K>
 inline
diff --git a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
index a7aac04faed..e35346a9d5a 100644
--- a/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
+++ b/Intersections_3/include/CGAL/Intersections_3/internal/Triangle_3_Sphere_3_do_intersect.h
@@ -40,7 +40,7 @@ do_intersect(const typename K::Sphere_3 &sp,
 {
   typedef typename K::RT RT;
   RT num, den;
-  
+
   CGAL::internal::squared_distance_RT(sp.center(), tr, num, den, k);
   return ! (compare_quotients<RT>(num, den,
                                   Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
@@ -54,7 +54,7 @@ do_intersect(const typename K::Triangle_3 &tr,
              const typename K::Sphere_3 &sp,
              const K & k)
 {
-  return do_intersect(sp, tr, k); 
+  return do_intersect(sp, tr, k);
 }
 
 
@@ -67,14 +67,14 @@ do_intersect(const typename K::Line_3 &lin,
 {
   typedef typename K::RT RT;
   RT num, den;
-  
+
   CGAL::internal::squared_distance_RT(sp.center(), lin, num, den, k);
   return ! (compare_quotients<RT>(num, den,
                                   Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
                                   Rational_traits<typename K::FT>().denominator(sp.squared_radius())) == LARGER);
 
 }
-  
+
 template <class K>
 inline
 typename K::Boolean
@@ -96,7 +96,7 @@ do_intersect(const typename K::Sphere_3 &sp,
 {
   typedef typename K::RT RT;
   RT num, den;
-  
+
   CGAL::internal::squared_distance_RT(sp.center(), ray, num, den, k);
   return ! (compare_quotients<RT>(num, den,
                                   Rational_traits<typename K::FT>().numerator(sp.squared_radius()),
@@ -123,7 +123,7 @@ do_intersect(const typename K::Sphere_3 &sp,
 {
   typedef typename K::RT RT;
   RT num, den;
-  
+
   CGAL::internal::squared_distance_RT(sp.center(), seg, num, den, k);
   return ! (compare_quotients<RT>(num, den,
                                   Rational_traits<typename K::FT>().numerator(sp.squared_radius()),

From 4f7703fd94844d66e928454722d5a8db6f119c12 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Thu, 6 May 2021 08:07:39 +0200
Subject: [PATCH 06/11] trailing whitespace

---
 .../test/Intersections_3/test_intersections_3.cpp      | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/Intersections_3/test/Intersections_3/test_intersections_3.cpp b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
index 13fb6846007..191d5a0f7aa 100644
--- a/Intersections_3/test/Intersections_3/test_intersections_3.cpp
+++ b/Intersections_3/test/Intersections_3/test_intersections_3.cpp
@@ -174,23 +174,23 @@ struct Test {
     assert(! do_intersect(sph, Tr(r,s,t)));
 
     assert(! do_intersect(sph, Tr(s,s2,s3)));
-    
+
     assert(do_intersect(sph, S(p,r)));
     assert(do_intersect(sph, S(q,r)));
     assert(do_intersect(sph, S(s,q)));
     assert(! do_intersect(sph, S(s,r)));
-    
+
     assert(do_intersect(sph, L(p,r)));
     assert(do_intersect(sph, L(q,r)));
     assert(do_intersect(sph, L(s,q)));
     assert(! do_intersect(sph, L(s,r)));
-    
+
     assert(do_intersect(sph, R(p,r)));
     assert(do_intersect(sph, R(q,r)));
     assert(do_intersect(sph, R(s,q)));
     assert(! do_intersect(sph, R(s,r)));
-    
-  }  
+
+  }
 
 
   void Cub_Cub()

From 599e8a74d05f2e81d9c02c81fd7f0ffee7bea9b8 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Thu, 6 May 2021 09:02:51 +0200
Subject: [PATCH 07/11] include <CGAL/Rational_traits.h>

---
 Distance_3/include/CGAL/squared_distance_3_0.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Distance_3/include/CGAL/squared_distance_3_0.h b/Distance_3/include/CGAL/squared_distance_3_0.h
index 7ec5e64337e..fe4b49edbfd 100644
--- a/Distance_3/include/CGAL/squared_distance_3_0.h
+++ b/Distance_3/include/CGAL/squared_distance_3_0.h
@@ -26,6 +26,7 @@
 #include <CGAL/Weighted_point_3.h>
 #include <CGAL/Vector_3.h>
 #include <CGAL/number_utils.h>
+#include <CGAL/Rational_traits.h>
 
 namespace CGAL {
 

From c46bccf4ac21c66a82c4e0d131d992151b02a3f8 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Thu, 6 May 2021 12:43:31 +0200
Subject: [PATCH 08/11] Add point/triangle distance tests

---
 Distance_3/test/Distance_3/test_distance_3.cpp | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/Distance_3/test/Distance_3/test_distance_3.cpp b/Distance_3/test/Distance_3/test_distance_3.cpp
index 3a1f72fe786..0632b699765 100644
--- a/Distance_3/test/Distance_3/test_distance_3.cpp
+++ b/Distance_3/test/Distance_3/test_distance_3.cpp
@@ -112,6 +112,17 @@ struct Test {
   {
     std::cout << "Point - Triangle\n";
     check_squared_distance (p(0, 1, 2), T(p(0, 0, 0), p( 2, 0, 0), p( 0, 2, 0)), 4);
+
+    T t(p(0,0,0), p(3,0,0), p(3,3,0));
+    check_squared_distance (p(-1, -1, 0), t, 2);
+    check_squared_distance (p(-1,  0, 0), t, 1);
+    check_squared_distance (p(0, 0, 0), t, 0);
+    check_squared_distance (p(1, 0, 0), t, 0);
+    check_squared_distance (p(4, 0, 0), t, 1);
+    check_squared_distance (p(1, -1, 0), t, 1);
+    check_squared_distance (p(1, 1, 1), t, 1);
+    check_squared_distance (p(1, 0, 1), t, 1);
+    check_squared_distance (p(0, 0, 1), t, 1);
   }
 
   void P_Tet()

From f29e7ca937fe63e844db905316d14153f4f11cb6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Fri, 7 May 2021 12:08:17 +0200
Subject: [PATCH 09/11] Factorize RT/non_RT and Cartesian/Homogeneous
 implementations

Note that the function squared_distance_to_triangle is modified,
see PR https://github.com/CGAL/cgal/pull/5647
---
 .../include/CGAL/squared_distance_3_0.h       |  89 ++---
 .../include/CGAL/squared_distance_3_1.h       | 328 ++++++------------
 .../include/CGAL/squared_distance_3_2.h       | 215 ++++++------
 3 files changed, 251 insertions(+), 381 deletions(-)

diff --git a/Distance_3/include/CGAL/squared_distance_3_0.h b/Distance_3/include/CGAL/squared_distance_3_0.h
index fe4b49edbfd..835c2b97b66 100644
--- a/Distance_3/include/CGAL/squared_distance_3_0.h
+++ b/Distance_3/include/CGAL/squared_distance_3_0.h
@@ -188,22 +188,6 @@ squared_distance(const typename K::Point_3 & pt1,
   return k.compute_squared_distance_3_object()(pt1, pt2);
 }
 
-
-template <class K>
-typename K::FT
-squared_distance_to_plane(const typename K::Vector_3 & normal,
-                          const typename K::Vector_3 & diff,
-                          const K& k)
-{
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-  RT dot, squared_length;
-  dot = wdot(normal, diff, k);
-  squared_length = wdot(normal, normal, k);
-  return FT(dot*dot) /
-    FT(wmult((K*)0, squared_length, diff.hw(), diff.hw()));
-}
-
 template <class K>
 void
 squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
@@ -213,63 +197,24 @@ squared_distance_to_plane_RT(const typename K::Vector_3 & normal,
                              const K& k)
 {
   typedef typename K::RT RT;
-  typedef typename K::FT FT;
   RT dot, squared_length;
   dot = wdot(normal, diff, k);
   squared_length = wdot(normal, normal, k);
-  num =  square(dot);
+  num = square(dot);
   den = wmult((K*)0, squared_length, diff.hw(), diff.hw());
 }
 
-
 template <class K>
 typename K::FT
-squared_distance_to_line(const typename K::Vector_3 & dir,
-                         const typename K::Vector_3 & diff,
-                         const K& k)
+squared_distance_to_plane(const typename K::Vector_3 & normal,
+                          const typename K::Vector_3 & diff,
+                          const K& k)
 {
-  typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
-  Vector_3 wcr = wcross(dir, diff, k);
-  return FT(wcr*wcr)/FT(wmult(
-                              (K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw()));
-}
-
-template <class K>
-void
-squared_distance_to_line_RT(const typename K::Vector_3 & dir,
-                            const typename K::Vector_3 & diff,
-                            typename K::RT& num,
-                            typename K::RT& den,
-                            const K& k,
-                            Cartesian_tag)
-{
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-  Vector_3 wcr = wcross(dir, diff, k);
-  num = wcr*wcr;
-  den = wmult((K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw());
-}
-
-template <class K>
-void
-squared_distance_to_line_RT(const typename K::Vector_3 & dir,
-                            const typename K::Vector_3 & diff,
-                            typename K::RT& num,
-                            typename K::RT& den,
-                            const K& k,
-                            Homogeneous_tag)
-{
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-  Vector_3 wcr = wcross(dir, diff, k);
-  FT ft = FT(wcr*wcr)/FT(wmult(
-                               (K*)0, RT(wdot(dir, dir, k)), diff.hw(), diff.hw()));
-  num = Rational_traits<FT>().numerator(ft);
-  den = Rational_traits<FT>().denominator(ft);
+  RT num, den;
+  squared_distance_to_plane_RT(normal, diff, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
 }
 
 template <class K>
@@ -280,9 +225,23 @@ squared_distance_to_line_RT(const typename K::Vector_3 & dir,
                             typename K::RT& den,
                             const K& k)
 {
-  typedef typename K::Kernel_tag Tag;
-  Tag tag;
-  squared_distance_to_line_RT(dir,diff,num,den,k,tag);
+  typedef typename K::Vector_3 Vector_3;
+  Vector_3 wcr = wcross(dir, diff, k);
+  num = wdot(wcr, wcr, k);
+  den = wmult((K*)0, wdot(dir, dir, k), diff.hw(), diff.hw());
+}
+
+template <class K>
+typename K::FT
+squared_distance_to_line(const typename K::Vector_3 & dir,
+                         const typename K::Vector_3 & diff,
+                         const K& k)
+{
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT num, den;
+  squared_distance_to_line_RT(dir, diff, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
 }
 
 template <class K>
diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h
index 1a883a35245..617db695eb4 100644
--- a/Distance_3/include/CGAL/squared_distance_3_1.h
+++ b/Distance_3/include/CGAL/squared_distance_3_1.h
@@ -30,20 +30,6 @@ namespace CGAL {
 
 namespace internal {
 
-template <class K>
-typename K::FT
-squared_distance(
-    const typename K::Point_3 &pt,
-    const typename K::Line_3 &line,
-    const K& k)
-{
-  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);
-  return internal::squared_distance_to_line(dir, diff, k);
-}
-
 template <class K>
 void
 squared_distance_RT(
@@ -60,6 +46,20 @@ squared_distance_RT(
   return internal::squared_distance_to_line_RT(dir, diff, num, den, k);
 }
 
+template <class K>
+typename K::FT
+squared_distance(
+    const typename K::Point_3 &pt,
+    const typename K::Line_3 &line,
+    const K& k)
+{
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT num, den;
+  squared_distance_RT(pt, line, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
+}
+
 template <class K>
 inline
 typename K::FT
@@ -72,6 +72,31 @@ squared_distance(
 }
 
 
+template <class K>
+void
+squared_distance_RT(
+    const typename K::Point_3 &pt,
+    const typename K::Ray_3 &ray,
+    typename K::RT& num,
+    typename K::RT& den,
+    const K& k)
+{
+  typename K::Construct_vector_3 construct_vector;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+
+  Vector_3 diff = construct_vector(ray.source(), pt);
+  const Vector_3 &dir = ray.direction().vector();
+  if (!is_acute_angle(dir,diff, k) )
+  {
+    num = wdot(diff, diff, k);
+    den = wmult((K*)0, RT(1), diff.hw(), diff.hw());
+    return;
+  }
+
+  squared_distance_to_line_RT(dir, diff, num, den, k);
+}
+
 template <class K>
 typename K::FT
 squared_distance(
@@ -79,6 +104,14 @@ squared_distance(
     const typename K::Ray_3 &ray,
     const K& k)
 {
+// #define CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
+#ifdef CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT num, den;
+  squared_distance_RT(pt, ray, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
+#else
   typename K::Construct_vector_3 construct_vector;
   typedef typename K::Vector_3 Vector_3;
 
@@ -87,73 +120,9 @@ squared_distance(
     if (!is_acute_angle(dir,diff, k) )
         return (typename K::FT)(diff*diff);
     return squared_distance_to_line(dir, diff, k);
+#endif
 }
 
-template <class K>
-void
-squared_distance_RT(
-    const typename K::Point_3 &pt,
-    const typename K::Ray_3 &ray,
-    typename K::RT& num,
-    typename K::RT& den,
-    const K& k,
-    Homogeneous_tag&)
-{
-  typename K::Construct_vector_3 construct_vector;
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-
-    Vector_3 diff = construct_vector(ray.source(), pt);
-    const Vector_3 &dir = ray.direction().vector();
-    if (!is_acute_angle(dir,diff, k) ){
-      FT ft = diff*diff;
-      num = Rational_traits<FT>().numerator(ft);
-      den = Rational_traits<FT>().denominator(ft);
-      return;
-    }
-    squared_distance_to_line_RT(dir, diff, num, den, k);
-}
-
-template <class K>
-void
-squared_distance_RT(
-    const typename K::Point_3 &pt,
-    const typename K::Ray_3 &ray,
-    typename K::RT& num,
-    typename K::RT& den,
-    const K& k,
-    Cartesian_tag& )
-{
-  typename K::Construct_vector_3 construct_vector;
-  typedef typename K::Vector_3 Vector_3;
-  typedef typename K::RT RT;
-
-    Vector_3 diff = construct_vector(ray.source(), pt);
-    const Vector_3 &dir = ray.direction().vector();
-    if (!is_acute_angle(dir,diff, k) ){
-         num = diff*diff;
-         den =  RT(1);
-         return;
-    }
-    squared_distance_to_line_RT(dir, diff, num, den, k);
-}
-
-template <class K>
-void
-squared_distance_RT(
-    const typename K::Point_3 &pt,
-    const typename K::Ray_3 &ray,
-    typename K::RT& num,
-    typename K::RT& den,
-    const K& k)
-{
-  typedef typename K::Kernel_tag Tag;
-  Tag tag;
-  squared_distance_RT(pt,ray,num,den,k,tag);
-}
-
-
 template <class K>
 inline
 typename K::FT
@@ -165,134 +134,7 @@ squared_distance(
     return squared_distance(pt, ray, k);
 }
 
-
-
-
 template <class K>
-typename K::FT
-squared_distance(
-    const typename K::Point_3 &pt,
-    const typename K::Segment_3 &seg,
-    const K& k,
-    const Homogeneous_tag)
-{
-    typename K::Construct_vector_3 construct_vector;
-    typedef typename K::Vector_3 Vector_3;
-    typedef typename K::RT RT;
-    typedef typename K::FT FT;
-    // 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());
-    RT d = wdot(diff,segvec, k);
-    if (d <= (RT)0)
-        return (FT(diff*diff));
-    RT e = wdot(segvec,segvec, k);
-    if ( (d * segvec.hw()) > (e * diff.hw()))
-        return squared_distance(pt, seg.target(), k);
-
-    Vector_3 wcr = wcross(segvec, diff, k);
-    return FT(wcr*wcr)/FT(e * diff.hw() * diff.hw());
-}
-
-template <class K>
-void
-squared_distance_RT(
-    const typename K::Point_3 &pt,
-    const typename K::Segment_3 &seg,
-    typename K::RT& num,
-    typename K::RT& den,
-    const K& k,
-    const Homogeneous_tag)
-{
-    typename K::Construct_vector_3 construct_vector;
-    typedef typename K::Vector_3 Vector_3;
-    typedef typename K::RT RT;
-    typedef typename K::FT FT;
-    // 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());
-    RT d = wdot(diff,segvec, k);
-    if (d <= (RT)0){
-      FT ft = FT(diff*diff);
-      num = Rational_traits<FT>().numerator(ft);
-      den = Rational_traits<FT>().denominator(ft);
-      return;
-    }
-    RT e = wdot(segvec,segvec, k);
-    if ( (d * segvec.hw()) > (e * diff.hw())){
-      FT sd = squared_distance(pt, seg.target(), k);
-      num = Rational_traits<FT>().numerator(sd);
-      den = Rational_traits<FT>().denominator(sd);
-      return;
-    }
-    Vector_3 wcr = wcross(segvec, diff, k);
-    FT ft((wcr*wcr)/FT(e * diff.hw() * diff.hw()));
-    num = Rational_traits<FT>().numerator(ft);
-    den = Rational_traits<FT>().denominator(ft);
-}
-
-template <class K>
-typename K::FT
-squared_distance(
-    const typename K::Point_3 &pt,
-    const typename K::Segment_3 &seg,
-    const K& k,
-    const Cartesian_tag&)
-{
-    typename K::Construct_vector_3 construct_vector;
-    typedef typename K::Vector_3 Vector_3;
-    typedef typename K::RT RT;
-    typedef typename K::FT FT;
-    // 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());
-    RT d = wdot(diff,segvec, k);
-    if (d <= (RT)0)
-        return (FT(diff*diff));
-    RT e = wdot(segvec,segvec, k);
-    if (d > e)
-        return squared_distance(pt, seg.target(), k);
-
-    Vector_3 wcr = wcross(segvec, diff, k);
-    return FT(wcr*wcr)/e;
-}
-
-template <class K>
-void
-squared_distance_RT(
-    const typename K::Point_3 &pt,
-    const typename K::Segment_3 &seg,
-    typename K::RT& num,
-    typename K::RT& den,
-    const K& k,
-    const Cartesian_tag&)
-{
-    typename K::Construct_vector_3 construct_vector;
-    typedef typename K::Vector_3 Vector_3;
-    typedef typename K::RT RT;
-    typedef typename K::FT FT;
-    // 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());
-    RT d = wdot(diff,segvec, k);
-    if (d <= (RT)0){
-      num = diff*diff;
-      den = RT(1);
-      return;
-    }
-    RT e = wdot(segvec,segvec, k);
-    if (d > e){
-      num = squared_distance(pt, seg.target(), k);
-      den = RT(1);
-      return;
-    }
-    Vector_3 wcr = wcross(segvec, diff, k);
-    num = wcr*wcr;
-    den = e;
-}
-
-template <class K>
-inline
 void
 squared_distance_RT(
     const typename K::Point_3 &pt,
@@ -301,24 +143,74 @@ squared_distance_RT(
     typename K::RT& den,
     const K& k)
 {
-  typedef typename K::Kernel_tag Tag;
-  Tag tag;
-  return squared_distance_RT(pt, seg, num, den, k, tag);
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+
+  typename K::Construct_vector_3 construct_vector;
+
+  // assert that the segment is valid (non zero length).
+  const Vector_3 diff_s = construct_vector(seg.source(), pt);
+  const Vector_3 segvec = construct_vector(seg.source(), seg.target());
+
+  const RT d = wdot(diff_s, segvec, k);
+  if (d <= RT(0))
+  {
+    // this is squared_distance(pt, seg.source())
+    num = wdot(diff_s, diff_s, k);
+    den = wmult((K*)0, RT(1), diff_s.hw(), diff_s.hw());
+    return;
+  }
+
+  const RT e = wdot(segvec, segvec, k);
+  if (wmult((K*)0 ,d, segvec.hw()) > wmult((K*)0, e, diff_s.hw()))
+  {
+    // this is squared_distance(pt, seg.target())
+    const Vector_3 diff_t = construct_vector(seg.target(), pt);
+    num = wdot(diff_t, diff_t, k);
+    den = wmult((K*)0, RT(1), diff_t.hw(), diff_t.hw());
+    return;
+  }
+
+  // This is an expanded call to squared_distance_to_line_RT() to avoid recomputing 'e'
+  const Vector_3 wcr = wcross(segvec, diff_s, k);
+  num = wdot(wcr, wcr, k);
+  den = wmult((K*)0, e, diff_s.hw(), diff_s.hw());
 }
 
 template <class K>
-inline
 typename K::FT
 squared_distance(
     const typename K::Point_3 &pt,
     const typename K::Segment_3 &seg,
     const K& k)
 {
-  typedef typename K::Kernel_tag Tag;
-  Tag tag;
-  return squared_distance(pt, seg, k, tag);
-}
+#ifdef CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT num, den;
+  squared_distance_RT(pt, seg, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
+#else
+  typename K::Construct_vector_3 construct_vector;
+  typedef typename K::Vector_3 Vector_3;
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  // 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());
+  RT d = wdot(diff,segvec, k);
+  if (d <= RT(0))
+    return (FT(diff*diff));
 
+  RT e = wdot(segvec,segvec, k);
+  if (wmult((K*)0 ,d, segvec.hw()) > wmult((K*)0, e, diff.hw()))
+    return squared_distance(pt, seg.target(), k);
+
+  // This is an expanded call to squared_distance_to_line() to avoid recomputing 'e'
+  Vector_3 wcr = wcross(segvec, diff, k);
+  return FT(wcr*wcr) / wmult((K*)0, e, diff.hw(), diff.hw());
+#endif
+}
 
 template <class K>
 inline
@@ -331,8 +223,6 @@ squared_distance(
     return squared_distance(pt, seg, k);
 }
 
-
-
 template <class K>
 typename K::Comparison_result
 compare_distance_pssC3(
@@ -352,7 +242,7 @@ compare_distance_pssC3(
       Vector_3 diff = construct_vector(seg1.source(), pt);
       Vector_3 segvec = construct_vector(seg1.source(), seg1.target());
       RT d = wdot(diff,segvec, k);
-      if (d <= (RT)0){
+      if (d <= RT(0)){
         d1 = (FT(diff*diff));
       }else{
         RT e = wdot(segvec,segvec, k);
@@ -369,7 +259,7 @@ compare_distance_pssC3(
       Vector_3 diff = construct_vector(seg2.source(), pt);
       Vector_3 segvec = construct_vector(seg2.source(), seg2.target());
       RT d = wdot(diff,segvec, k);
-      if (d <= (RT)0){
+      if (d <= RT(0)){
         d2 = (FT(diff*diff));
       }else{
         RT e = wdot(segvec,segvec, k);
@@ -405,7 +295,7 @@ compare_distance_ppsC3(
       Vector_3 diff = construct_vector(seg.source(), pt);
       Vector_3 segvec = construct_vector(seg.source(), seg.target());
       RT d = wdot(diff,segvec, k);
-      if (d <= (RT)0){
+      if (d <= RT(0)){
         d2 = (FT(diff*diff));
       }else{
         RT e = wdot(segvec,segvec, k);
diff --git a/Distance_3/include/CGAL/squared_distance_3_2.h b/Distance_3/include/CGAL/squared_distance_3_2.h
index 6a9acd00153..1150516494f 100644
--- a/Distance_3/include/CGAL/squared_distance_3_2.h
+++ b/Distance_3/include/CGAL/squared_distance_3_2.h
@@ -218,49 +218,6 @@ on_left_of_triangle_edge(const typename K::Point_3 & pt,
   return result;
 }
 
-template <class K>
-inline typename K::FT
-squared_distance_to_triangle(
-    const typename K::Point_3 & pt,
-    const typename K::Point_3 & t0,
-    const typename K::Point_3 & t1,
-    const typename K::Point_3 & t2,
-    bool & inside,
-    const K& k)
-{
-  typename K::Construct_vector_3 vector;
-  typedef typename K::Vector_3 Vector_3;
-  const Vector_3 e1 = vector(t0, t1);
-  const Vector_3 oe3 = vector(t0, t2);
-  const Vector_3 normal = wcross(e1, oe3, k);
-
-  if(normal != NULL_VECTOR
-     && on_left_of_triangle_edge(pt, normal, t0, t1, k)
-     && on_left_of_triangle_edge(pt, normal, t1, t2, k)
-     && on_left_of_triangle_edge(pt, normal, t2, t0, k))
-      {
-        // the projection of pt is inside the triangle
-        inside = true;
-        return squared_distance_to_plane(normal, vector(t0, pt), k);
-      }
-      else {
-        // The case normal==NULL_VECTOR covers the case when the triangle
-        // is colinear, or even more degenerate. In that case, we can
-        // simply take also the distance to the three segments.
-        typename K::FT d1 = squared_distance(pt,
-                                             typename K::Segment_3(t2, t0),
-                                             k);
-        typename K::FT d2 = squared_distance(pt,
-                                             typename K::Segment_3(t1, t2),
-                                             k);
-        typename K::FT d3 = squared_distance(pt,
-                                             typename K::Segment_3(t0, t1),
-                                             k);
-
-        return (std::min)( (std::min)(d1, d2), d3);
-      }
-}
-
 template <class K>
 inline void
 squared_distance_to_triangle_RT(
@@ -279,62 +236,50 @@ squared_distance_to_triangle_RT(
   const Vector_3 oe3 = vector(t0, t2);
   const Vector_3 normal = wcross(e1, oe3, k);
 
-  if(normal != NULL_VECTOR
-     && on_left_of_triangle_edge(pt, normal, t0, t1, k)
-     && on_left_of_triangle_edge(pt, normal, t1, t2, k)
-     && on_left_of_triangle_edge(pt, normal, t2, t0, k))
-      {
-        // the projection of pt is inside the triangle
-        inside = true;
-        squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
-      }
-      else {
-        // The case normal==NULL_VECTOR covers the case when the triangle
-        // is colinear, or even more degenerate. In that case, we can
-        // simply take also the distance to the three segments.
+  if(normal == NULL_VECTOR) {
+    // The case normal==NULL_VECTOR covers the case when the triangle
+    // is colinear, or even more degenerate. In that case, we can
+    // simply take also the distance to the three segments.
+    squared_distance_RT(pt, typename K::Segment_3(t2, t0), num, den, k);
 
-        squared_distance_RT(pt,
-                           typename K::Segment_3(t2, t0),
-                           num,
-                           den,
-                           k);
-        typename K::RT num2, den2;
-        squared_distance_RT(pt,
-                            typename K::Segment_3(t1, t2),
-                            num2,
-                            den2,
-                            k);
-        if(compare_quotients(num2,den2,num,den) == SMALLER){
-          num =  num2;
-          den = den2;
-        }
-        squared_distance_RT(pt,
-                            typename K::Segment_3(t0, t1),
-                            num2,
-                            den2,
-                            k);
-        if(compare_quotients(num2,den2,num,den) == SMALLER){
-          num =  num2;
-          den = den2;
-        }
-      }
-}
+    typename K::RT num2, den2;
+    squared_distance_RT(pt, typename K::Segment_3(t1, t2), num2, den2, k);
+    if(compare_quotients(num2,den2, num,den) == SMALLER) {
+      num = num2;
+      den = den2;
+    }
 
-template <class K>
-inline typename K::FT
-squared_distance(
-    const typename K::Point_3 & pt,
-    const typename K::Triangle_3 & t,
-    const K& k)
-{
-  typename K::Construct_vertex_3 vertex;
-  bool inside = false;
-  return squared_distance_to_triangle(pt,
-                                      vertex(t, 0),
-                                      vertex(t, 1),
-                                      vertex(t, 2),
-                                      inside,
-                                      k);
+    // should not be needed since at most 2 edges cover the full triangle in the degenerate case
+    squared_distance_RT(pt, typename K::Segment_3(t0, t1), num2, den2, k);
+    if(compare_quotients(num2,den2, num,den) == SMALLER){
+      num = num2;
+      den = den2;
+    }
+
+    return;
+  }
+
+  const bool b01 = on_left_of_triangle_edge(pt, normal, t0, t1, k);
+  if(!b01) {
+    squared_distance_RT(pt, typename K::Segment_3(t0, t1), num, den, k);
+    return;
+  }
+
+  const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k);
+  if(!b12) {
+    squared_distance_RT(pt, typename K::Segment_3(t1, t2), num, den, k);
+    return;
+  }
+
+  const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k);
+  if(!b20) {
+    squared_distance_RT(pt, typename K::Segment_3(t2, t0), num, den, k);
+    return;
+  }
+
+  // The projection of pt is inside the triangle
+  inside = true;
+  squared_distance_to_plane_RT(normal, vector(t0, pt), num, den, k);
 }
 
 template <class K>
@@ -358,6 +303,82 @@ squared_distance_RT(
                                   k);
 }
 
+template <class K>
+inline typename K::FT
+squared_distance_to_triangle(
+    const typename K::Point_3 & pt,
+    const typename K::Point_3 & t0,
+    const typename K::Point_3 & t1,
+    const typename K::Point_3 & t2,
+    bool & inside,
+    const K& k)
+{
+  typename K::Construct_vector_3 vector;
+  typedef typename K::Vector_3 Vector_3;
+  const Vector_3 e1 = vector(t0, t1);
+  const Vector_3 oe3 = vector(t0, t2);
+  const Vector_3 normal = wcross(e1, oe3, k);
+
+  if(normal == NULL_VECTOR) {
+    // The case normal==NULL_VECTOR covers the case when the triangle
+    // is colinear, or even more degenerate. In that case, we can
+    // simply take also the distance to the three segments.
+    typename K::FT d1 = squared_distance(pt, typename K::Segment_3(t2, t0), k);
+    typename K::FT d2 = squared_distance(pt, typename K::Segment_3(t1, t2), k);
+
+    // should not be needed since at most 2 edges cover the full triangle in the degenerate case
+    typename K::FT d3 = squared_distance(pt, typename K::Segment_3(t0, t1), k);
+
+    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, typename K::Segment_3(t0, t1), k);
+  }
+
+  const bool b12 = on_left_of_triangle_edge(pt, normal, t1, t2, k);
+  if(!b12) {
+    return internal::squared_distance(pt, typename K::Segment_3(t1, t2), k);
+  }
+
+  const bool b20 = on_left_of_triangle_edge(pt, normal, t2, t0, k);
+  if(!b20) {
+    return internal::squared_distance(pt, typename K::Segment_3(t2, t0), k);
+  }
+
+  // The projection of pt is inside the triangle
+  inside = true;
+  return squared_distance_to_plane(normal, vector(t0, pt), k);
+}
+
+template <class K>
+inline typename K::FT
+squared_distance(
+    const typename K::Point_3 & pt,
+    const typename K::Triangle_3 & t,
+    const K& k)
+{
+  // @todo if the factorized code is used in the end, then the function above, i.e.
+  // squared_distance_to_triangle(), should be removed
+#ifdef CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
+  typedef typename K::RT RT;
+  typedef typename K::FT FT;
+  RT num, den;
+  squared_distance_RT(pt, t, num, den, k);
+  return Rational_traits<FT>().make_rational(num, den);
+#else
+  typename K::Construct_vertex_3 vertex;
+  bool inside = false;
+  return squared_distance_to_triangle(pt,
+                                      vertex(t, 0),
+                                      vertex(t, 1),
+                                      vertex(t, 2),
+                                      inside,
+                                      k);
+#endif
+}
+
 } // namespace internal
 
 

From 430defe1104d3427d1e5b3e599287ae021579506 Mon Sep 17 00:00:00 2001
From: Andreas Fabri <andreas.fabri@geometryfactory.com>
Date: Tue, 11 May 2021 11:11:08 +0200
Subject: [PATCH 10/11] Add dependencies introduced by including
 Rational_traits

---
 Distance_3/package_info/Distance_3/dependencies | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Distance_3/package_info/Distance_3/dependencies b/Distance_3/package_info/Distance_3/dependencies
index 2372c4ab8bb..d8f0e86a70c 100644
--- a/Distance_3/package_info/Distance_3/dependencies
+++ b/Distance_3/package_info/Distance_3/dependencies
@@ -2,7 +2,9 @@ Algebraic_foundations
 Distance_2
 Distance_3
 Installation
+Interval_support
 Kernel_23
+Modular_arithmetic
 Number_types
 Profiling_tools
 STL_Extension

From 2820a9430eee5bd7b54d9eb2f35ee7eecc17e759 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?=
 <mael.rouxel.labbe@geometryfactory.com>
Date: Tue, 25 May 2021 10:46:29 +0200
Subject: [PATCH 11/11] Do not use RT functions in FT functions as it causes a
 slowdown

---
 .../include/CGAL/squared_distance_3_1.h       | 43 ++++++++-----------
 .../include/CGAL/squared_distance_3_2.h       | 12 +-----
 2 files changed, 20 insertions(+), 35 deletions(-)

diff --git a/Distance_3/include/CGAL/squared_distance_3_1.h b/Distance_3/include/CGAL/squared_distance_3_1.h
index 617db695eb4..fdc02382dca 100644
--- a/Distance_3/include/CGAL/squared_distance_3_1.h
+++ b/Distance_3/include/CGAL/squared_distance_3_1.h
@@ -104,23 +104,24 @@ squared_distance(
     const typename K::Ray_3 &ray,
     const K& k)
 {
-// #define CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
-#ifdef CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-  RT num, den;
-  squared_distance_RT(pt, ray, num, den, k);
-  return Rational_traits<FT>().make_rational(num, den);
-#else
-  typename K::Construct_vector_3 construct_vector;
+  // This duplicates code from the _RT functions, but it is a slowdown to do something like:
+  //
+  //   RT num, den;
+  //   squared_distance_RT(pt, ray, num, den, k);
+  //   return Rational_traits<FT>().make_rational(num, den);
+  //
+  // See https://github.com/CGAL/cgal/pull/5680
+
   typedef typename K::Vector_3 Vector_3;
 
-    Vector_3 diff = construct_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 squared_distance_to_line(dir, diff, k);
-#endif
+  typename K::Construct_vector_3 construct_vector = k.construct_vector_3_object();
+
+  Vector_3 diff = construct_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 squared_distance_to_line(dir, diff, k);
 }
 
 template <class K>
@@ -184,17 +185,12 @@ squared_distance(
     const typename K::Segment_3 &seg,
     const K& k)
 {
-#ifdef CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-  RT num, den;
-  squared_distance_RT(pt, seg, num, den, k);
-  return Rational_traits<FT>().make_rational(num, den);
-#else
-  typename K::Construct_vector_3 construct_vector;
   typedef typename K::Vector_3 Vector_3;
   typedef typename K::RT RT;
   typedef typename K::FT FT;
+
+  typename K::Construct_vector_3 construct_vector = k.construct_vector_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());
@@ -209,7 +205,6 @@ squared_distance(
   // This is an expanded call to squared_distance_to_line() to avoid recomputing 'e'
   Vector_3 wcr = wcross(segvec, diff, k);
   return FT(wcr*wcr) / wmult((K*)0, e, diff.hw(), diff.hw());
-#endif
 }
 
 template <class K>
diff --git a/Distance_3/include/CGAL/squared_distance_3_2.h b/Distance_3/include/CGAL/squared_distance_3_2.h
index 1150516494f..c5cd5cf1694 100644
--- a/Distance_3/include/CGAL/squared_distance_3_2.h
+++ b/Distance_3/include/CGAL/squared_distance_3_2.h
@@ -359,16 +359,7 @@ squared_distance(
     const typename K::Triangle_3 & t,
     const K& k)
 {
-  // @todo if the factorized code is used in the end, then the function above, i.e.
-  // squared_distance_to_triangle(), should be removed
-#ifdef CGAL_DISTANCE_3_USE_RT_FUNCTIONS_IN_FT_FUNCTIONS
-  typedef typename K::RT RT;
-  typedef typename K::FT FT;
-  RT num, den;
-  squared_distance_RT(pt, t, num, den, k);
-  return Rational_traits<FT>().make_rational(num, den);
-#else
-  typename K::Construct_vertex_3 vertex;
+  typename K::Construct_vertex_3 vertex = k.construct_vertex_3_object();
   bool inside = false;
   return squared_distance_to_triangle(pt,
                                       vertex(t, 0),
@@ -376,7 +367,6 @@ squared_distance(
                                       vertex(t, 2),
                                       inside,
                                       k);
-#endif
 }
 
 } // namespace internal