Merge pull request #5532 from lrineau/Kernel_23-Do_interesect__without_division-GF

Remove operator/ from CGAL::Mpzf and implement Do_intersect_3 with RT only (without division)

Algebraic_foundations/include/CGAL/number_utils.h

@@ -20,6 +20,7 @@
 #include <CGAL/number_type_config.h>
 #include <CGAL/Algebraic_structure_traits.h>
 #include <CGAL/Real_embeddable_traits.h>
+#include <CGAL/Kernel/Same_uncertainty.h>
 
 namespace CGAL {
 CGAL_NTS_BEGIN_NAMESPACE
@@ -324,6 +325,36 @@ NT approximate_sqrt(const NT& nt)
   return approximate_sqrt(nt, Sqrt());
 }
 
+template <class NT>
+typename Same_uncertainty_nt<Comparison_result, NT>::type
+compare_quotients(const NT& xnum, const NT& xden,
+                  const NT& ynum, const NT& yden)
+{
+  // No assumptions on the sign of  den  are made
+
+  // code assumes that SMALLER == - 1;
+  CGAL_precondition( SMALLER == static_cast<Comparison_result>(-1) );
+
+  int xsign =  sign(xnum) *  sign(xden) ;
+  int ysign =  sign(ynum) *  sign(yden) ;
+  if (xsign == 0) return static_cast<Comparison_result>(-ysign);
+  if (ysign == 0) return static_cast<Comparison_result>(xsign);
+  // now (x != 0) && (y != 0)
+  int diff = xsign - ysign;
+  if (diff == 0)
+    {
+      int msign =  sign(xden) *  sign(yden);
+      NT leftop  = NT(xnum * yden * msign);
+      NT rightop = NT(ynum * xden * msign);
+      return  CGAL::compare(leftop, rightop);
+    }
+  else
+    {
+      return (xsign < ysign) ? SMALLER : LARGER;
+    }
+}
+
+
 CGAL_NTS_END_NAMESPACE
 } //namespace CGAL
 

Bounding_volumes/include/CGAL/Min_sphere_of_spheres_d/Min_sphere_of_spheres_d_impl.h

