Merge pull request #5937 from danston/Number_types-add_benchmarks-danston

Robust Quotient with cpp_int

Arithmetic_kernel/include/CGAL/BOOST_MP_arithmetic_kernel.h

@@ -9,8 +9,8 @@
 //
 // Author: Marc Glisse <marc.glisse@inria.fr>
 
-#ifndef CGAL_GMPXX_ARITHMETIC_KERNEL_H
+#ifndef CGAL_BOOST_MP_ARITHMETIC_KERNEL_H
-#define CGAL_GMPXX_ARITHMETIC_KERNEL_H
+#define CGAL_BOOST_MP_ARITHMETIC_KERNEL_H
 
 #include <CGAL/Arithmetic_kernel/Arithmetic_kernel_base.h>
 #include <CGAL/Get_arithmetic_kernel.h>

NewKernel_d/test/NewKernel_d/Epick_d.cpp

@@ -1,6 +1,5 @@
 #include <boost/config.hpp>
 
-//#define BOOST_RESULT_OF_USE_DECLTYPE 1
 #include <CGAL/Epick_d.h>
 #include <CGAL/Epeck_d.h>
 #include <typeinfo>

Number_types/include/CGAL/Exact_integer.h

@@ -50,30 +50,30 @@ typedef unspecified_type Exact_integer;
 
 #else // not DOXYGEN_RUNNING
 
-#if CGAL_USE_GMPXX
+#if ( (defined(CGAL_TEST_SUITE) && CGAL_VERSION_NR == 1050500900) || defined(CGAL_FORCE_USE_BOOST_MP))\
-
+    && BOOST_VERSION > 107800 && defined(CGAL_USE_BOOST_MP)
+// use boost-mp by default in the testsuite until 5.5-beta is out
+typedef BOOST_cpp_arithmetic_kernel::Integer Exact_integer;
+#else // BOOST_VERSION > 107800
+#ifdef CGAL_USE_GMPXX
 typedef mpz_class Exact_integer;
-
+#elif defined(CGAL_USE_GMP)
-#elif CGAL_USE_GMP
+#if defined(CGAL_USE_BOOST_MP)
-# ifdef CGAL_USE_BOOST_MP
+typedef BOOST_gmp_arithmetic_kernel::Integer Exact_integer;
-typedef boost::multiprecision::mpz_int Exact_integer;
+#else
-# else
 typedef Gmpz Exact_integer;
-# endif
+#endif
-
+#elif defined(CGAL_USE_LEDA)
-#elif CGAL_USE_LEDA
-
 typedef leda_integer Exact_integer;
-
+#elif defined(CGAL_USE_BOOST_MP)
-#elif CGAL_USE_CORE
+typedef BOOST_cpp_arithmetic_kernel::Integer Exact_integer;
-
+#elif defined(CGAL_USE_CORE)
 typedef CORE::BigInt Exact_integer;
-
+#else
-#elif defined CGAL_USE_BOOST_MP
+#error "ERROR: Cannot determine a BigInt type!"
-
-typedef boost::multiprecision::cpp_int Exact_integer;
-
 #endif // CGAL_USE_CORE
+#endif // BOOST_VERSION > 107800
+
 #endif // not DOXYGEN_RUNNING
 
 } /* end namespace CGAL */

Number_types/include/CGAL/Number_types/internal/Exact_type_selector.h

@@ -40,7 +40,7 @@
 #  include <CGAL/leda_real.h>
 #endif
 #ifdef CGAL_USE_CORE
-// #  include <CGAL/CORE_Expr.h>
+// # include <CGAL/CORE_Expr.h>
 namespace CORE {
 class Expr;
 }
