diff --git a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
index dfabb50a0b0..e7fd18b8842 100644
--- a/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
+++ b/Principal_component_analysis/include/CGAL/PCA_util_Eigen.h
@@ -17,11 +17,18 @@
 #include <CGAL/Eigen_diagonalize_traits.h>
 #include <CGAL/Default_diagonalize_traits.h>
 #include <CGAL/Dimension.h>
+#include <CGAL/Subiterator.h>
 
 namespace CGAL {
 
 namespace internal {
 
+template <typename FT>
+FT approximate_cbrt (const FT& x)
+{
+  return static_cast<FT>(std::cbrt (CGAL::to_double(x)));
+}
+
 // assemble covariance matrix from a triangle set
 template < typename InputIterator,
            typename K >
@@ -49,9 +56,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 1.0/12.0, 1.0/24.0, 1.0/24.0,
-            1.0/24.0, 1.0/12.0, 1.0/24.0,
-            1.0/24.0, 1.0/24.0, 1.0/12.0;
+  moment << FT(1.0/12.0), FT(1.0/24.0), FT(1.0/24.0),
+            FT(1.0/24.0), FT(1.0/12.0), FT(1.0/24.0),
+            FT(1.0/24.0), FT(1.0/24.0), FT(1.0/12.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -67,7 +74,7 @@ assemble_covariance_matrix_3(InputIterator first,
                       t[0].y(), t[1].y(), t[2].y(),
                       t[0].z(), t[1].z(), t[2].z();
 
-    FT area = std::sqrt(t.squared_area());
+    FT area = CGAL::approximate_sqrt(t.squared_area());
 
     // skip zero measure primitives
     if(area == (FT)0.0)
@@ -89,14 +96,16 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 
@@ -141,18 +150,25 @@ assemble_covariance_matrix_3(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT x0 = t[0].x();
-    FT y0 = t[0].y();
-    FT z0 = t[0].z();
-    FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0,
-                   t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
-                   t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
-    Matrix transformation (delta);
-    FT volume = t.volume();
+    FT x0 = t.xmin();
+    FT y0 = t.ymin();
+    FT z0 = t.zmin();
+
+    FT x1 = t.xmax();
+    FT y1 = t.ymax();
+    FT z1 = t.zmax();
+
+    Matrix transformation;
+    transformation << x1 - x0, 0      , 0      ,
+                      0      , y1 - y0, 0      ,
+                      0      , 0      , z1 - z0;
+
+    FT volume = (x1-x0) * (y1-y0) * (z1-z0);
 
     // skip zero measure primitives
     if(volume == (FT)0.0)
       continue;
+    CGAL_assertion(volume > 0.0);
 
     // Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
 
@@ -160,9 +176,9 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation = volume * transformation * moment * transformation.transpose();
 
     // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
-    FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
-    FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
-    FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
+    FT xav0 = (x1 - x0) / FT(2.0);
+    FT yav0 = (y1 - y0) / FT(2.0);
+    FT zav0 = (z1 - z0) / FT(2.0);
 
     // and add to covariance matrix
     covariance[0] += transformation(0,0) + volume * (2*x0*xav0 + x0*x0);
@@ -175,6 +191,8 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += volume;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
   covariance[0] += mass * (- c.x() * c.x());
@@ -198,6 +216,32 @@ assemble_covariance_matrix_3(InputIterator first,
                              const CGAL::Dimension_tag<2>&,
                              const Eigen_diagonalize_traits<typename K::FT, 3>&)
 {
+#if 1
+  typedef typename K::FT              FT;
+  typedef typename K::Iso_cuboid_3    Iso_cuboid;
+  typedef typename K::Triangle_3      Triangle;
+  auto converter = [](const Iso_cuboid& c, int idx) -> Triangle
+    {
+      // Decomposition of 6 faces of the cuboid into 12 triangles
+      static constexpr std::array<std::array<int, 3>, 12 > indices
+      = {{ { 0, 1, 2 }, { 0, 2, 3 }, { 2, 3, 4 }, { 2, 4, 7 },
+           { 3, 4, 5 }, { 3, 5, 0 }, { 4, 5, 6 }, { 4, 6, 7 },
+           { 5, 6, 1 }, { 5, 1, 0 }, { 6, 7, 2 }, { 6, 2, 1 } }};
+      return Triangle (c[indices[idx][0]], c[indices[idx][1]], c[indices[idx][2]]);
+    };
+
+  assemble_covariance_matrix_3
+    (make_subiterator<Triangle, 12> (first, converter),
+     make_subiterator<Triangle, 12> (beyond),
+     covariance, c, K(), (Triangle*)nullptr, CGAL::Dimension_tag<2>(),
+     Eigen_diagonalize_traits<FT, 3>());
+
+#else
+  // This variant uses the standard formulas but seems to be broken
+  // (line/plane estimated appear to be wrong). In the absence of a
+  // reliable fix so far, the above workaround applying PCA to a
+  // decomposition of the cuboid into triangles is used.
+
   typedef typename K::FT FT;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
   typedef typename Eigen::Matrix<FT, 3, 3> Matrix;
@@ -225,26 +269,25 @@ assemble_covariance_matrix_3(InputIterator first,
     const Iso_cuboid& t = *it;
 
     // defined for convenience.
-    FT x0 = t[0].x();
-    FT y0 = t[0].y();
-    FT z0 = t[0].z();
-    FT delta[9] = {t[1].x()-x0, t[3].x()-x0, t[5].x()-x0,
-                   t[1].y()-y0, t[3].y()-y0, t[5].y()-y0,
-                   t[1].z()-z0, t[3].z()-z0, t[5].z()-z0};
-    Matrix transformation (delta);
-    FT area = std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
-                  delta[6]*delta[6],1/3.0)*std::pow(delta[1]*delta[1] +
-                  delta[4]*delta[4] + delta[7]*delta[7],1/3.0)*2 +
-                  std::pow(delta[0]*delta[0] + delta[3]*delta[3] +
-                  delta[6]*delta[6],1/3.0)*std::pow(delta[2]*delta[2] +
-                  delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2 +
-                  std::pow(delta[1]*delta[1] + delta[4]*delta[4] +
-                  delta[7]*delta[7],1/3.0)*std::pow(delta[2]*delta[2] +
-                  delta[5]*delta[5] + delta[8]*delta[8],1/3.0)*2;
+    FT x0 = t.xmin();
+    FT y0 = t.ymin();
+    FT z0 = t.zmin();
+
+    FT x1 = t.xmax();
+    FT y1 = t.ymax();
+    FT z1 = t.zmax();
+
+    Matrix transformation;
+    transformation << x1 - x0, 0      , 0      ,
+                      0      , y1 - y0, 0      ,
+                      0      , 0      , z1 - z0;
+
+    FT area = FT(2) * ((x1-x0)*(y1-y0) + (x1-x0)*(z1-z0) + (y1-y0)*(z1-z0));
 
     // skip zero measure primitives
     if(area == (FT)0.0)
       continue;
+    CGAL_assertion(area > 0.0);
 
     // Find the 2nd order moment for the cuboid wrt to the origin by an affine transformation.
 
@@ -252,9 +295,9 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation = area * transformation * moment * transformation.transpose();
 
     // Translate the 2nd order moment to the minimum corner (x0,y0,z0) of the cuboid.
-    FT xav0 = (delta[0] + delta[1] + delta[2])/4.0;
-    FT yav0 = (delta[3] + delta[4] + delta[5])/4.0;
-    FT zav0 = (delta[6] + delta[7] + delta[8])/4.0;
+    FT xav0 = (x1 - x0) / (2.0);
+    FT yav0 = (y1 - y0) / (2.0);
+    FT zav0 = (z1 - z0) / (2.0);
 
     // and add to covariance matrix
     covariance[0] += transformation(0,0) + area * (2*x0*xav0 + x0*x0);
@@ -267,15 +310,17 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
-
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
+#endif
 }
 
 // assemble covariance matrix from a sphere set
@@ -305,9 +350,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 4.0/15.0, 0.0,      0.0,
-            0.0,      4.0/15.0, 0.0,
-            0.0,      0.0,      4.0/15.0;
+  moment << FT(4.0/15.0), FT(0.0),      FT(0.0),
+            FT(0.0),      FT(4.0/15.0), FT(0.0),
+            FT(0.0),      FT(0.0),      FT(4.0/15.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -318,7 +363,7 @@ assemble_covariance_matrix_3(InputIterator first,
     const Sphere& t = *it;
 
     // defined for convenience.
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     Matrix transformation;
     transformation << radius, 0.0, 0.0,
                       0.0, radius, 0.0,
@@ -350,14 +395,16 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += volume;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 // assemble covariance matrix from a sphere set
@@ -387,9 +434,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 4.0/3.0, 0.0,     0.0,
-            0.0,     4.0/3.0, 0.0,
-            0.0,     0.0,     4.0/3.0;
+  moment << FT(4.0/3.0), FT(0.0),     FT(0.0),
+            FT(0.0),     FT(4.0/3.0), FT(0.0),
+            FT(0.0),     FT(0.0),     FT(4.0/3.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -401,7 +448,7 @@ assemble_covariance_matrix_3(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     Matrix transformation;
     transformation << radius, 0.0,    0.0,
                       0.0,    radius, 0.0,
@@ -433,14 +480,16 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 
@@ -471,9 +520,9 @@ assemble_covariance_matrix_3(InputIterator first,
   //     5
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 1.0/60.0,  1.0/120.0, 1.0/120.0,
-            1.0/120.0, 1.0/60.0,  1.0/120.0,
-            1.0/120.0, 1.0/120.0, 1.0/60.0;
+  moment << FT(1.0/60.0),  FT(1.0/120.0), FT(1.0/120.0),
+            FT(1.0/120.0), FT(1.0/60.0),  FT(1.0/120.0),
+            FT(1.0/120.0), FT(1.0/120.0), FT(1.0/60.0);
 
   Matrix accum; // zero by default
   accum << 0, 0, 0, 0, 0, 0, 0, 0, 0;
@@ -548,9 +597,9 @@ assemble_covariance_matrix_3(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   Matrix moment;
-  moment << 1.0/3.0, 0.5/3.0, 0.0,
-            0.5/3.0, 1.0/3.0, 0.0,
-            0.0,     0.0,     0.0;
+  moment << FT(1.0/3.0), FT(0.5/3.0), FT(0.0),
+            FT(0.5/3.0), FT(1.0/3.0), FT(0.0),
+            FT(0.0),     FT(0.0),     FT(0.0);
 
   for(InputIterator it = first;
       it != beyond;
@@ -566,8 +615,7 @@ assemble_covariance_matrix_3(InputIterator first,
     transformation << t[0].x(), t[1].x(), 0.0,
                       t[0].y(), t[1].y(), 0.0,
                       t[0].z(), t[1].z(), 1.0;
-    using std::sqrt;
-    FT length = sqrt(t.squared_length());
+    FT length = CGAL::approximate_sqrt(t.squared_length());
 
     // skip zero measure primitives
     if(length == (FT)0.0)
@@ -589,14 +637,16 @@ assemble_covariance_matrix_3(InputIterator first,
     mass += length;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.z() * c.x());
-  covariance[3] += mass * (-1.0 * c.y() * c.y());
-  covariance[4] += mass * (-1.0 * c.z() * c.y());
-  covariance[5] += mass * (-1.0 * c.z() * c.z());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.z() * c.x());
+  covariance[3] += -mass * (c.y() * c.y());
+  covariance[4] += -mass * (c.z() * c.y());
+  covariance[5] += -mass * (c.z() * c.z());
 
 }
 
diff --git a/Principal_component_analysis/include/CGAL/Subiterator.h b/Principal_component_analysis/include/CGAL/Subiterator.h
new file mode 100644
index 00000000000..186776918be
--- /dev/null
+++ b/Principal_component_analysis/include/CGAL/Subiterator.h
@@ -0,0 +1,99 @@
+// Copyright (c) 2020  GeometryFactory (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+// Author(s) : Simon Giraudot
+//
+
+#ifndef CGAL_PCA_SUBITERATOR_H
+#define CGAL_PCA_SUBITERATOR_H
+
+#include <boost/iterator/iterator_facade.hpp>
+
+namespace CGAL
+{
+
+template <typename InputIterator, typename ValueType, int Size>
+class Subiterator
+  : public boost::iterator_facade<Subiterator<InputIterator, ValueType, Size>,
+                                  ValueType,
+                                  std::input_iterator_tag>
+{
+public:
+  using Self = Subiterator<InputIterator, ValueType, Size>;
+  using Facade = boost::iterator_facade<Self, ValueType, std::input_iterator_tag>;
+
+  using Input_type = typename std::iterator_traits<InputIterator>::value_type;
+  using Output_type = ValueType;
+
+  using Converter = std::function<Output_type(const Input_type&, int)>;
+
+private:
+
+  Converter m_converter;
+  int m_next_index;
+  InputIterator m_base;
+  mutable Output_type m_current;
+
+public:
+
+  Subiterator() { }
+
+  Subiterator(InputIterator begin, const Converter& converter)
+    : m_converter(converter), m_next_index(0), m_base(begin)
+  { }
+
+  Subiterator(InputIterator end)
+    : m_next_index(0), m_base(end)
+  { }
+
+private:
+
+  friend class boost::iterator_core_access;
+
+  void increment()
+  {
+    ++ m_next_index;
+    if (m_next_index == Size)
+    {
+      ++ m_base;
+      m_next_index = 0;
+    }
+  }
+
+  bool equal(const Self& other) const
+  {
+    return this->m_base == other.m_base && this->m_next_index == other.m_next_index;
+  }
+
+  Output_type& dereference() const
+  {
+    m_current = m_converter (*m_base, m_next_index);
+    return const_cast<Output_type&>(m_current);
+  }
+};
+
+template <typename ValueType, int Size, typename InputIterator>
+Subiterator<InputIterator, ValueType, Size>
+make_subiterator (InputIterator begin,
+                  const typename Subiterator<InputIterator, ValueType, Size>::Converter& converter)
+{
+  return Subiterator<InputIterator, ValueType, Size>(begin, converter);
+}
+
+template <typename ValueType, int Size, typename InputIterator>
+Subiterator<InputIterator, ValueType, Size>
+make_subiterator (InputIterator end)
+{
+  return Subiterator<InputIterator, ValueType, Size>(end);
+}
+
+
+} // namespace CGAL
+
+#endif // CGAL_PCA_SUBITERATOR_H
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
index dc45b98497e..d6abf3d1482 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_circles_2.h
@@ -83,7 +83,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     FT delta[4] = {radius, 0.0,
                    0.0, radius};
     Matrix transformation = init_matrix<FT>(2,delta);
@@ -107,11 +107,13 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
@@ -133,7 +135,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 
@@ -187,7 +189,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
     // defined for convenience.
     // FT example = CGAL::to_double(t[0].x());
-    FT radius = std::sqrt(t.squared_radius());
+    FT radius = CGAL::approximate_sqrt(t.squared_radius());
     FT delta[4] = {radius, 0.0,
                    0.0, radius};
     Matrix transformation = init_matrix<FT>(2,delta);
@@ -197,7 +199,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // Find the 2nd order moment for the circle wrt to the origin by an affine transformation.
 
     // Transform the standard 2nd order moment using the transformation matrix
-    transformation = 0.5 * length * transformation * moment * LA::transpose(transformation);
+    transformation = FT(0.5) * length * transformation * moment * LA::transpose(transformation);
 
     // Translate the 2nd order moment to the center of the circle.
     FT x0 = t.center().x();
@@ -211,11 +213,13 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += length;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
@@ -237,7 +241,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for circle set with 1D tag
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
index 23d80a5839f..c0071304d6a 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_cuboids_3.h
@@ -20,6 +20,7 @@
 #include <CGAL/PCA_util.h>
 #include <CGAL/linear_least_squares_fitting_points_3.h>
 #include <CGAL/linear_least_squares_fitting_segments_3.h>
+#include <CGAL/Subiterator.h>
 
 #include <list>
 #include <iterator>
@@ -112,29 +113,22 @@ linear_least_squares_fitting_3(InputIterator first,
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[3],t[0]));
-    segments.push_back(Segment(t[4],t[5]));
-    segments.push_back(Segment(t[5],t[6]));
-    segments.push_back(Segment(t[6],t[7]));
-    segments.push_back(Segment(t[7],t[4]));
-    segments.push_back(Segment(t[5],t[0]));
-    segments.push_back(Segment(t[1],t[6]));
-    segments.push_back(Segment(t[2],t[7]));
-    segments.push_back(Segment(t[3],t[4]));
-  }
+  auto converter = [](const Iso_cuboid& c, int idx) -> Segment
+    {
+      if (idx < 7)
+        return Segment (c[idx], c[idx+1]);
+      if (idx < 10)
+        return Segment (c[idx - 6], c[(idx-1)%8]);
+      if (idx == 10)
+        return Segment (c[5], c[0]);
+      CGAL_assertion (idx == 11);
+      return Segment (c[7], c[4]);
+    };
 
-  // compute fitting plane
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 12> (first, converter),
+     make_subiterator<Segment, 12> (beyond),
+     plane,c,(Segment*)nullptr,k,tag,diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
@@ -154,29 +148,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Point_3 Point;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  auto converter = [](const Iso_cuboid& c, int idx) -> Point { return c[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
-    points.push_back(t[4]);
-    points.push_back(t[5]);
-    points.push_back(t[6]);
-    points.push_back(t[7]);
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 8> (first, converter),
+     make_subiterator<Point, 8> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
@@ -259,33 +240,26 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3 Segment;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  auto converter = [](const Iso_cuboid& c, int idx) -> Segment
+    {
+      if (idx < 7)
+        return Segment (c[idx], c[idx+1]);
+      if (idx < 10)
+        return Segment (c[idx - 6], c[(idx-1)%8]);
+      if (idx == 10)
+        return Segment (c[5], c[0]);
+      CGAL_assertion (idx == 11);
+      return Segment (c[7], c[4]);
+    };
+
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[3],t[0]));
-    segments.push_back(Segment(t[4],t[5]));
-    segments.push_back(Segment(t[5],t[6]));
-    segments.push_back(Segment(t[6],t[7]));
-    segments.push_back(Segment(t[7],t[4]));
-    segments.push_back(Segment(t[5],t[0]));
-    segments.push_back(Segment(t[1],t[6]));
-    segments.push_back(Segment(t[2],t[7]));
-    segments.push_back(Segment(t[3],t[4]));
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 12> (first, converter),
+     make_subiterator<Segment, 12> (beyond),
+     line,c,(Segment*)nullptr,k,tag,diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
@@ -305,29 +279,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Point_3 Point;
   typedef typename K::Iso_cuboid_3 Iso_cuboid;
+  auto converter = [](const Iso_cuboid& c, int idx) -> Point { return c[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_cuboid& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
-    points.push_back(t[4]);
-    points.push_back(t[5]);
-    points.push_back(t[6]);
-    points.push_back(t[7]);
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 8> (first, converter),
+     make_subiterator<Point, 8> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_cuboids_3
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
index 7aad5b86554..3bec5e8d4d5 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_rectangles_2.h
@@ -19,6 +19,7 @@
 #include <CGAL/centroid.h>
 #include <CGAL/Linear_algebraCd.h>
 #include <CGAL/PCA_util.h>
+#include <CGAL/Subiterator.h>
 
 #include <iterator>
 #include <list>
@@ -67,8 +68,8 @@ linear_least_squares_fitting_2(InputIterator first,
   typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0. }};
 
   // assemble 2nd order moment about the origin.
-  FT temp[4] = {1/3.0, 0.25,
-                0.25,  1/3.0};
+  FT temp[4] = {FT(1/3.0), FT(0.25),
+                FT(0.25),  FT(1/3.0)};
   Matrix moment = init_matrix<FT>(2,temp);
 
   for(InputIterator it = first;
@@ -100,8 +101,8 @@ linear_least_squares_fitting_2(InputIterator first,
     transformation = area * transformation * moment * LA::transpose(transformation);
 
     // Translate the 2nd order moment to the center of the rectangle.
-    FT xav0 = (x1-x0)/2.0;
-    FT yav0 = (y2-y0)/2.0;
+    FT xav0 = (x1-x0)/FT(2);
+    FT yav0 = (y2-y0)/FT(2);
     // and add to covariance matrix
     covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);
     covariance[1] += transformation[0][1] + area * (x0*yav0 + xav0*y0 + x0*y0);
@@ -110,11 +111,13 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
@@ -136,7 +139,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for rectangle set with 2D tag
@@ -155,25 +158,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Iso_rectangle_2 Iso_rectangle;
   typedef typename K::Segment_2         Segment_2;
+  auto converter = [](const Iso_rectangle& r, int idx) -> Segment_2 { return Segment_2(r[idx], r[(idx+1)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment_2> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_rectangle& t = *it;
-    segments.push_back(Segment_2(t[0],t[1]));
-    segments.push_back(Segment_2(t[1],t[2]));
-    segments.push_back(Segment_2(t[2],t[3]));
-    segments.push_back(Segment_2(t[3],t[0]));
-  }
-
-  return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,
-                                        (Segment_2*)nullptr, K(),tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Segment_2, 4> (first, converter),
+     make_subiterator<Segment_2, 4> (beyond),
+     line,c,(Segment_2*)nullptr,K(),tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for rectangle set with 1D tag
 
@@ -194,25 +188,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Iso_rectangle_2 Iso_rectangle;
   typedef typename K::Point_2         Point_2;
+  auto converter = [](const Iso_rectangle& r, int idx) -> Point_2 { return r[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point_2> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Iso_rectangle& t = *it;
-    points.push_back(Point_2(t[0]));
-    points.push_back(Point_2(t[1]));
-    points.push_back(Point_2(t[2]));
-    points.push_back(Point_2(t[3]));
-  }
-
-  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,
-                                        (Point_2*)nullptr, K(),tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Point_2, 4> (first, converter),
+     make_subiterator<Point_2, 4> (beyond),
+     line,c,(Point_2*)nullptr,K(),tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for rectangle set with 0D tag
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
index f290e741eef..a6d581d97df 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_2.h
@@ -19,6 +19,7 @@
 #include <CGAL/centroid.h>
 #include <CGAL/Linear_algebraCd.h>
 #include <CGAL/PCA_util.h>
+#include <CGAL/Subiterator.h>
 
 #include <iterator>
 #include <list>
@@ -67,7 +68,7 @@ linear_least_squares_fitting_2(InputIterator first,
 
   // assemble 2nd order moment about the origin.
   FT temp[4] = {1.0, 0.5, 0.5, 1.0};
-  Matrix moment = (1.0/3.0) * init_matrix<FT>(2,temp);
+  Matrix moment = FT(1.0/3.0) * init_matrix<FT>(2,temp);
 
   for(InputIterator it = first;
       it != beyond;
@@ -83,7 +84,7 @@ linear_least_squares_fitting_2(InputIterator first,
                    t[0].y(), t[1].y()};
     Matrix transformation = init_matrix<FT>(2,delta);
     using std::sqrt;
-    FT length = sqrt(t.squared_length());
+    FT length = CGAL::approximate_sqrt(t.squared_length());
     CGAL_assertion(length != 0.0);
 
     // Find the 2nd order moment for the segment wrt to the origin by an affine transformation.
@@ -99,11 +100,13 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += length;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * ( c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   // solve for eigenvalues and eigenvectors.
   // eigen values are sorted in ascending order,
@@ -125,7 +128,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for segment set with 1D tag
@@ -144,21 +147,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename K::Point_2  Point;
   typedef typename K::Segment_2 Segment;
+  auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Segment& s = *it;
-    points.push_back(s[0]);
-    points.push_back(s[1]);
-  }
-  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Point, 2> (first, converter),
+     make_subiterator<Point, 2> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for segment set with 1D tag
 
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
index 79a22b08cec..b685d8d1f41 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_segments_3.h
@@ -19,6 +19,7 @@
 #include <CGAL/centroid.h>
 #include <CGAL/PCA_util.h>
 #include <CGAL/linear_least_squares_fitting_points_3.h>
+#include <CGAL/Subiterator.h>
 
 #include <list>
 #include <iterator>
@@ -74,24 +75,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3  Segment;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Segment& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 2> (first, converter),
+     make_subiterator<Point, 2> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_segments_3
 
 // fits a line to a 3D segment set
@@ -141,24 +134,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Segment_3  Segment;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Segment& s, int idx) -> Point { return s[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Segment& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 2> (first, converter),
+     make_subiterator<Point, 2> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_segments_3
 
 } // end namespace internal
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
index 4875556d73e..cfe029d5f0c 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_tetrahedra_3.h
@@ -22,6 +22,7 @@
 #include <CGAL/linear_least_squares_fitting_points_3.h>
 #include <CGAL/linear_least_squares_fitting_segments_3.h>
 #include <CGAL/linear_least_squares_fitting_triangles_3.h>
+#include <CGAL/Subiterator.h>
 
 #include <iterator>
 
@@ -75,26 +76,17 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Triangle_3 Triangle;
+  auto converter = [](const Tetrahedron& t, int idx) -> Triangle
+    { return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Triangle> triangles;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    triangles.push_back(Triangle(t[0],t[1],t[2]));
-    triangles.push_back(Triangle(t[0],t[2],t[3]));
-    triangles.push_back(Triangle(t[0],t[3],t[1]));
-    triangles.push_back(Triangle(t[3],t[1],t[2]));
- }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),plane,c,(Triangle*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Triangle, 4> (first, converter),
+     make_subiterator<Triangle, 4> (beyond),
+     plane,c,(Triangle*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a plane to a 3D tetrahedron set
@@ -113,28 +105,21 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Segment_3 Segment;
+  auto converter = [](const Tetrahedron& t, int idx) -> Segment
+    {
+      if (idx < 4)
+        return Segment (t[idx], t[(idx+1)%4]);
+      return Segment (t[idx-4], t[idx-2]);
+    };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[1],t[3]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[0],t[2]));
-    segments.push_back(Segment(t[0],t[3]));
- }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 6> (first, converter),
+     make_subiterator<Segment, 6> (beyond),
+     plane,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_tetrahedrons_3
 
@@ -154,26 +139,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Point_3 Point;
+  auto converter = [](const Tetrahedron& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
- }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 4> (first, converter),
+     make_subiterator<Point, 4> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a line to a 3D tetrahedron set
@@ -223,26 +198,17 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Triangle_3 Triangle;
+  auto converter = [](const Tetrahedron& t, int idx) -> Triangle
+    { return Triangle(t[idx], t[(idx+1)%4], t[(idx+2)%4]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Triangle> triangles;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    triangles.push_back(Triangle(t[0],t[1],t[2]));
-    triangles.push_back(Triangle(t[0],t[2],t[3]));
-    triangles.push_back(Triangle(t[0],t[3],t[1]));
-    triangles.push_back(Triangle(t[3],t[1],t[2]));
- }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(triangles.begin(),triangles.end(),line,c,(Triangle*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Triangle, 4> (first, converter),
+     make_subiterator<Triangle, 4> (beyond),
+     line,c,(Triangle*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a line to a 3D tetrahedron set
@@ -261,28 +227,21 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Segment_3 Segment;
+  auto converter = [](const Tetrahedron& t, int idx) -> Segment
+    {
+      if (idx < 4)
+        return Segment (t[idx], t[(idx+1)%4]);
+      return Segment (t[idx-4], t[idx-2]);
+    };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[1],t[3]));
-    segments.push_back(Segment(t[2],t[3]));
-    segments.push_back(Segment(t[0],t[2]));
-    segments.push_back(Segment(t[0],t[3]));
- }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 6> (first, converter),
+     make_subiterator<Segment, 6> (beyond),
+     line,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedrons_3
 
 // fits a line to a 3D tetrahedron set
@@ -301,26 +260,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Tetrahedron_3    Tetrahedron;
   typedef typename K::Point_3 Point;
+  auto converter = [](const Tetrahedron& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Tetrahedron& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-    points.push_back(t[3]);
- }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 4> (first, converter),
+     make_subiterator<Point, 4> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_tetrahedra_3
 
 } // end namespace internal
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
index 2cff2b8081e..ab96ad6d408 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_2.h
@@ -68,8 +68,8 @@ linear_least_squares_fitting_2(InputIterator first,
   typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0. }};
 
   // assemble the 2nd order moment about the origin.
-  FT temp[4] = {1/12.0, 1/24.0,
-                1/24.0, 1/12.0};
+  FT temp[4] = {FT(1/12.0), FT(1/24.0),
+                FT(1/24.0), FT(1/12.0)};
 
   Matrix moment = init_matrix<FT>(2,temp);
 
@@ -88,7 +88,7 @@ linear_least_squares_fitting_2(InputIterator first,
                    t[1].y() - y0, t[2].y() - y0};
 
     Matrix transformation = init_matrix<FT>(2,delta);
-    FT area = 0.5 * std::abs(LA::determinant(transformation));
+    FT area = FT(0.5) * CGAL::abs(LA::determinant(transformation));
     CGAL_assertion(area!=0);
 
     // Find the 2nd order moment for the triangle wrt to the origin by an affine transformation.
@@ -97,8 +97,8 @@ linear_least_squares_fitting_2(InputIterator first,
     transformation = 2 * area * transformation * moment * LA::transpose(transformation);
 
     // Translate the 2nd order moment to (x0,y0).
-    FT xav0 = (delta[0]+delta[1])/3.0;
-    FT yav0 = (delta[2]+delta[3])/3.0;
+    FT xav0 = (delta[0]+delta[1])/FT(3);
+    FT yav0 = (delta[2]+delta[3])/FT(3);
 
     // and add to the covariance matrix
     covariance[0] += transformation[0][0] + area * (x0*xav0*2 + x0*x0);
@@ -108,11 +108,13 @@ linear_least_squares_fitting_2(InputIterator first,
     mass += area;
   }
 
+  CGAL_assertion_msg (mass != FT(0), "Can't compute PCA of null measure.");
+
   // Translate the 2nd order moment calculated about the origin to
   // the center of mass to get the covariance.
-  covariance[0] += mass * (-1.0 * c.x() * c.x());
-  covariance[1] += mass * (-1.0 * c.x() * c.y());
-  covariance[2] += mass * (-1.0 * c.y() * c.y());
+  covariance[0] += -mass * (c.x() * c.x());
+  covariance[1] += -mass * (c.x() * c.y());
+  covariance[2] += -mass * (c.y() * c.y());
 
   //  std::cout<<"cov: "<<covariance[0]*covariance[2]<<" =? "<<covariance[1]*covariance[1]<<std::endl;
 
@@ -136,7 +138,7 @@ linear_least_squares_fitting_2(InputIterator first,
     // isotropic case (infinite number of directions)
     // by default: assemble a line that goes through
     // the centroid and with a default horizontal vector.
-    line = Line(c, Vector(1.0, 0.0));
+    line = Line(c, Vector(FT(1), FT(0)));
     return (FT)0.0;
   }
 } // end linear_least_squares_fitting_2 for triangle set with 2D tag
@@ -157,23 +159,16 @@ linear_least_squares_fitting_2(InputIterator first,
   // types
   typedef typename Kernel::Triangle_2 Triangle;
   typedef typename Kernel::Segment_2  Segment;
+  auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[0]));
-  }
-
-  return linear_least_squares_fitting_2(segments.begin(),segments.end(),line,c,tag,Kernel(),
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_2
+    (make_subiterator<Segment, 3> (first, converter),
+     make_subiterator<Segment, 3> (beyond),
+     line,c,(Segment*)nullptr,Kernel(),tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_2 for triangle set with 1D tag
 
@@ -194,24 +189,16 @@ linear_least_squares_fitting_2(InputIterator first,
 
   typedef typename Kernel::Triangle_2 Triangle;
   typedef typename Kernel::Point_2 Point;
+  auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    points.push_back(Point(t[0]));
-    points.push_back(Point(t[1]));
-    points.push_back(Point(t[2]));
-  }
-
-  return linear_least_squares_fitting_2(points.begin(),points.end(),line,c,tag,Kernel(),
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_2
+    (make_subiterator<Point, 3> (first, converter),
+     make_subiterator<Point, 3> (beyond),
+     line,c,(Point*)nullptr,Kernel(),tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_2 for triangle set with 0D tag
 
 } // end namespace internal
diff --git a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
index 2337428e1a8..a2a1e96bec4 100644
--- a/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
+++ b/Principal_component_analysis/include/CGAL/linear_least_squares_fitting_triangles_3.h
@@ -18,6 +18,7 @@
 #include <CGAL/basic.h>
 #include <CGAL/centroid.h>
 #include <CGAL/PCA_util.h>
+#include <CGAL/Subiterator.h>
 
 #include <list>
 #include <iterator>
@@ -73,24 +74,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Segment_3  Segment;
+  auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[0]));
-  }
-
-  // compute fitting plane
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),plane,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 3> (first, converter),
+     make_subiterator<Segment, 3> (beyond),
+     plane,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_triangles_3
 
@@ -110,23 +103,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-  }
 
-  // compute fitting plane
-  return linear_least_squares_fitting_3(points.begin(),points.end(),plane,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 3> (first, converter),
+     make_subiterator<Point, 3> (beyond),
+     plane,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 
 } // end linear_least_squares_fitting_triangles_3
 
@@ -177,25 +163,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Segment_3  Segment;
+  auto converter = [](const Triangle& t, int idx) -> Segment { return Segment(t[idx], t[(idx+1)%3]); };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
 
-  std::list<Segment> segments;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    segments.push_back(Segment(t[0],t[1]));
-    segments.push_back(Segment(t[1],t[2]));
-    segments.push_back(Segment(t[2],t[0]));
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(segments.begin(),segments.end(),line,c,(Segment*)nullptr,k,tag,
-                                        diagonalize_traits);
-
+  return linear_least_squares_fitting_3
+    (make_subiterator<Segment, 3> (first, converter),
+     make_subiterator<Segment, 3> (beyond),
+     line,c,(Segment*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_triangles_3
 
 // fits a line to a 3D triangle set
@@ -214,24 +191,16 @@ linear_least_squares_fitting_3(InputIterator first,
 {
   typedef typename K::Triangle_3  Triangle;
   typedef typename K::Point_3  Point;
+  auto converter = [](const Triangle& t, int idx) -> Point { return t[idx]; };
 
   // precondition: at least one element in the container.
   CGAL_precondition(first != beyond);
-  std::list<Point> points;
-  for(InputIterator it = first;
-      it != beyond;
-      it++)
-  {
-    const Triangle& t = *it;
-    points.push_back(t[0]);
-    points.push_back(t[1]);
-    points.push_back(t[2]);
-  }
-
-  // compute fitting line
-  return linear_least_squares_fitting_3(points.begin(),points.end(),line,c,(Point*)nullptr,k,tag,
-                                        diagonalize_traits);
 
+  return linear_least_squares_fitting_3
+    (make_subiterator<Point, 3> (first, converter),
+     make_subiterator<Point, 3> (beyond),
+     line,c,(Point*)nullptr,k,tag,
+     diagonalize_traits);
 } // end linear_least_squares_fitting_triangles_3
 
 } // end namespace internal
diff --git a/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp b/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
new file mode 100644
index 00000000000..606f3abbc56
--- /dev/null
+++ b/Principal_component_analysis/test/Principal_component_analysis/test_kernels.cpp
@@ -0,0 +1,107 @@
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/linear_least_squares_fitting_2.h>
+#include <CGAL/linear_least_squares_fitting_3.h>
+
+template <typename K>
+typename K::Point_2 point_2 (int x, int y)
+{ return typename K::Point_2 (typename K::FT(x), typename K::FT(y)); }
+template <typename K>
+typename K::Point_3 point_3 (int x, int y, int z)
+{ return typename K::Point_3 (typename K::FT(x), typename K::FT(y), typename K::FT(z)); }
+
+// Generate default objects so that the test compiles *AND* runs fine
+template <typename K, typename O> O default_object(const O&) { return O(); }
+template <typename K> typename K::Point_2 default_object(const typename K::Point_2&)
+{ return point_2<K>(0,0); }
+template <typename K> typename K::Segment_2 default_object(const typename K::Segment_2&)
+{ return typename K::Segment_2(point_2<K>(0,0), point_2<K>(0,1)); }
+template <typename K> typename K::Circle_2 default_object(const typename K::Circle_2&)
+{ return typename K::Circle_2(point_2<K>(0,0), typename K::FT(1.0)); }
+template <typename K> typename K::Triangle_2 default_object(const typename K::Triangle_2&)
+{ return typename K::Triangle_2(point_2<K>(0,0), point_2<K>(0,1), point_2<K>(1,0)); }
+template <typename K> typename K::Iso_rectangle_2 default_object(const typename K::Iso_rectangle_2&)
+{ return typename K::Iso_rectangle_2(point_2<K>(0,0), point_2<K>(1,1)); }
+template <typename K> typename K::Point_3 default_object(const typename K::Point_3&)
+{ return point_3<K>(0,0,0); }
+template <typename K> typename K::Segment_3 default_object(const typename K::Segment_3&)
+{ return typename K::Segment_3(point_3<K>(0,0,0), point_3<K>(0,0,1)); }
+template <typename K> typename K::Sphere_3 default_object(const typename K::Sphere_3&)
+{ return typename K::Sphere_3(point_3<K>(0,0,0), typename K::FT(1.0)); }
+template <typename K> typename K::Triangle_3 default_object(const typename K::Triangle_3&)
+{ return typename K::Triangle_3(point_3<K>(0,0,0), point_3<K>(0,0,1), point_3<K>(0,1,0)); }
+template <typename K> typename K::Tetrahedron_3 default_object(const typename K::Tetrahedron_3&)
+{ return typename K::Tetrahedron_3(point_3<K>(0,0,0), point_3<K>(0,0,1),
+                                   point_3<K>(0,1,0), point_3<K>(1,0,0)); }
+template <typename K> typename K::Iso_cuboid_3 default_object(const typename K::Iso_cuboid_3&)
+{ return typename K::Iso_cuboid_3(point_3<K>(0,0,0), point_3<K>(1,1,1)); }
+
+template <typename Kernel, typename Object, int dim>
+void test_2d()
+{
+  std::array<Object, 1> dummy = { default_object<Kernel>(Object()) };
+  typename Kernel::Line_2 line;
+  typename Kernel::Point_2 centroid;
+  CGAL::linear_least_squares_fitting_2 (dummy.begin(), dummy.end(), line, centroid,
+                                        CGAL::Dimension_tag<dim>(), Kernel(),
+                                        CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 2>());
+}
+
+template <typename Kernel, typename Object, int dim>
+void test_3d()
+{
+  std::array<Object, 1> dummy = { default_object<Kernel>(Object()) };
+  typename Kernel::Line_3 line;
+  typename Kernel::Plane_3 plane;
+  typename Kernel::Point_3 centroid;
+  CGAL::linear_least_squares_fitting_3 (dummy.begin(), dummy.end(), line, centroid,
+                                        CGAL::Dimension_tag<dim>(), Kernel(),
+                                        CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>());
+  CGAL::linear_least_squares_fitting_3 (dummy.begin(), dummy.end(), plane, centroid,
+                                        CGAL::Dimension_tag<dim>(), Kernel(),
+                                        CGAL::Eigen_diagonalize_traits<typename Kernel::FT, 3>());
+}
+
+template <typename Kernel>
+void test_kernel()
+{
+  test_2d<Kernel, typename Kernel::Point_2, 0>();
+  test_2d<Kernel, typename Kernel::Segment_2, 1>();
+  test_2d<Kernel, typename Kernel::Segment_2, 0>();
+  test_2d<Kernel, typename Kernel::Circle_2, 2>();
+  test_2d<Kernel, typename Kernel::Circle_2, 1>();
+  test_2d<Kernel, typename Kernel::Triangle_2, 2>();
+  test_2d<Kernel, typename Kernel::Triangle_2, 1>();
+  test_2d<Kernel, typename Kernel::Triangle_2, 0>();
+  test_2d<Kernel, typename Kernel::Iso_rectangle_2, 2>();
+  test_2d<Kernel, typename Kernel::Iso_rectangle_2, 1>();
+  test_2d<Kernel, typename Kernel::Iso_rectangle_2, 0>();
+  test_3d<Kernel, typename Kernel::Point_3, 0>();
+  test_3d<Kernel, typename Kernel::Segment_3, 1>();
+  test_3d<Kernel, typename Kernel::Segment_3, 0>();
+  test_3d<Kernel, typename Kernel::Sphere_3, 3>();
+  test_3d<Kernel, typename Kernel::Sphere_3, 2>();
+  test_3d<Kernel, typename Kernel::Triangle_3, 2>();
+  test_3d<Kernel, typename Kernel::Triangle_3, 1>();
+  test_3d<Kernel, typename Kernel::Triangle_3, 0>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 3>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 2>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 1>();
+  test_3d<Kernel, typename Kernel::Tetrahedron_3, 0>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 3>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 2>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 1>();
+  test_3d<Kernel, typename Kernel::Iso_cuboid_3, 0>();
+}
+
+int main()
+{
+  test_kernel<CGAL::Simple_cartesian<float> >();
+  test_kernel<CGAL::Simple_cartesian<double> >();
+  test_kernel<CGAL::Exact_predicates_inexact_constructions_kernel>();
+  test_kernel<CGAL::Exact_predicates_exact_constructions_kernel>();
+
+  return EXIT_SUCCESS;
+}
diff --git a/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp b/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp
index 26bddf76fe7..c9a4e42ee61 100644
--- a/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp
+++ b/Principal_component_analysis/test/Principal_component_analysis/test_validity.cpp
@@ -225,7 +225,7 @@ void assert_quality (const std::vector<Point_3>& points, const Fitted& fitted)
   std::cerr << "mean distance = " << mean_dist << std::endl;
 
   CGAL_assertion_code
-    (double limit = 1e-3 * std::sqrt (CGAL::squared_distance (points.front(), points.back())));
+    (double limit = 1e-5 * std::sqrt (CGAL::squared_distance (points.front(), points.back())));
   CGAL_assertion (mean_dist < limit);
 }
 
diff --git a/Principal_component_analysis_LGPL/include/CGAL/centroid.h b/Principal_component_analysis_LGPL/include/CGAL/centroid.h
index 63fae0354ce..5d6bee250d2 100644
--- a/Principal_component_analysis_LGPL/include/CGAL/centroid.h
+++ b/Principal_component_analysis_LGPL/include/CGAL/centroid.h
@@ -120,9 +120,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Segment& s = *it;
-    using std::abs;
-    using std::sqrt;
-    FT length = sqrt(abs(s.squared_length()));
+    FT length = CGAL::approximate_sqrt(CGAL::abs(s.squared_length()));
     Point c = K().construct_midpoint_2_object()(s[0],s[1]);
     v = v + length * (c - ORIGIN);
     sum_lengths += length;
@@ -216,7 +214,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Triangle& triangle = *it;
-    FT unsigned_area = std::abs(triangle.area());
+    FT unsigned_area = CGAL::abs(triangle.area());
     Point c = K().construct_centroid_2_object()(triangle[0],triangle[1],triangle[2]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
@@ -252,7 +250,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Circle& s = *it;
-    FT radius = std::sqrt(s.squared_radius());
+    FT radius = CGAL::approximate_sqrt(s.squared_radius());
     Point c = s.center();
     v = v + radius * (c - ORIGIN);
     sum_lengths += radius;
@@ -381,7 +379,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Iso_rectangle& r = *it;
-    FT unsigned_area = std::abs(r.area());
+    FT unsigned_area = CGAL::abs(r.area());
     Point c = K().construct_centroid_2_object()(r[0],r[1],r[2],r[3]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
@@ -442,8 +440,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Segment& s = *it;
-    using std::sqrt;
-    FT length = sqrt(s.squared_length());
+    FT length = CGAL::approximate_sqrt(s.squared_length());
     Point c = CGAL::midpoint(s.source(),s.target());
     // Point c = K().construct_midpoint_3_object()(s[0],s[1]);
     //Point c = Point((s[0][0] + s[1][0])/2.0, (s[0][1] + s[1][1])/2.0, (s[0][2] + s[1][2])/2.0);
@@ -539,7 +536,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Triangle& triangle = *it;
-    FT unsigned_area = std::sqrt(triangle.squared_area());
+    FT unsigned_area = CGAL::approximate_sqrt(triangle.squared_area());
     Point c = K().construct_centroid_3_object()(triangle[0],triangle[1],triangle[2]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
@@ -608,7 +605,7 @@ centroid(InputIterator begin,
       it++)
   {
     const Sphere& sphere = *it;
-    FT unsigned_volume = sphere.squared_radius() * std::sqrt(sphere.squared_radius());
+    FT unsigned_volume = sphere.squared_radius() * CGAL::approximate_sqrt(sphere.squared_radius());
     Point c = sphere.center();
     v = v + unsigned_volume * (c - ORIGIN);
     sum_volumes += unsigned_volume;
@@ -717,7 +714,7 @@ centroid(InputIterator begin,
   {
     const Iso_cuboid& cuboid = *it;
     FT unsigned_area = 2 * ((cuboid.xmax()-cuboid.xmin())*(cuboid.ymax()-cuboid.ymin()) + (cuboid.xmax()-cuboid.xmin())*(cuboid.zmax()-cuboid.zmin()) + (cuboid.ymax()-cuboid.ymin())*(cuboid.zmax()-cuboid.zmin()));
-    Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
+    Point c = K().construct_midpoint_3_object()(cuboid[0],cuboid[7]);
     v = v + unsigned_area * (c - ORIGIN);
     sum_areas += unsigned_area;
   }
@@ -750,7 +747,7 @@ centroid(InputIterator begin,
   {
     const Iso_cuboid& cuboid = *it;
     FT unsigned_volume = cuboid.volume();
-    Point c = K().construct_centroid_3_object()(cuboid[0],cuboid[1],cuboid[3],cuboid[5]);
+    Point c = K().construct_midpoint_3_object()(cuboid[0],cuboid[7]);
     v = v + unsigned_volume * (c - ORIGIN);
     sum_volumes += unsigned_volume;
   }
diff --git a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
index 60aee088ce1..38ee3c7c4c7 100644
--- a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
+++ b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
@@ -26,10 +26,20 @@
 
 //#define DO_NOT_USE_EIGEN_COMPUTEDIRECT_FOR_DIAGONALIZATION
 
-#include <CGAL/array.h>
+#include <CGAL/number_utils.h>
+#include <CGAL/double.h>
 
 namespace CGAL {
 
+namespace internal {
+
+template <typename FT>
+struct Restricted_FT { typedef double type; };
+template <>
+struct Restricted_FT<float> { typedef float type; };
+
+}
+
 /// \ingroup PkgSolverInterfaceRef
 ///
 /// The class `Eigen_diagonalize_traits` provides an interface to the
@@ -47,14 +57,15 @@ namespace CGAL {
 template <typename FT, unsigned int dim = 3>
 class Eigen_diagonalize_traits
 {
+  typedef typename internal::Restricted_FT<FT>::type RFT;
 public:
   typedef std::array<FT, dim>                  Vector;
   typedef std::array<FT, dim*dim>              Matrix;
   typedef std::array<FT, (dim * (dim+1) / 2)>  Covariance_matrix;
 
 private:
-  typedef Eigen::Matrix<FT, dim, dim>            EigenMatrix;
-  typedef Eigen::Matrix<FT, dim, 1>              EigenVector;
+  typedef Eigen::Matrix<RFT, dim, dim>            EigenMatrix;
+  typedef Eigen::Matrix<RFT, dim, 1>              EigenVector;
 
   /// Construct the covariance matrix
   static EigenMatrix construct_covariance_matrix(const Covariance_matrix& cov)
@@ -65,7 +76,7 @@ private:
     {
       for(std::size_t j=i; j<dim; ++j)
       {
-        m(i,j) = static_cast<float>(cov[(dim * i) + j - ((i * (i+1)) / 2)]);
+        m(i,j) = static_cast<RFT>(CGAL::to_double(cov[(dim * i) + j - ((i * (i+1)) / 2)]));
 
         if(i != j)
           m(j,i) = m(i,j);