@@ -30,6 +30,8 @@
 
 namespace CGAL_MINIBALL_NAMESPACE {
 
+namespace  Bounding_volumes {
+
   template<typename FT>
   inline bool compare(const FT& a,const FT& b,
                       const FT& ap,const FT& bp) {
@@ -56,6 +58,8 @@ namespace CGAL_MINIBALL_NAMESPACE {
     }
   }
 
+} // namespace Bounding_volumes
+
   template<class Traits>
   void Min_sphere_of_spheres_d<Traits>::update(LP_algorithm) {
     using namespace Min_sphere_of_spheres_d_impl;
@@ -168,7 +172,7 @@ namespace CGAL_MINIBALL_NAMESPACE {
                                          t.center_cartesian_begin(*l[k]),FT(0),std::plus<FT>(),
         Subtract_and_square<FT>());
 
-      if (compare(max,maxp,t.radius(*l[k]),dist)) {
+      if (Bounding_volumes::compare(max,maxp,t.radius(*l[k]),dist)) {
         max   = t.radius(*l[k]);
         maxp  = dist;
         i     = k;
@@ -203,7 +207,7 @@ namespace CGAL_MINIBALL_NAMESPACE {
          Subtract_and_square_to_double<FT>());
 
       const double r = CGAL_MINIBALL_NTS to_double(t.radius(*l[k]));
-      if (compare(max,maxp,r,dist)) {
+      if (Bounding_volumes::compare(max,maxp,r,dist)) {
         max   = r;
         maxp  = dist;
         i     = k;

Cartesian_kernel/include/CGAL/Cartesian/function_objects.h

@@ -253,18 +253,18 @@ namespace CartesianKernelFunctors {
     result_type
     operator()( const Tetrahedron_3& t, const Point_3& p) const
     {
-      FT alpha, beta, gamma;
+      FT alpha, beta, gamma, denom;
 
       Cartesian_internal::solve(t.vertex(1)-t.vertex(0),
                                 t.vertex(2)-t.vertex(0),
                                 t.vertex(3)-t.vertex(0),
-                                p - t.vertex(0), alpha, beta, gamma);
+                                p - t.vertex(0), alpha, beta, gamma, denom);
       if (   (alpha < 0) || (beta < 0) || (gamma < 0)
-          || (alpha + beta + gamma > 1) )
+          || (alpha + beta + gamma > denom) )
           return ON_UNBOUNDED_SIDE;
 
       if (   (alpha == 0) || (beta == 0) || (gamma == 0)
-          || (alpha+beta+gamma == 1) )
+          || (alpha+beta+gamma == denom) )
         return ON_BOUNDARY;
 
       return ON_BOUNDED_SIDE;
@@ -3860,10 +3860,10 @@ namespace CartesianKernelFunctors {
                v1 = t.vertex(1)-o,
                v2 = t.vertex(2)-o;
 
-      FT alpha, beta, gamma;
+      FT alpha, beta, gamma, denum;
-      Cartesian_internal::solve(v0, v1, v2, p-o, alpha, beta, gamma);
+      Cartesian_internal::solve(v0, v1, v2, p-o, alpha, beta, gamma, denum);
       return (alpha >= FT(0)) && (beta >= FT(0)) && (gamma >= FT(0))
-          && ((alpha+beta+gamma == FT(1)));
+          && ((alpha+beta+gamma == denum));
     }
 
     result_type

Cartesian_kernel/include/CGAL/Cartesian/solve_3.h

@@ -24,6 +24,20 @@ namespace CGAL {
 
 namespace Cartesian_internal {
 
+template <class R>
+void solve (const VectorC3<R> &v0,
+            const VectorC3<R> &v1,
+            const VectorC3<R> &v2,
+            const VectorC3<R> &d,
+            typename R::FT &alpha, typename R::FT &beta, typename R::FT &gamma, typename R::FT &denom)
+{
+  CGAL::solve(v0.x(), v0.y(), v0.z(),
+              v1.x(), v1.y(), v1.z(),
+              v2.x(), v2.y(), v2.z(),
+              d.x(),  d.y(),  d.z(),
+              alpha, beta, gamma, denom);
+}
+
 template <class R>
 void solve (const VectorC3<R> &v0,
             const VectorC3<R> &v1,
@@ -32,10 +46,10 @@ void solve (const VectorC3<R> &v0,
             typename R::FT &alpha, typename R::FT &beta, typename R::FT &gamma)
 {
   CGAL::solve(v0.x(), v0.y(), v0.z(),
-        v1.x(), v1.y(), v1.z(),
+              v1.x(), v1.y(), v1.z(),
-        v2.x(), v2.y(), v2.z(),
+              v2.x(), v2.y(), v2.z(),
-        d.x(),  d.y(),  d.z(),
+              d.x(),  d.y(),  d.z(),
-        alpha, beta, gamma);
+              alpha, beta, gamma);
 }
 
 } // namespace Cartesian_internal

Filtered_kernel/include/CGAL/internal/Static_filters/Do_intersect_3.h

@@ -78,6 +78,13 @@ public:
     return Intersections::internal::do_intersect(t,s, SFK());
   }
 
+  result_type
+  operator()(const Triangle_3 &t0, const Triangle_3& t1) const
+  {
+    return Intersections::internal::do_intersect(t0,t1, SFK());
+  }
+
+
   result_type
   operator()(const Bbox_3& b, const Segment_3 &s) const
   {

Intersections_3/include/CGAL/Intersections_3/internal/intersection_3_1_impl.h

@@ -26,7 +26,7 @@
 #include <CGAL/Intersections_3/Iso_cuboid_3_Line_3.h>
 #include <CGAL/utils_classes.h>
 #include <CGAL/squared_distance_3.h>
-
+#include <CGAL/rank.h>
 #include <CGAL/Intersections_3/internal/bbox_intersection_3.h>
 namespace CGAL {
 
@@ -239,7 +239,7 @@ intersection_point(const typename K::Plane_3 &plane1,
 
   const FT den = minor_0*m22 - minor_1*m12 + minor_2*m02; // determinant of M
 
-  if(den == FT(0)){
+  if(is_zero(den)){
     return boost::none;
   }
 
@@ -379,7 +379,7 @@ intersection(const typename K::Line_3 &l1,
   const Vector_3 v3v2 = cross_product(v3,v2);
   const Vector_3 v1v2 = cross_product(v1,v2);
   const FT sl = v1v2.squared_length();
-  if(certainly(sl == FT(0)))
+  if(certainly(is_zero(sl)))
     return intersection_return<typename K::Intersect_3, typename K::Line_3, typename K::Line_3>();
   const FT t = ((v3v2.x()*v1v2.x()) + (v3v2.y()*v1v2.y()) + (v3v2.z()*v1v2.z())) / sl;
 
@@ -523,6 +523,7 @@ do_intersect(const typename K::Segment_3  &s1,
   return false;
 }
 
+
 template <class K>
 typename Intersection_traits<K, typename K::Line_3, typename K::Segment_3>::result_type
 intersection(const typename K::Line_3 &l,
@@ -853,9 +854,9 @@ do_intersect(const typename K::Plane_3 &p,
   typedef typename K::FT FT;
   const FT d2 = CGAL::square(p.a()*s.center().x() +
                              p.b()*s.center().y() +
-                             p.c()*s.center().z() + p.d()) /
+                             p.c()*s.center().z() + p.d());
-      (square(p.a()) + square(p.b()) + square(p.c()));
+
-  return d2 <= s.squared_radius();
+  return d2 <= s.squared_radius() * (square(p.a()) + square(p.b()) + square(p.c()));
 }
 
 template <class K>
@@ -967,20 +968,15 @@ template <class K>
 bool
 do_intersect(const typename K::Plane_3 &plane,
              const typename K::Ray_3 &ray,
-             const K& k)
+             const K& )
 {
-    typedef typename K::Point_3 Point_3;
+    typename K::Oriented_side_3 oriented_side_3;
 
-    typename Intersection_traits<K, typename K::Plane_3, typename K::Line_3>
+    Oriented_side os = oriented_side_3(plane,ray.source());
-      ::result_type
+    if(os == ON_ORIENTED_BOUNDARY){
-      line_intersection = internal::intersection(plane, ray.supporting_line(), k);
+      return true;
-
+    }
-    if(!line_intersection)
+    return sign(ray.to_vector()* plane.orthogonal_vector()) * os == -1;
-        return false;
-    if(const Point_3 *isp = intersect_get<Point_3>(line_intersection))
-        return ray.collinear_has_on(*isp);
-
-    return true;
 }
 
 
@@ -1642,38 +1638,171 @@ template <class R>
 inline bool
 do_intersect(const Plane_3<R>& plane1, const Plane_3<R>& plane2, const R&)
 {
-  return bool(intersection(plane1, plane2));
+    typedef typename R::RT RT;
+    const RT &a = plane1.a();
+    const RT &b = plane1.b();
+    const RT &c = plane1.c();
+    const RT &d = plane1.d();
+    const RT &p = plane2.a();
+    const RT &q = plane2.b();
+    const RT &r = plane2.c();
+    const RT &s = plane2.d();
+
+    RT det = a*q-p*b;
+    if (det != 0) {
+      return true;
+    }
+    det = a*r-p*c;
+    if (det != 0) {
+      return true;
+    }
+    det = b*r-c*q;
+    if (det != 0) {
+      return true;
+    }
+// degenerate case
+    if (a!=0 || p!=0) {
+      return (a*s == p*d);
+    }
+    if (b!=0 || q!=0) {
+      return (b*s == q*d);
+    }
+    if (c!=0 || r!=0) {
+      return (c*s == r*d);
+    }
+    return true;
 }
 
 
 template <class R>
 inline bool
 do_intersect(const Plane_3<R> &plane1, const Plane_3<R> &plane2,
-             const Plane_3<R> &plane3, const R& r)
+             const Plane_3<R> &plane3, const R&)
 {
-  return bool(intersection(plane1, plane2, plane3, r));
+  typedef typename R::RT RT;
+
+  if(! is_zero(determinant(plane1.a(), plane1.b(), plane1.c(),
+                           plane2.a(), plane2.b(), plane2.c(),
+                           plane3.a(), plane3.b(), plane3.c()))){
+    return true;
+  }
+
+  int pcount = 0;
+  bool b12, b13,b23;
+  if((b12 = parallel(plane1,plane2))) pcount++;
+  if((b13 = parallel(plane1,plane3))) pcount++;
+  if((b23 = parallel(plane2,plane3))) pcount++;
+
+  if(pcount == 3){
+    return (( (plane1 == plane2) || (plane1 == plane2.opposite())) && ( (plane1 == plane3) || (plane1 == plane3.opposite())));
+  }
+
+  if(pcount == 1){
+    if(b12 && ((plane1 == plane2)||(plane1 == plane2.opposite()  ))) return true;
+    if(b13 && ((plane1 == plane3)||(plane1 == plane3.opposite()  ))) return true;
+    if(b23 && ((plane2 == plane3)||(plane2 == plane3.opposite()  ))) return true;
+  }
+
+  int rd = rank_34<RT>(plane1.a(), plane1.b(), plane1.c(), plane1.d(),
+                       plane2.a(), plane2.b(), plane2.c(), plane2.d(),
+                       plane3.a(), plane3.b(), plane3.c(), plane3.d());
+
+  return rd == 2;
 }
 
 
 template <class R>
 inline bool
-do_intersect(const Iso_cuboid_3<R> &i, const Iso_cuboid_3<R> &j, const R&)
+do_intersect(const Iso_cuboid_3<R> &icub1, const Iso_cuboid_3<R> &icub2, const R&)
 {
-  return bool(CGAL::intersection(i, j));
+    typedef typename R::Point_3 Point_3;
+
+    Point_3 min_points[2];
+    Point_3 max_points[2];
+    min_points[0] = (icub1.min)();
+    min_points[1] = (icub2.min)();
+    max_points[0] = (icub1.max)();
+    max_points[1] = (icub2.max)();
+    const int DIM = 3;
+    int min_idx[DIM];
+    int max_idx[DIM];
+    Point_3 newmin;
+    Point_3 newmax;
+    for (int dim = 0; dim < DIM; ++dim) {
+        min_idx[dim] =
+          min_points[0].cartesian(dim) >= min_points[1].cartesian(dim) ? 0 : 1;
+        max_idx[dim] =
+          max_points[0].cartesian(dim) <= max_points[1].cartesian(dim) ? 0 : 1;
+        if (min_idx[dim] != max_idx[dim]
+                && max_points[max_idx[dim]].cartesian(dim)
+                   < min_points[min_idx[dim]].cartesian(dim))
+          return false;
+    }
+    return true;
+
 }
 
 template <class R>
 inline bool
-do_intersect(const Line_3<R> &l, const Iso_cuboid_3<R> &j, const R&)
+do_intersect(const Line_3<R> &line, const Iso_cuboid_3<R> &box, const R&)
 {
-  return bool(CGAL::intersection(l, j));
+    typedef typename R::Point_3 Point_3;
+    typedef typename R::Vector_3 Vector_3;
+    typedef typename R::FT FT;
+    bool all_values = true;
+    FT _min = 0, _max = 0; // initialization to stop compiler warning
+    FT _denum;
+    Point_3 const & _ref_point=line.point();
+    Vector_3 const & _dir=line.direction().vector();
+    Point_3 const & _iso_min=(box.min)();
+    Point_3 const & _iso_max=(box.max)();
+    for (int i=0; i< _ref_point.dimension(); i++) {
+        if (_dir.homogeneous(i) == 0) {
+            if (_ref_point.cartesian(i) < _iso_min.cartesian(i)) {
+              return false;
+            }
+            if (_ref_point.cartesian(i) > _iso_max.cartesian(i)) {
+              return false;
+            }
+        } else {
+            FT newmin, newmax;
+            FT newdenum = _dir.cartesian(i);
+            if (_dir.homogeneous(i) > 0) {
+              newmin = (_iso_min.cartesian(i) - _ref_point.cartesian(i));
+              newmax = (_iso_max.cartesian(i) - _ref_point.cartesian(i));
+            } else {
+              newmin = (_iso_max.cartesian(i) - _ref_point.cartesian(i));
+              newmax = (_iso_min.cartesian(i) - _ref_point.cartesian(i));
+            }
+            if (all_values) {
+                _min = newmin;
+                _max = newmax;
+                _denum = newdenum;
+            } else {
+
+              if (compare_quotients(newmin, newdenum, _min, _denum) == LARGER)
+                    _min = newmin;
+              if (compare_quotients(newmax, newdenum, _max, _denum) == LARGER)
+                    _max = newmax;
+              if (compare_quotients(_max, _denum, _min, _denum) == SMALLER) {
+                  return false;
+                }
+                _denum = newdenum;
+
+            }
+            all_values = false;
+        }
+    }
+    CGAL_kernel_assertion(!all_values);
+    return true;
 }
 
+
 template <class R>
 inline bool
-do_intersect(const Iso_cuboid_3<R> &j, const Line_3<R> &l, const R&)
+do_intersect(const Iso_cuboid_3<R> &j, const Line_3<R> &l, const R& r)
 {
-  return bool(CGAL::intersection(l, j));
+  return do_intersect(l, j, r);
 }
 } // namespace internal
 } // namespace Intersections

Intersections_3/test/Intersections_3/test_intersections_3.cpp

@@ -1,5 +1,8 @@
 // 3D intersection tests.
 
+// We want to check that no division is performed for interface macro Do_intersect_3_RT
+#define CGAL_NO_MPZF_DIVISION_OPERATOR
+
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Homogeneous.h>
 #include <CGAL/MP_Float.h>

Kernel_23/include/CGAL/Kernel/interface_macros.h

@@ -496,8 +496,8 @@ CGAL_Kernel_pred(Counterclockwise_in_between_2,
                  counterclockwise_in_between_2_object)
 CGAL_Kernel_pred(Do_intersect_2,
                  do_intersect_2_object)
-CGAL_Kernel_pred(Do_intersect_3,
+CGAL_Kernel_pred_RT(Do_intersect_3,
-                 do_intersect_3_object)
+                    do_intersect_3_object)
 CGAL_Kernel_pred(Equal_xy_3,
                  equal_xy_3_object)
 CGAL_Kernel_pred(Equal_x_2,

Kernel_23/include/CGAL/Kernel/solve.h

@@ -58,6 +58,37 @@ void solve (const FT &a1, const FT &a2, const FT &a3,
 }
 
 
+template <class FT>
+void solve (const FT &a1, const FT &a2, const FT &a3,
+            const FT &b1, const FT &b2, const FT &b3,
+            const FT &c1, const FT &c2, const FT &c3,
+            const FT &d1, const FT &d2, const FT &d3,
+            FT &x, FT &y, FT &z, FT& denom)
+{
+  FT ab23 = a3*b2 - a2*b3;
+  FT ab13 = a3*b1 - a1*b3;
+  FT ab12 = a2*b1 - a1*b2;
+
+  denom = ab23*c1 - ab13*c2 + ab12*c3;
+
+  FT cd23 = c3*d2 - c2*d3;
+  FT cd13 = c3*d1 - c1*d3;
+  FT cd12 = c2*d1 - c1*d2;
+
+  x = (b3*cd12 - b2*cd13 + b1*cd23);
+
+  y = (a2*cd13 - cd12*a3 - cd23*a1);
+
+  z = (ab23*d1 + ab12*d3 - ab13*d2);
+  if(denom < 0){
+    denom = -denom;
+    x = -x;
+    y = -y;
+    z = -z;
+  }
+}
+
+
 // this is for a parabola c1, c2, c3 are equal to 1
 template <class FT>
 void solve_quadratic (const FT &a1, const FT &a2, const FT &a3,

Kernel_23/include/CGAL/rank.h

@@ -0,0 +1,164 @@
+// Copyright (c) 2021
+// GeometryFactory (France),
+//
+// This file is part of CGAL (www.cgal.org)
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Sebastien Loriot
+
+#ifndef CGAL_RANK_H
+#define CGAL_RANK_H
+
+namespace CGAL {
+
+template <class RT>
+int rank_11(const RT& a0)
+{
+  return a0!=0 ? 1 : 0;
+}
+
+template <class RT>
+int rank_21(const RT& a0, const RT& a1)
+{
+  return (a0!=0 || a1 !=0) ? 1 : 0;
+}
+
+template <class RT>
+int rank_12(const RT& a0, const RT& a1)
+{
+  return (a0!=0 || a1 !=0) ? 1 : 0;
+}
+
+template <class RT>
+int rank_31(const RT& a0, const RT& a1, const RT& a2)
+{
+  return (a0!=0 || a1 !=0 || a2 !=0) ? 1 : 0;
+}
+
+template <class RT>
+int rank_32(const RT& a0, const RT& b0,
+            const RT& a1, const RT& b1,
+            const RT& a2, const RT& b2)
+{
+  if (a0==0)
+  {
+    if (a1==0)
+    {
+      if (a2==0)
+      {
+        return rank_31<RT>(b0,b1,b2);
+      }
+      else
+      {
+        return 1 + rank_21<RT>(b0, b1);
+      }
+    }
+    else
+    {
+      return 1 + rank_21<RT>(b0, a1*b2-a2*b1);
+    }
+  }
+  else
+  {
+    return 1 + rank_21<RT>(a0*b1-a1*b0, a0*b2-a2*b0);
+  }
+}
+
+template <class RT>
+int rank_22(const RT& a0, const RT& b0,
+            const RT& a1, const RT& b1)
+{
+  if (a0==0)
+  {
+    if (a1==0)
+    {
+      return rank_21<RT>(b0,b1);
+    }
+    else
+      return 1 + rank_11<RT>(b0);
+  }
+  return 1 + rank_11<RT>(a0*b1-a1*b0);
+}
+
+template <class RT>
+int rank_33(const RT& a0, const RT& b0, const RT& c0,
+            const RT& a1, const RT& b1, const RT& c1,
+            const RT& a2, const RT& b2, const RT& c2)
+{
+  if (a0==0)
+  {
+    if (a1==0)
+    {
+      if (a2==0)
+      {
+        return rank_32<RT>(b0, c0, b1, c1, b2, c2);
+      }
+      else
+      {
+        return 1 + rank_22<RT>(b0, c0, b1, c1);
+      }
+    }
+    else
+      return 1 + rank_22<RT>(b0, c0, a1*b2-a2*b1, a1*c2-a2*c1);
+  }
+  else
+  {
+    return 1 + rank_22<RT>(a0*b1-a1*b0, a0*c1-a1*c0, a0*b2-a2*b0, a0*c2-a2*c0);
+  }
+}
+
+
+template <class RT>
+int rank_23(const RT& a0, const RT& b0, const RT& c0,
+            const RT& a1, const RT& b1, const RT& c1)
+{
+  if (a0==0)
+  {
+    if (a1==0)
+    {
+      return rank_22<RT>(b0, c0, b1, c1);
+    }
+    else
+      return 1 + rank_12<RT>(b0,c0);
+  }
+  else
+  {
+    return 1 + rank_12<RT>(a0*b1-a1*b0,a0*c1-a1*c0);
+  }
+}
+
+
+template <class RT>
+int rank_34(const RT& a0, const RT& b0, const RT& c0, const RT& d0,
+            const RT& a1, const RT& b1, const RT& c1, const RT& d1,
+            const RT& a2, const RT& b2, const RT& c2, const RT& d2)
+{
+  if (a0==0)
+  {
+    if (a1==0)
+    {
+      if (a2==0)
+      {
+        return rank_33<RT>(b0, c0, d0, b1, c1, d1, b2, c2, d2);
+      }
+      else
+      {
+        return 1 + rank_23<RT>(b0, c0, d0, b1, c1, d1);
+      }
+    }
+    else
+      return 1 + rank_23<RT>(b0, c0, d0,a1*b2-a2*b1, a1*c2-a2*c1, a1*d2 - a2*d1);
+  }
+  else
+  {
+    return 1 + rank_23<RT>(a0*b1-a1*b0, a0*c1-a1*c0, a0*d1-a1*d0, a0*b2-a2*b0, a0*c2-a2*c0, a0*d2-a2*d0);
+  }
+}
+
+} // namespace CGAL
+
+#endif // CGAL_RANK_H

Kernel_23/test/Kernel_23/rank.cpp

@@ -0,0 +1,115 @@
+#include <CGAL/rank.h>
+#include <cassert>
+#include <iostream>
+
+typedef int FT;
+
+
+void test_rank_33(FT a0, FT b0, FT c0,
+                  FT a1, FT b1, FT c1,
+                  FT a2, FT b2, FT c2,
+                  int expected)
+{
+  std::cout << "testing:\n"
+            << "\t " << a0 << "\t" << b0 << "\t" << c0 << "\n"
+            << "\t " << a1 << "\t" << b1 << "\t" << c1 << "\n"
+            << "\t " << a2 << "\t" << b2 << "\t" << c2 << "\n";
+
+  assert(CGAL::rank_33<FT>(a0,b0,c0,a1,b1,c1,a2,b2,c2)==expected);
+  assert(CGAL::rank_33<FT>(a0,b0,c0,a2,b2,c2,a1,b1,c1)==expected);
+  assert(CGAL::rank_33<FT>(a1,b1,c1,a0,b0,c0,a2,b2,c2)==expected);
+  assert(CGAL::rank_33<FT>(a1,b1,c1,a2,b2,c2,a0,b0,c0)==expected);
+  assert(CGAL::rank_33<FT>(a2,b2,c2,a0,b0,c0,a1,b1,c1)==expected);
+  assert(CGAL::rank_33<FT>(a2,b2,c2,a1,b1,c1,a0,b0,c0)==expected);
+}
+
+void test_rank_34(FT a0, FT b0, FT c0, FT d0,
+                  FT a1, FT b1, FT c1, FT d1,
+                  FT a2, FT b2, FT c2, FT d2,
+                  int expected)
+{
+  std::cout << "testing:\n"
+            << "\t " << a0 << "\t" << b0 << "\t" << c0 << "\t" << d0 << "\n"
+            << "\t " << a1 << "\t" << b1 << "\t" << c1 << "\t" << d1 << "\n"
+            << "\t " << a2 << "\t" << b2 << "\t" << c2 << "\t" << d2 << "\n";
+
+  assert(CGAL::rank_34<FT>(a0,b0,c0,d0,a1,b1,c1,d1,a2,b2,c2,d2)==expected);
+  assert(CGAL::rank_34<FT>(a0,b0,c0,d0,a2,b2,c2,d2,a1,b1,c1,d1)==expected);
+  assert(CGAL::rank_34<FT>(a1,b1,c1,d1,a0,b0,c0,d0,a2,b2,c2,d2)==expected);
+  assert(CGAL::rank_34<FT>(a1,b1,c1,d1,a2,b2,c2,d2,a0,b0,c0,d0)==expected);
+  assert(CGAL::rank_34<FT>(a2,b2,c2,d2,a0,b0,c0,d0,a1,b1,c1,d1)==expected);
+  assert(CGAL::rank_34<FT>(a2,b2,c2,d2,a1,b1,c1,d1,a0,b0,c0,d0)==expected);
+}
+
+
+
+int main()
+{
+  test_rank_33(1,0,0,
+               0,1,0,
+               0,0,1, 3);
+  test_rank_33(1,0,0,
+               0,1,0,
+               1,0,1, 3);
+  test_rank_33(1,0,0,
+               1,1,1,
+               1,0,1, 3);
+  test_rank_33(1,0,0,
+               1,1,1,
+               1,0,1, 3);
+  test_rank_33(1,0,0,
+               1,1,0,
+               1,5,0, 2);
+  test_rank_33(0,0,0,
+               0,0,0,
+               0,0,0, 0);
+  test_rank_33(1,2,3,
+               1,2,3,
+               1,2,3, 1);
+  test_rank_33(1,2,3,
+               2,4,6,
+               4,8,12, 1);
+  test_rank_33(1,2,3,
+               2,4,6,
+               4,8,11, 2);
+  test_rank_33(1,0,1,
+               1,0,1,
+               1,0,1, 1);
+  test_rank_33(1,1,1,
+               1,0,1,
+               1,0,1, 2);
+
+  test_rank_33(0,1,1,
+               0,0,1,
+               0,0,1, 2);
+
+  test_rank_34(1,0,0,0,
+               0,1,0,0,
+               0,0,1,0, 3);
+
+  test_rank_34(1,0,0,0,
+               0,1,0,0,
+               1,1,0,0, 2);
+
+  test_rank_34(1,0,0,0,
+               0,1,0,0,
+               1,1,0,1, 3);
+
+  test_rank_34(0,1,0,0,
+               0,0,1,0,
+               0,1,1,0, 2);
+
+  test_rank_34(0,1,0,0,
+               0,0,1,0,
+               0,1,1,1, 3);
+
+  test_rank_34(1,0,0,0,
+               0,0,1,0,
+               1,0,1,0, 2);
+
+  test_rank_34(1,0,0,0,
+               0,0,1,0,
+               1,0,1,1, 3);
+
+  return 0;
+}

NewKernel_d/include/CGAL/NewKernel_d/utils.h

@@ -94,25 +94,6 @@ struct Has_type_different_from <T, No, true>
                         enum { value=true };
         };
 
-        // like std::multiplies but allows mixing types
-        // in C++11 in doesn't need to be a template
-        template < class Ret >
-        struct multiplies {
-                template<class A,class B>
-                decltype(auto) operator()(A&&a,B&&b)const
-                {
-                        return std::forward<A>(a)*std::forward<B>(b);
-                }
-        };
-        template < class Ret >
-        struct division {
-                template<class A,class B>
-                decltype(auto) operator()(A&&a,B&&b)const
-                {
-                        return std::forward<A>(a)/std::forward<B>(b);
-                }
-        };
-
         using std::decay;
 
         template<class T,class U> struct Type_copy_ref { typedef U type; };

Number_types/include/CGAL/Mpzf.h

@@ -10,6 +10,10 @@
 //
 // Author(s)        :  Marc Glisse
 
+#ifndef CGAL_NO_MPZF_DIVISION_OPERATOR
+#define CGAL_MPZF_DIVISION_OPERATOR 1
+#endif
+
 #ifndef CGAL_MPZF_H
 #define CGAL_MPZF_H
 #include <cstdlib>
@@ -774,7 +778,11 @@ struct Mpzf {
     return res;
   }
 
+#ifndef CGAL_MPZF_DIVISION_OPERATOR
+  friend Mpzf division(Mpzf const&a, Mpzf const&b){
+#else // CGAL_MPZF_DIVISION_OPERATOR
   friend Mpzf operator/(Mpzf const&a, Mpzf const&b){
+#endif // CGAL_MPZF_DIVISION_OPERATOR
     // FIXME: Untested
     int asize=std::abs(a.size);
     int bsize=std::abs(b.size);
@@ -909,7 +917,11 @@ struct Mpzf {
   Mpzf& operator+=(Mpzf const&x){ *this=*this+x; return *this; }
   Mpzf& operator-=(Mpzf const&x){ *this=*this-x; return *this; }
   Mpzf& operator*=(Mpzf const&x){ *this=*this*x; return *this; }
+#ifdef CGAL_MPZF_DIVISION_OPERATOR
   Mpzf& operator/=(Mpzf const&x){ *this=*this/x; return *this; }
+#else // not CGAL_MPZF_DIVISION_OPERATOR
+  Mpzf& operator/=(Mpzf const&x){ *this=division(*this,x); return *this; }
+#endif // not CGAL_MPZF_DIVISION_OPERATOR
 
   bool is_canonical () const {
     if (size == 0) return true;
@@ -1096,7 +1108,11 @@ std::istream& operator>> (std::istream& is, Mpzf& a)
           Type operator()(
               const Type& x,
               const Type& y ) const {
+#ifdef CGAL_MPZF_DIVISION_OPERATOR
             return x / y;
+#else // not CGAL_MPZF_DIVISION_OPERATOR
+            return division(x, y);
+#endif // not CGAL_MPZF_DIVISION_OPERATOR
           }
         };