@@ -54,23 +54,39 @@ namespace CGAL { namespace internal {
 // It should support the built-in types.
 template < typename >
 struct Exact_field_selector
+
+#if ( (defined(CGAL_TEST_SUITE) && CGAL_VERSION_NR == 1050500900) || defined(CGAL_FORCE_USE_BOOST_MP))\
+    && BOOST_VERSION > 107800 && defined(CGAL_USE_BOOST_MP)
+// use boost-mp by default in the testsuite until 5.5-beta is out
+// Boost
+{ typedef BOOST_cpp_arithmetic_kernel::Rational Type; };
+#else // BOOST_VERSION > 107800
 #ifdef CGAL_USE_GMPXX
 { typedef mpq_class Type; };
 #elif defined(CGAL_USE_GMP)
-# if defined(CGAL_USE_BOOST_MP)
+#if defined(CGAL_USE_BOOST_MP)
-{ typedef boost::multiprecision::mpq_rational Type; };
+{ typedef BOOST_gmp_arithmetic_kernel::Rational Type; };
-# else
+#else
 { typedef Gmpq Type; };
-# endif
+#endif
 #elif defined(CGAL_USE_LEDA)
 { typedef leda_rational Type; };
-#elif 0 && defined(CGAL_USE_BOOST_MP)
+#elif defined(CGAL_USE_BOOST_MP)
 // See the discussion in https://github.com/CGAL/cgal/pull/3614
 // This is disabled for now because cpp_rational is even slower than Quotient<MP_Float>. Quotient<cpp_int> will be a good candidate after some polishing.
+// In fact, the new version of cpp_rational from here: https://github.com/boostorg/multiprecision/pull/366
+// is much better than Quotient<cpp_int> because it is using smart gcd and is well-supported
+// while Quotient does not. Though, we can still use it if needed.
+#if BOOST_VERSION <= 107800
+// See this comment: https://github.com/CGAL/cgal/pull/5937#discussion_r721533675
+{ typedef Quotient<boost::multiprecision::cpp_int> Type; };
+#else
 { typedef BOOST_cpp_arithmetic_kernel::Rational Type; };
+#endif
 #else
 { typedef Quotient<MP_Float> Type; };
 #endif
+#endif // BOOST_VERSION > 107800
 
 // By default, a field is a safe choice of ring.
 template < typename T >

Number_types/include/CGAL/Quotient.h

@@ -33,6 +33,9 @@
 #include <CGAL/Kernel/mpl.h>
 
 #include <boost/operators.hpp>
+#ifdef CGAL_USE_BOOST_MP
+#include <boost/multiprecision/cpp_int.hpp>
+#endif // CGAL_USE_BOOST_MP
 
 namespace CGAL {
 
@@ -46,6 +49,15 @@ template < typename NT >
 inline void
 simplify_quotient(NT &, NT &) {}
 
+#ifdef CGAL_USE_BOOST_MP
+inline void
+simplify_quotient(boost::multiprecision::cpp_int & a, boost::multiprecision::cpp_int & b) {
+  const boost::multiprecision::cpp_int r = boost::multiprecision::gcd(a, b);
+  a /= r;
+  b /= r;
+}
+#endif // CGAL_USE_BOOST_MP
+
 // This one should be replaced by some functor or tag.
 // Meanwhile, the class is specialized for Gmpz, mpz_class, leda_integer.
 template < typename NT >
@@ -687,7 +699,7 @@ template < class NT > class Real_embeddable_traits_quotient_base< Quotient<NT> >
       : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
       public:
         std::pair<double, double> operator()( const Type& x ) const {
-          Interval_nt<> quot =
+          const Interval_nt<> quot =
                           Interval_nt<>(CGAL_NTS to_interval(x.numerator())) /
                           Interval_nt<>(CGAL_NTS to_interval(x.denominator()));
           return std::make_pair(quot.inf(), quot.sup());

Number_types/include/CGAL/boost_mp.h

@@ -13,6 +13,7 @@
 #define CGAL_BOOST_MP_H
 
 #include <CGAL/config.h>
+#include <CGAL/number_utils.h>
 
 // It is easier to disable this number type completely for old versions.
 // Before 1.63, I/O is broken.  Again, disabling the whole file is just the
@@ -23,6 +24,7 @@
     (!defined _MSC_VER || BOOST_VERSION >= 107000)
 #define CGAL_USE_BOOST_MP 1
 
+#include <CGAL/Quotient.h>
 #include <CGAL/functional.h> // *ary_function
 #include <CGAL/number_type_basic.h>
 #include <CGAL/Modular_traits.h>
@@ -139,6 +141,227 @@ struct Algebraic_structure_traits<boost::multiprecision::detail::expression<T1,T
 
 // Real_embeddable_traits
 
+namespace Boost_MP_internal {
+
+  // here we know that `intv` contains int64 numbers such that their msb is std::numeric_limits<double>::digits-1
+  // TODO: possibly return denormals sometimes...
+  inline
+  std::pair<double,double> shift_positive_interval( const std::pair<double,double>& intv, const int e ) {
+    CGAL_assertion(intv.first > 0.0);
+    CGAL_assertion(intv.second > 0.0);
+
+    CGAL_assertion_code(
+    union {
+#ifdef CGAL_LITTLE_ENDIAN
+      struct { uint64_t man:52; uint64_t exp:11; uint64_t sig:1; } s;
+#else /* CGAL_BIG_ENDIAN */
+      //WARNING: untested!
+      struct { uint64_t sig:1; uint64_t exp:11; uint64_t man:52; } s;
+#endif
+      double d;
+    } conv;
+    conv.d = intv.first;
+    )
+    // Check that the exponent of intv.inf is 52, which corresponds to a 53 bit integer
+    CGAL_assertion(conv.s.exp - ((1 << (11 - 1)) - 1) == std::numeric_limits<double>::digits - 1);
+
+    typedef std::numeric_limits<double> limits;
+
+    // warning: min_exponent and max_exponent are 1 more than what the name suggests
+    if (e < limits::min_exponent - limits::digits)
+      return {0, (limits::min)()};
+    if (e > limits::max_exponent - limits::digits)
+      return {(limits::max)(), limits::infinity()}; // intv is positive
+
+    const double scale = std::ldexp(1.0, e); // ldexp call is exact
+    return { scale * intv.first, scale * intv.second }; // cases that would require a rounding mode have been handled above
+  }
+
+  // This function checks if the computed interval is correct and if it is tight.
+  template<typename Type>
+  bool are_bounds_correct( const double l, const double u, const Type& x ) {
+    typedef std::numeric_limits<double> limits;
+
+    const double inf = std::numeric_limits<double>::infinity();
+    if ( u!=l && (l==-inf || u==inf
+                          || (u==0 && l >=  -(limits::min)())
+                          || (l==0 && u <= (limits::min)())) )
+    {
+      return x >  Type((limits::max)()) ||
+             x < Type(-(limits::max)()) ||
+             (x > Type(-(limits::min)()) && x < Type((limits::min)()));
+    }
+
+    if (!(u == l || u == std::nextafter(l, +inf))) return false;
+    //TODO: Type(nextafter(l,inf))>x && Type(nextafter(u,-inf))<x
+
+    const Type lb(l), ub(u);
+    const bool are_bounds_respected = (lb <= x && x <= ub);
+    return are_bounds_respected;
+  }
+
+  // This one returns zero length interval that is inf = sup.
+  inline
+  std::pair<double, double> get_0ulp_interval( const int shift, const uint64_t p ) {
+
+    const double pp_dbl = static_cast<double>(p);
+    const std::pair<double,double> intv(pp_dbl, pp_dbl);
+
+    return shift_positive_interval(intv, -shift);
+  }
+
+  // This one returns 1 unit length interval.
+  inline
+  std::pair<double, double> get_1ulp_interval( const int shift, const uint64_t p ) {
+
+    const double pp_dbl = static_cast<double>(p);
+    const double qq_dbl = pp_dbl+1;
+    const std::pair<double,double> intv(pp_dbl, qq_dbl);
+    return shift_positive_interval(intv, -shift);
+  }
+
+  template<typename ET>
+  std::pair<double, double> to_interval( ET x, int extra_shift = 0 );
+
+  // This is a version of to_interval that converts a rational type into a
+  // double tight interval.
+  template<typename Type, typename ET>
+  std::pair<double, double> to_interval( ET xnum, ET xden ) {
+
+    CGAL_assertion(!CGAL::is_zero(xden));
+    CGAL_assertion_code(const Type input(xnum, xden));
+    double l = 0.0, u = 0.0;
+    if (CGAL::is_zero(xnum)) { // return [0.0, 0.0]
+      CGAL_assertion(are_bounds_correct(l, u, input));
+      return std::make_pair(l, u);
+    }
+    CGAL_assertion(!CGAL::is_zero(xnum));
+
+    // Handle signs.
+    bool change_sign = false;
+    const bool is_num_pos = CGAL::is_positive(xnum);
+    const bool is_den_pos = CGAL::is_positive(xden);
+    if (!is_num_pos && !is_den_pos) {
+      xnum = -xnum;
+      xden = -xden;
+    } else if (!is_num_pos && is_den_pos) {
+      change_sign = true;
+      xnum = -xnum;
+    } else if (is_num_pos && !is_den_pos) {
+      change_sign = true;
+      xden = -xden;
+    }
+    CGAL_assertion(CGAL::is_positive(xnum) && CGAL::is_positive(xden));
+
+    const int64_t num_dbl_digits = std::numeric_limits<double>::digits - 1;
+    const int64_t msb_num = static_cast<int64_t>(boost::multiprecision::msb(xnum));
+    const int64_t msb_den = static_cast<int64_t>(boost::multiprecision::msb(xden));
+
+#if 0 // Optimisation for the case of input that are double
+    // An alternative strategy would be to convert numerator and denominator to
+    // intervals, then divide. However, this would require setting the rounding
+    // mode (and dividing intervals is not completely free). An important
+    // special case is when the rational is exactly equal to a double
+    // (fit_in_double). Then the denominator is a power of 2, so we can skip
+    // the division and it becomes unnecessary to set the rounding mode, we
+    // just need to modify the exponent correction for the denominator.
+    if(msb_den == static_cast<int64_t>(lsb(xden))) {
+      std::tie(l,u)=to_interval(xnum, msb_den);
+      if (change_sign) {
+        CGAL_assertion(are_bounds_correct(-u, -l, input));
+        return {-u, -l};
+      }
+      CGAL_assertion(are_bounds_correct(l, u, input));
+      return {u, l};
+    }
+#endif
+
+    const int64_t msb_diff = msb_num - msb_den;
+    // Shift so the division result has at least 53 (and at most 54) bits
+    int shift = static_cast<int>(num_dbl_digits - msb_diff + 1);
+    CGAL_assertion(shift == num_dbl_digits - msb_diff + 1);
+
+    if (shift > 0) {
+      xnum <<= +shift;
+    } else if (shift < 0) {
+      xden <<= -shift;
+    }
+    CGAL_assertion(num_dbl_digits + 1 ==
+      static_cast<int64_t>(boost::multiprecision::msb(xnum)) -
+      static_cast<int64_t>(boost::multiprecision::msb(xden)));
+
+    ET p, r;
+    boost::multiprecision::divide_qr(xnum, xden, p, r);
+    uint64_t uip = static_cast<uint64_t>(p);
+    const int64_t p_bits = static_cast<int64_t>(boost::multiprecision::msb(p));
+    bool exact = r.is_zero();
+
+    if (p_bits > num_dbl_digits) { // case 54 bits
+      exact &= ((uip & 1) == 0);
+      uip>>=1;
+      --shift;
+    }
+    std::tie(l, u) = exact ? get_0ulp_interval(shift, uip) : get_1ulp_interval(shift, uip);
+
+    if (change_sign) {
+      const double t = l;
+      l = -u;
+      u = -t;
+    }
+
+    CGAL_assertion(are_bounds_correct(l, u, input));
+    return std::make_pair(l, u);
+  }
+
+  // This is a version of to_interval that converts an integer type into a
+  // double tight interval.
+  template<typename ET>
+  std::pair<double, double> to_interval( ET x, int extra_shift) {
+
+    CGAL_assertion_code(const ET input = x);
+    double l = 0.0, u = 0.0;
+    if (CGAL::is_zero(x)) { // return [0.0, 0.0]
+      CGAL_assertion(are_bounds_correct(l, u, input));
+      return std::make_pair(l, u);
+    }
+    CGAL_assertion(!CGAL::is_zero(x));
+
+    bool change_sign = false;
+    const bool is_pos = CGAL::is_positive(x);
+    if (!is_pos) {
+      change_sign = true;
+      x = -x;
+    }
+    CGAL_assertion(CGAL::is_positive(x));
+
+    const int64_t n = static_cast<int64_t>(boost::multiprecision::msb(x)) + 1;
+    const int64_t num_dbl_digits = std::numeric_limits<double>::digits;
+
+    if (n > num_dbl_digits) {
+      const int64_t mindig = static_cast<int64_t>(boost::multiprecision::lsb(x));
+      int e = static_cast<int>(n - num_dbl_digits);
+      x >>= e;
+      if (n - mindig > num_dbl_digits)
+        std::tie(l, u) = get_1ulp_interval(-e+extra_shift, static_cast<uint64_t>(x));
+      else
+        std::tie(l, u) = get_0ulp_interval(-e+extra_shift, static_cast<uint64_t>(x));
+    } else {
+      l = u = extra_shift==0 ? static_cast<double>(static_cast<uint64_t>(x))
+                             : std::ldexp(static_cast<double>(static_cast<uint64_t>(x)),-extra_shift);
+    }
+
+    if (change_sign) {
+      const double t = l;
+      l = -u;
+      u = -t;
+    }
+
+    CGAL_assertion(extra_shift != 0 || are_bounds_correct(l, u, input));
+    return std::make_pair(l, u);
+  }
+
+} // Boost_MP_internal
+
 template <class NT>
 struct RET_boost_mp_base
     : public INTERN_RET::Real_embeddable_traits_base< NT , CGAL::Tag_true > {
@@ -192,24 +415,41 @@ struct RET_boost_mp_base
 
     struct To_interval
         : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
+
         std::pair<double, double>
         operator()(const Type& x) const {
-            // See if https://github.com/boostorg/multiprecision/issues/108 suggests anything better
+
-            // assume the conversion is within 1 ulp
+          // See if https://github.com/boostorg/multiprecision/issues/108 suggests anything better
-            // adding IA::smallest() doesn't work because inf-e=inf, even rounded down.
+          // assume the conversion is within 1 ulp
-            double i = x.template convert_to<double>();
+          // adding IA::smallest() doesn't work because inf-e=inf, even rounded down.
-            double s = i;
+
-            double inf = std::numeric_limits<double>::infinity();
+          // We must use to_nearest here.
-            int cmp = x.compare(i);
+          double i;
-            if (cmp > 0) {
+          const double inf = std::numeric_limits<double>::infinity();
-              s = nextafter(s, +inf);
+          {
-              CGAL_assertion(x.compare(s) < 0);
+            Protect_FPU_rounding<true> P(CGAL_FE_TONEAREST);
+            i = static_cast<double>(x);
+            if (i == +inf) {
+              return std::make_pair((std::numeric_limits<double>::max)(), i);
+            } else if (i == -inf) {
+              return std::make_pair(i, std::numeric_limits<double>::lowest());
             }
-            else if (cmp < 0) {
+          }
-              i = nextafter(i, -inf);
+          double s = i;
-              CGAL_assertion(x.compare(i) > 0);
+          CGAL_assertion(CGAL::abs(i) != inf && CGAL::abs(s) != inf);
-            }
+
-            return std::pair<double, double> (i, s);
+          // Throws uncaught exception: Cannot convert a non-finite number to an integer.
+          // We can catch it earlier by using the CGAL_assertion() one line above.
+          const int cmp = x.compare(i);
+          if (cmp > 0) {
+            s = nextafter(s, +inf);
+            CGAL_assertion(x.compare(s) < 0);
+          }
+          else if (cmp < 0) {
+            i = nextafter(i, -inf);
+            CGAL_assertion(x.compare(i) > 0);
+          }
+          return std::pair<double, double>(i, s);
         }
     };
 };
@@ -219,11 +459,30 @@ struct RET_boost_mp;
 
 template <class NT>
 struct RET_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_integer> >
-    : RET_boost_mp_base <NT> {};
+    : RET_boost_mp_base <NT> {
+    typedef NT Type;
+    struct To_interval
+        : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
+
+        std::pair<double, double> operator()( const Type& x ) const {
+          return Boost_MP_internal::to_interval(x);
+        }
+    };
+};
 
 template <class NT>
 struct RET_boost_mp <NT, boost::mpl::int_<boost::multiprecision::number_kind_rational> >
-    : RET_boost_mp_base <NT> {};
+    : RET_boost_mp_base <NT> {
+    typedef NT Type;
+    struct To_interval
+        : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
+
+        std::pair<double, double> operator()( const Type& x ) const {
+          return Boost_MP_internal::to_interval<Type>(
+            boost::multiprecision::numerator(x), boost::multiprecision::denominator(x));
+        }
+    };
+};
 
 #ifdef CGAL_USE_MPFR
 // Because of these full specializations, things get instantiated more eagerly. Make them artificially partial if necessary.
@@ -607,6 +866,25 @@ namespace internal {
 #endif
 } // namespace internal
 
+#ifdef CGAL_USE_BOOST_MP
+
+template< > class Real_embeddable_traits< Quotient<boost::multiprecision::cpp_int> >
+    : public INTERN_QUOTIENT::Real_embeddable_traits_quotient_base< Quotient<boost::multiprecision::cpp_int> > {
+
+  public:
+    typedef Quotient<boost::multiprecision::cpp_int> Type;
+
+    class To_interval
+      : public CGAL::cpp98::unary_function< Type, std::pair< double, double > > {
+      public:
+        std::pair<double, double> operator()( const Type& x ) const {
+          return Boost_MP_internal::to_interval<Type>(x.num, x.den);
+        }
+    };
+};
+
+#endif // CGAL_USE_BOOST_MP
+
 } //namespace CGAL
 
 #include <CGAL/BOOST_MP_arithmetic_kernel.h>

Number_types/test/Number_types/CMakeLists.txt

@@ -88,3 +88,8 @@ endif()#NOT CGAL_DISABLE_GMP
 # all the programs below will be linked against MPFI in case it is present
 create_single_source_cgal_program("Quotient_new.cpp")
 create_single_source_cgal_program("test_nt_Coercion_traits.cpp")
+
+find_package(Boost)
+if(Boost_FOUND)
+create_single_source_cgal_program("to_interval_test_boost.cpp")
+endif()

Number_types/test/Number_types/to_interval_test_boost.cpp

@@ -0,0 +1,867 @@
+// STL.
+#include <array>
+#include <limits>
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <vector>
+
+// Boost.
+#include <boost/type_index.hpp>
+
+// Kernels.
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+// CGAL.
+#include <CGAL/Random.h>
+#include <CGAL/Quotient.h>
+#include <CGAL/Exact_integer.h>
+#include <CGAL/Exact_rational.h>
+#include <CGAL/Lazy_exact_nt.h>
+
+#ifdef CGAL_USE_BOOST_MP
+
+#if false // https://github.com/boostorg/multiprecision/issues/370
+
+void test_boost_eval_lehmer() {
+
+  const boost::multiprecision::cpp_int a("500000000000000052504760255204421627079393309355027816932345132815919505535709229444276879024105562954502314530690391078574434507015318513443905076213688875017942541908041275407131568575177172639474548726709751235383681696449966404295647940685784470144122251803020020951078103818191513659921807053133698549053838430992170843235673537548059987539601671975279280846041564435631581262016246808786828637048154067265620710396778995313534536353760281048487250661054626168637371167135426013683337484254647996964562455566714879467026196409913165805735073230830136024016362543811769017875638974011487488573436");
+  const boost::multiprecision::cpp_int b("1500000000000000157514280765613264881238179928065083450797035398447758516607127688332830637072316688863506943592071173235723303521045955540331715228641066625053827625724123826221394705725531517918423646180129253706151045089349899212886943822057353410432366755409060062853234311454574540979765421159386595647161515292215193506006556519037965168192736708179557957863203557666055574947146355487693991882510747766220045897624670399027877365714431356466054500731862264092476764347207739651025585146903094168986610767496468412336047796468657032646893153521091155634158263410282629846280069312485301157888001");
+  const boost::multiprecision::cpp_int r = boost::multiprecision::gcd(a, b);
+}
+
+#endif // issue 370
+
+void test_minimal_boost_gcd() {
+
+  boost::multiprecision::cpp_int u = 1;
+  for (unsigned i = 1; i <= 50; ++i) {
+    u *= i;
+  }
+  std::cout << "u: " << u << std::endl;
+
+  boost::multiprecision::cpp_int v = 1;
+  for (unsigned i = 1; i <= 100; ++i) {
+    v *= i;
+  }
+  std::cout << "v: " << v << std::endl;
+
+  // const auto r = boost::multiprecision::gcd(u, v); // fail
+  const boost::multiprecision::cpp_int r = boost::multiprecision::gcd(u, v); // pass
+  std::cout << "r: " << r << std::endl;
+
+  u = u / r;
+  v = v / r;
+
+  std::cout << "new u: " << u << std::endl;
+  std::cout << "new v: " << v << std::endl;
+}
+
+#if false // _MM_ROUND_UP is not available on all platforms
+
+void test_minimal_nextafter() {
+
+  _MM_SET_ROUNDING_MODE(_MM_ROUND_UP); // fail
+  // _MM_SET_ROUNDING_MODE(_MM_ROUND_NEAREST); // pass
+
+  const boost::multiprecision::cpp_int x("1312729512902970206056841780066779136");
+
+  double i = x.template convert_to<double>();
+  double s = i;
+
+  const double inf = std::numeric_limits<double>::infinity();
+  assert(i != inf && s != inf);
+  const int cmp = x.compare(i);
+  if (cmp > 0) {
+    s = nextafter(s, +inf);
+    assert(x.compare(s) < 0);
+  } else if (cmp < 0) {
+    i = nextafter(i, -inf);
+    assert(x.compare(i) > 0);
+  }
+}
+
+#endif // test_minimal_nextafter
+
+void test_to_interval_boost() {
+
+  using NT = boost::multiprecision::cpp_int;
+  using Quotient = CGAL::Quotient<NT>;
+  using Traits = CGAL::Real_embeddable_traits<Quotient>;
+  using Interval = typename Traits::To_interval;
+
+  NT n, d;
+  Quotient x;
+  double i, s;
+
+  n = NT("-15284404573383541");
+  d = NT("4503599627370496");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: -3.3938195750112902793" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: -3.3938195750112898352" << std::endl;
+  std::cout << std::endl;
+
+  // Results for current tight using cpp_rational.
+  assert(i == -3.3938195750112902793);
+  assert(s == -3.3938195750112898352);
+}
+
+// In assert, we use values from impl2.
+void test_to_interval_tight_rational_1() {
+
+  // In green, we compare to impl2.
+  #define TESTCASE_RAT_10 // pass all three
+  #define TESTCASE_RAT_11 // impl1: fails tightness (sup is larger, s = 9.3488310472396616291)
+  #define TESTCASE_RAT_12 // pass all three
+  #define TESTCASE_RAT_13 // pass all three
+  #define TESTCASE_RAT_14 // pass all three
+  #define TESTCASE_RAT_15 // pass all three
+  #define TESTCASE_RAT_16 // pass all three
+  #define TESTCASE_RAT_17 // pass all three
+  #define TESTCASE_RAT_18 // pass all three
+  #define TESTCASE_RAT_19 // pass all three
+
+  using NT = boost::multiprecision::cpp_int;
+  using Quotient = CGAL::Quotient<NT>;
+  using Traits = CGAL::Real_embeddable_traits<Quotient>;
+  using Interval = typename Traits::To_interval;
+
+  // std::cout << std::endl;
+  // std::cout << boost::typeindex::type_id<Quotient>() << std::endl;
+  // std::cout << std::endl;
+  // std::cout << boost::typeindex::type_id<typename Traits::Type>() << std::endl;
+  // std::cout << std::endl;
+
+  NT n, d;
+  Quotient x;
+  double i, s;
+
+  std::cout << std::endl;
+  std::cout << "- T1 testing tight interval for rationals ..." << std::endl;
+  std::cout << std::endl;
+
+  #ifdef TESTCASE_RAT_10 // small numbers
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE0 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = NT("39792587355159975");
+  d = NT("140737488355328");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 282.7433388230813307" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 282.74333882308138755" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 282.7433388230813307);
+  assert(s == 282.74333882308138755);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_11 // large numbers
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE1 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = NT("772537196160711547532081795586792063331305895970601529435744397743492241616327030886637827664482971614281724796166908515292029740442872965475211471498392497954317530347232852540146110053764627070672243390766540271554856759037331142360111552286202392826786995364211101723592791550906796165626083442695020580821188398298798456115881346136681033873");
+  d = NT("82634630175374856683315372867724319098240552701588533218371381248009342768269285501674184091886435054368116496214846441734481770666205690731018817430937185570378353100803926136323598244976110318516454816403989543192819758059431171537258117598056453283568595627159988837663160716950017789671313834717457946818990093589809113731838629064768225280");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 9.3488310472396563" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 9.3488310472396580764" << std::endl;
+  std::cout << std::endl;
+
+  // Results for current tight using cpp_rational.
+  assert(i == 9.3488310472396563);
+  assert(s == 9.3488310472396580764);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_12
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE2 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = NT("772537196160711547532081795586792063331305895970601529435744397743492241616327030886637827664482971614281724796166908515292029740442872965475211471498392497954317530347232852540146110053764627070672243390766540271554856759037331142360111552286202392826786995364211101723592791550906796165626083442695020580821188398298798456115881346136681033873");
+  d = NT("1");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 1.7976931348623157081e+308" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: inf" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == std::numeric_limits<double>::max());
+  assert(s == std::numeric_limits<double>::infinity());
+
+  #endif
+
+  #ifdef TESTCASE_RAT_13
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE3 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = NT("1");
+  d = NT("772537196160711547532081795586792063331305895970601529435744397743492241616327030886637827664482971614281724796166908515292029740442872965475211471498392497954317530347232852540146110053764627070672243390766540271554856759037331142360111552286202392826786995364211101723592791550906796165626083442695020580821188398298798456115881346136681033873");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.0 or higher" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.0 or higher" << std::endl;
+  std::cout << std::endl;
+
+  assert(i >= 0.0 && i <= std::numeric_limits<double>::min() * 2.0);
+  assert(s >= 0.0 && s <= std::numeric_limits<double>::min() * 2.0);
+  assert(i <= s);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_14
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE4 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = NT("10");
+  d = NT("10");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 1" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 1" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 1.0);
+  assert(s == 1.0);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_15
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE5 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = NT("1");
+  d = NT("6");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.16666666666666665741" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.16666666666666668517" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 0.16666666666666665741);
+  assert(s == 0.16666666666666668517);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_16
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE6 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = +NT("6");
+  d = +NT("3");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: " << 6.0 / 3.0 << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: " << 6.0 / 3.0 << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 2.0);
+  assert(s == 2.0);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_17
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE7 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = -NT("1");
+  d = -NT("2");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: " << 1.0 / 2.0 << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: " << 1.0 / 2.0 << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 1.0 / 2.0);
+  assert(s == 1.0 / 2.0);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_18
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE8 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = -NT("1");
+  d = +NT("3");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: -0.33333333333333337034" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: -0.33333333333333331483" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == -0.33333333333333337034);
+  assert(s == -0.33333333333333331483);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_19 // small numbers, num > 0 and den < 0
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE9 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  n = +NT("39792587355159975");
+  d = -NT("140737488355328");
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: -282.74333882308138755" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: -282.7433388230813307" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == -282.74333882308138755);
+  assert(s == -282.7433388230813307);
+
+  #endif
+
+  std::cout << "---SUCCESS ALL---" << std::endl;
+  std::cout << std::endl;
+}
+
+// In assert, we use values from impl2.
+void test_to_interval_tight_rational_2() {
+
+  // In green, we compare to impl2.
+  #define TESTCASE_RAT_20  // pass all three
+  #define TESTCASE_RAT_21  // impl1: fails tightness (sup is larger, s = 0.43464565325999998668)
+  #define TESTCASE_RAT_22  // pass all three
+  #define TESTCASE_RAT_23  // pass all three
+  #define TESTCASE_RAT_24  // pass all three
+  #define TESTCASE_RAT_25  // pass all three
+  #define TESTCASE_RAT_26  // pass all three
+  #define TESTCASE_RAT_27  // pass all three
+  #define TESTCASE_RAT_28  // pass all three
+  #define TESTCASE_RAT_29  // pass all three
+  #define TESTCASE_RAT_210 // impl1, impl3: fails tightness (sup is larger, s = 5.9425938166208590782e+26)
+  #define TESTCASE_RAT_211 // impl1, impl3: fails tightness (inf is smaller, i = 3602879701896396.5, sup is larger, s = 3602879701896398)
+
+  using NT = boost::multiprecision::cpp_int;
+  using Quotient = CGAL::Quotient<NT>;
+  using Traits = CGAL::Real_embeddable_traits<Quotient>;
+  using Interval = typename Traits::To_interval;
+
+  NT n, d;
+  Quotient x;
+  double i, s;
+
+  std::cout << std::endl;
+  std::cout << "- T2 testing tight interval for rationals ..." << std::endl;
+  std::cout << std::endl;
+
+  #ifdef TESTCASE_RAT_20
+
+  std::cout << "TEST 0" << std::endl; // num, case 2
+
+  n = NT("-15284404573383541");
+  d = NT("4503599627370496");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: -3.3938195750112902793" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: -3.3938195750112898352" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == -3.3938195750112902793);
+  assert(s == -3.3938195750112898352);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_21
+
+  std::cout << "TEST 1" << std::endl; // num, case 4
+
+  const double nn = 0.43464565326;
+  x = Quotient(nn);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.43464565325999998668" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.43464565325999998668" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 0.43464565325999998668);
+  assert(s == 0.43464565325999998668);
+  assert(i == s);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_22
+
+  std::cout << "TEST 2" << std::endl; // num, case 4
+
+  n = NT("1");
+  d = NT("2");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.5" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.5" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 0.5);
+  assert(s == 0.5);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_23
+
+  std::cout << "TEST 3" << std::endl; // shift = 0
+
+  n = NT("7725371961607115");
+  d = NT("1");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 7725371961607115" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 7725371961607115" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 7725371961607115);
+  assert(s == 7725371961607115);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_24
+
+  std::cout << "TEST 4" << std::endl;
+
+  n = NT("772537196160711547532081795586792063331305895970601529435744397743492241616327030886637827664482971614281724796166908515292029740442872965475211471498392497954317530347232852540146110053764627070672243390766540271554856759037331142360111552286202392826786995364211101723592791550906796165626083442695020580821188398298798456115881346136681033873");
+  d = NT("1");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 1.7976931348623157081e+308" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: inf" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == std::numeric_limits<double>::max());
+  assert(s == std::numeric_limits<double>::infinity());
+
+  #endif
+
+  #ifdef TESTCASE_RAT_25
+
+  std::cout << "TEST 5" << std::endl;
+
+  n = NT("1");
+  d = NT("772537196160711547532081795586792063331305895970601529435744397743492241616327030886637827664482971614281724796166908515292029740442872965475211471498392497954317530347232852540146110053764627070672243390766540271554856759037331142360111552286202392826786995364211101723592791550906796165626083442695020580821188398298798456115881346136681033873");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.0 or higher" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.0 or higher" << std::endl;
+  std::cout << std::endl;
+
+  assert(i >= 0.0 && i <= std::numeric_limits<double>::min() * 2.0);
+  assert(s >= 0.0 && s <= std::numeric_limits<double>::min() * 2.0);
+  assert(i <= s);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_26
+
+  std::cout << "TEST 6" << std::endl; // case shift > 0 && p_bits = 51, subcase1
+
+  n = NT("1");
+  d = NT("10");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.099999999999999991673" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.10000000000000000555" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 0.099999999999999991673);
+  assert(s == 0.10000000000000000555);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_27
+
+  std::cout << "TEST 7" << std::endl; // non representable double
+
+  n = NT("1");
+  d = NT("3");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0.33333333333333331483" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0.33333333333333337034" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 0.33333333333333331483);
+  assert(s == 0.33333333333333337034);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_28
+
+  std::cout << "TEST 8" << std::endl; // fails assertion (q_bits == num_dbl_digits || r != 0)
+
+  n = NT("21");
+  d = NT("3");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 7" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 7" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 7);
+  assert(s == 7);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_29
+
+  std::cout << "TEST 9" << std::endl; // case shift > 0 && p_bits = 51, subcase3
+
+  n = NT("17");
+  d = NT("3");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 5.6666666666666660745" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 5.6666666666666669627" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 5.6666666666666660745);
+  assert(s == 5.6666666666666669627);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_210
+
+  std::cout << "TEST 10" << std::endl; // case shift < 0 && p_bits = 51
+
+  n = NT("7725371961607115475320817955");
+  d = NT("13");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 5.9425938166208577038e+26" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 5.942593816620858391e+26" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 5.9425938166208577038e+26);
+  assert(s == 5.942593816620858391e+26);
+
+  #endif
+
+  #ifdef TESTCASE_RAT_211
+
+  std::cout << "TEST 11" << std::endl; // case shift = 0 && p_bits = 51 && cmp = 0, subcase3
+
+  n = NT("36028797018963975");
+  d = NT("10");
+
+  x = Quotient(n, d);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 3602879701896397.5" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 3602879701896397.5" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 3602879701896397.5);
+  assert(s == 3602879701896397.5);
+
+  #endif
+
+  std::cout << "---SUCCESS ALL---" << std::endl;
+  std::cout << std::endl;
+}
+
+void test_to_interval_tight_integer() {
+
+  #define TESTCASE_INT_0
+  #define TESTCASE_INT_1
+  #define TESTCASE_INT_2
+  #define TESTCASE_INT_3
+
+  using NT = boost::multiprecision::cpp_int;
+  using Traits = CGAL::Real_embeddable_traits<NT>;
+  using Interval = typename Traits::To_interval;
+
+  // std::cout << std::endl;
+  // std::cout << boost::typeindex::type_id<NT>() << std::endl;
+  // std::cout << std::endl;
+  // std::cout << boost::typeindex::type_id<typename Traits::Type>() << std::endl;
+  // std::cout << std::endl;
+
+  NT x;
+  double i, s;
+
+  std::cout << std::endl;
+  std::cout << "- T testing tight interval for integers ..." << std::endl;
+  std::cout << std::endl;
+
+  #ifdef TESTCASE_INT_0
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE0 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  x = NT("0");
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 0" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 0" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 0.0);
+  assert(s == 0.0);
+
+  #endif
+
+  #ifdef TESTCASE_INT_1
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE1 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  x = NT("5");
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 5" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: 5" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == 5.0);
+  assert(s == 5.0);
+
+  #endif
+
+  #ifdef TESTCASE_INT_2
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE2 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  x = NT(-12);
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: -12" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: -12" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == -12.0);
+  assert(s == -12.0);
+
+  #endif
+
+  #ifdef TESTCASE_INT_3
+
+  std::cout << "=============" << std::endl;
+  std::cout << "CASE3 RESULT:" << std::endl;
+  std::cout << "=============" << std::endl;
+
+  x = NT("772537196160711547532081795586792063331305895970601529435744397743492241616327030886637827664482971614281724796166908515292029740442872965475211471498392497954317530347232852540146110053764627070672243390766540271554856759037331142360111552286202392826786995364211101723592791550906796165626083442695020580821188398298798456115881346136681033873");
+  std::tie(i, s) = Interval()(x);
+
+  std::cout << std::endl;
+  std::cout << "inf: " << i << std::endl;
+  std::cout << "ref: 1.7976931348623157081e+308" << std::endl;
+  std::cout << "sup: " << s << std::endl;
+  std::cout << "ref: inf" << std::endl;
+  std::cout << std::endl;
+
+  assert(i == std::numeric_limits<double>::max());
+  assert(s == std::numeric_limits<double>::infinity());
+
+  #endif
+
+  std::cout << "---SUCCESS ALL---" << std::endl;
+  std::cout << std::endl;
+}
+
+void test_shift_positive() {
+  {
+    double d = (1LL << 53) - 1;
+    auto shift = std::numeric_limits<double>::max_exponent - (std::numeric_limits<double>::digits);
+    auto r = CGAL::Boost_MP_internal::shift_positive_interval({d,d},shift);
+    d = ldexp(d,shift);
+    assert(!isinf(d));
+    assert(r.first == d && d == r.second);
+  }
+  {
+    double d = (1LL << 52);
+    auto shift = std::numeric_limits<double>::max_exponent - std::numeric_limits<double>::digits + 1;
+    auto r = CGAL::Boost_MP_internal::shift_positive_interval({d,d},shift);
+    d = ldexp(d,shift);
+    assert(d >= (std::numeric_limits<double>::max)());
+    assert(r.first <= d && d <= r.second);
+  }
+  {
+    double d = (1LL << 53) - 1;
+    auto shift = std::numeric_limits<double>::min_exponent - std::numeric_limits<double>::digits - 1;
+    auto r = CGAL::Boost_MP_internal::shift_positive_interval({d,d},shift);
+    d = ldexp(d,shift);
+    assert(d <= (std::numeric_limits<double>::min)());
+    assert(r.first <= d && d <= r.second);
+  }
+  {
+    double d = (1LL << 53) - 2;
+    auto shift = std::numeric_limits<double>::min_exponent - std::numeric_limits<double>::digits - 1;
+    auto r = CGAL::Boost_MP_internal::shift_positive_interval({d,d},shift);
+    d = ldexp(d,shift);
+    assert(d < (std::numeric_limits<double>::min)());
+    assert(r.first <= d && d <= r.second);
+  }
+  {
+    double d = (1LL << 52);
+    auto shift = std::numeric_limits<double>::min_exponent - std::numeric_limits<double>::digits;
+    auto r = CGAL::Boost_MP_internal::shift_positive_interval({d,d},shift);
+    d = ldexp(d,shift);
+    assert(d == (std::numeric_limits<double>::min)());
+    assert(r.first == d && d == r.second);
+  }
+}
+#endif // CGAL_USE_BOOST_MP
+
+int main() {
+#ifdef CGAL_USE_BOOST_MP
+  test_shift_positive();
+#endif
+
+  // Make sure we have the same seed.
+  CGAL::get_default_random() = CGAL::Random(0);
+  std::cout.precision(20);
+  #ifdef CGAL_USE_BOOST_MP
+
+  #if false
+  test_boost_eval_lehmer();
+  test_minimal_nextafter();
+  #endif
+  test_minimal_boost_gcd();
+  test_to_interval_boost();
+
+  // Test rational types.
+  test_to_interval_tight_rational_1();
+  test_to_interval_tight_rational_2();
+
+  // Test integer types.
+  test_to_interval_tight_integer();
+
+  #endif // CGAL_USE_BOOST_MP
+  return EXIT_SUCCESS;
+}