diff --git a/Number_types/include/CGAL/FPU.h b/Number_types/include/CGAL/FPU.h
index e4e45430f98..0d7bf61e60c 100644
--- a/Number_types/include/CGAL/FPU.h
+++ b/Number_types/include/CGAL/FPU.h
@@ -1,8 +1,3 @@
-// TODO: Remove
-#ifndef CGAL_ALWAYS_ROUND_TO_NEAREST
-#define CGAL_ALWAYS_ROUND_TO_NEAREST
-#endif
-
 // Copyright (c) 1998-2019
 // Utrecht University (The Netherlands),
 // ETH Zurich (Switzerland),
@@ -122,7 +117,6 @@ extern "C" {
 #if defined(CGAL_HAS_SSE2) && (defined(__x86_64__) || defined(_M_X64))
 #  define CGAL_USE_SSE2 1
 #endif
-#undef CGAL_USE_SSE2
 #ifdef CGAL_CFG_DENORMALS_COMPILE_BUG
 double& get_static_minimin(); // Defined in Interval_arithmetic_impl.h
 #endif
@@ -352,7 +346,17 @@ inline double IA_bug_sqrt(double d)
 // With GCC, we can do slightly better : test with __builtin_constant_p()
 // that both arguments are constant before stopping one of them.
 // Use inline functions instead ?
-inline double CGAL_IA_UP(double d){return nextafter(d,std::numeric_limits<double>::infinity());}
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+inline double CGAL_IA_UP(double d)
+{
+  return nextafter(d,std::numeric_limits<double>::infinity());
+}
+#else
+inline double CGAL_IA_UP(double d)
+{
+  return d;
+}
+#endif
 #define CGAL_IA_ADD(a,b) CGAL_IA_UP((a)+CGAL_IA_STOP_CPROP(b))
 #define CGAL_IA_SUB(a,b) CGAL_IA_UP(CGAL_IA_STOP_CPROP(a)-(b))
 #define CGAL_IA_MUL(a,b) CGAL_IA_UP(CGAL_IA_STOP_CPROP(a)*CGAL_IA_STOP_CPROP(b))
@@ -496,29 +500,15 @@ inline
 void
 FPU_set_cw (FPU_CW_t cw)
 {
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-  assert(cw == CGAL_FE_TONEAREST);
-  assert(FPU_get_cw() == CGAL_FE_TONEAREST);
-#endif
   CGAL_IA_SETFPCW(cw);
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-  assert(FPU_get_cw() == CGAL_FE_TONEAREST);
-#endif
 }
 
 inline
 FPU_CW_t
 FPU_get_and_set_cw (FPU_CW_t cw)
 {
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-  assert(cw == CGAL_FE_TONEAREST);
-  assert(FPU_get_cw() == CGAL_FE_TONEAREST);
-#endif
     FPU_CW_t old = FPU_get_cw();
     FPU_set_cw(cw);
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-  assert(FPU_get_cw() == CGAL_FE_TONEAREST);
-#endif
     return old;
 }
 
@@ -526,18 +516,18 @@ FPU_get_and_set_cw (FPU_CW_t cw)
 // A class whose constructor sets the FPU mode to +inf, saves a backup of it,
 // and whose destructor resets it back to the saved state.
 
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+#define CGAL_FE_PROTECTED CGAL_FE_TONEAREST
+#else
+#define CGAL_FE_PROTECTED CGAL_FE_UPWARD
+#endif
+
 template <bool Protected = true> struct Protect_FPU_rounding;
 
 template <>
 struct Protect_FPU_rounding<true>
 {
-  Protect_FPU_rounding(
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-    FPU_CW_t r = CGAL_FE_TONEAREST
-#else
-    FPU_CW_t r = CGAL_FE_UPWARD
-#endif
-  )
+  Protect_FPU_rounding(FPU_CW_t r = CGAL_FE_PROTECTED)
     : backup( FPU_get_and_set_cw(r) ) {}
 
   ~Protect_FPU_rounding()
@@ -553,11 +543,7 @@ template <>
 struct Protect_FPU_rounding<false>
 {
   Protect_FPU_rounding() {}
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-  Protect_FPU_rounding(FPU_CW_t /*= CGAL_FE_TONEAREST*/) {}
-#else
-  Protect_FPU_rounding(FPU_CW_t /*= CGAL_FE_UPWARD*/) {}
-#endif
+  Protect_FPU_rounding(FPU_CW_t /*= CGAL_FE_PROTECTED */) {}
 };
 
 
@@ -571,21 +557,13 @@ struct Checked_protect_FPU_rounding
 {
   Checked_protect_FPU_rounding()
   {
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_TONEAREST);
-#else
-    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_UPWARD);
-#endif
+    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_PROTECTED);
   }
 
   Checked_protect_FPU_rounding(FPU_CW_t r)
     : Protect_FPU_rounding<Protected>(r)
   {
-#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
-    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_TONEAREST);
-#else
-    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_UPWARD);
-#endif
+    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_PROTECTED);
   }
 };
 
@@ -596,6 +574,8 @@ struct Checked_protect_FPU_rounding
 // Its destructor restores the FPU state as it was previously.
 // Note that this affects "long double" as well, and other potential side effects.
 // And note that it does not (cannot) "fix" the same problem for the exponent.
+//
+// (How should this interact with ALWAYS_ROUND_TO_NEAREST?)
 
 struct Set_ieee_double_precision
 #ifdef CGAL_FPU_HAS_EXCESS_PRECISION
diff --git a/Number_types/test/Number_types/CMakeLists.txt b/Number_types/test/Number_types/CMakeLists.txt
index fdc0448365e..aa51590da93 100644
--- a/Number_types/test/Number_types/CMakeLists.txt
+++ b/Number_types/test/Number_types/CMakeLists.txt
@@ -31,6 +31,7 @@ create_single_source_cgal_program("Gmpz.cpp")
 create_single_source_cgal_program("Gmpzf_new.cpp")
 create_single_source_cgal_program("int.cpp")
 create_single_source_cgal_program("Interval_nt.cpp")
+create_single_source_cgal_program("Interval_nt_nearest.cpp")
 create_single_source_cgal_program("Interval_nt_new.cpp")
 create_single_source_cgal_program("ioformat.cpp")
 create_single_source_cgal_program("known_bit_size_integers.cpp")
diff --git a/Number_types/test/Number_types/Interval_nt.cpp b/Number_types/test/Number_types/Interval_nt.cpp
index d0a470c1320..4a56129d31f 100644
--- a/Number_types/test/Number_types/Interval_nt.cpp
+++ b/Number_types/test/Number_types/Interval_nt.cpp
@@ -36,41 +36,71 @@ bool spiral_test()
 
 // Tests for constant propagation through intervals.
 // This must not be performed otherwise rounding modes are ignored.
-// Non-inlined operators usually stop cprop (*, /, sqrt).
+// On the other hand, if we always round to nearest, then constant propagation
+// is desirable.
+// Note: Non-inlined operators usually stop cprop (*, /, sqrt).
 template < typename IA_nt >
 bool cprop_test()
 {
   // Testing cprop through +.
   IA_nt add = IA_nt(0.00001)+10.1;
   bool good_add = !add.is_point();
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  if (good_add)
+    std::cerr << "ERROR : No constant propagation through operator+." <<std::endl;
+#else
   if (!good_add)
     std::cerr << "ERROR : Constant propagation through operator+." <<std::endl;
+#endif
 
   // Testing cprop through -.
   IA_nt sub = IA_nt(0.00001)-10.1;
   bool good_sub = !sub.is_point();
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  if (good_sub)
+    std::cerr << "ERROR : No constant propagation through operator-." <<std::endl;
+#else
   if (!good_sub)
     std::cerr << "ERROR : Constant propagation through operator-." <<std::endl;
+#endif
 
   // Testing cprop through *.
   IA_nt mul = IA_nt(0.00001)*10.1;
   bool good_mul = !mul.is_point();
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  if (good_mul)
+    std::cerr << "ERROR : No constant propagation through operator*." <<std::endl;
+#else
   if (!good_mul)
     std::cerr << "ERROR : Constant propagation through operator*." <<std::endl;
+#endif
 
   // Testing cprop through /.
   IA_nt div = IA_nt(0.00001)/10.1;
   bool good_div = !div.is_point();
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  if (good_div)
+    std::cerr << "ERROR : No constant propagation through operator/." <<std::endl;
+#else
   if (!good_div)
     std::cerr << "ERROR : Constant propagation through operator/." <<std::endl;
+#endif
 
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  // We have no expectation of propagation through sqrt.
+#else
   // Testing cprop through sqrt.
   IA_nt sqrt2 = CGAL_NTS sqrt(IA_nt(2));
   bool good_sqrt = !sqrt2.is_point();
   if (!good_sqrt)
     std::cerr << "ERROR : Constant propagation through sqrt()." <<std::endl;
+#endif
 
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  return !(good_add || good_sub || good_mul || good_div);
+#else
   return good_add && good_sub && good_mul && good_div && good_sqrt;
+#endif
 }
 
 // Here we iteratively compute sqrt(interval), where interval is [0.5;1.5]
@@ -92,12 +122,32 @@ bool square_root_test()
     a = b;
   };
   a -= 1.0;
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  DEBUG (
+  std::cout << "i          = " << i << std::endl;
+  std::cout << "sup        : " << a.sup() << std::endl;
+  std::cout << "inf        : " << a.inf() << std::endl;
+  std::cout << "width ok ? : " << (-a.inf() == 1/(double(1<<30)*(1<<22))) << std::endl;
+  ) // DEBUG
+  if (i != 54) {
+    return false;
+  }
+  // When we round to nearest it doesn't quite converge.
+  if (a.sup() > 2/(double(1<<30)*(1<<22))) {
+    return false;
+  }
+  if (-2/(double(1<<30)*(1<<22)) > a.inf()) {
+    return false;
+  }
+  return true;
+#else
   DEBUG (
   std::cout << "i          = " << i << std::endl;
   std::cout << "sup = -inf : " << (a.sup() == -a.inf()) << std::endl;
   std::cout << "width ok ? : " << (-a.inf() == 1/(double(1<<30)*(1<<22))) << std::endl;
   ) // DEBUG
   return i==54 && a.sup() == - a.inf() && a.sup() == 1/(double(1<<30)*(1<<22));
+#endif
 }
 
 
@@ -139,6 +189,15 @@ bool overflow_test()
     DEBUG( std::cout << "f = " << f << std::endl; )
   }
 
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  return a.is_same(IA_nt(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())) &&
+         b.is_same(IA_nt(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())) &&
+         c.is_same(IA_nt(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())) &&
+         d.is_same(IA_nt(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())) &&
+         e.is_same(IA_nt(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())) &&
+         f.is_same(IA_nt(std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity())) &&
+         g.is_same(-f);
+#else
   return a.is_same(IA_nt(CGAL_IA_MAX_DOUBLE, std::numeric_limits<double>::infinity())) &&
          b.is_same(IA_nt(CGAL_IA_MAX_DOUBLE, std::numeric_limits<double>::infinity())) &&
          c.is_same(IA_nt::largest()) &&
@@ -146,6 +205,7 @@ bool overflow_test()
          e.is_same(IA_nt::largest()) &&
          f.is_same(IA_nt(CGAL_IA_MAX_DOUBLE, std::numeric_limits<double>::infinity())) &&
          g.is_same(-f);
+#endif
 }
 
 
@@ -164,9 +224,15 @@ bool underflow_test()
   for (i=0; i<20; i++) b = b * b;
   for (i=0; i<20; i++) c = CGAL_NTS square(c);
 
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  return a.is_same(IA_nt(0, 0))
+      && b.is_same(IA_nt(0, 0))
+      && c.is_same(IA_nt(0, 0));
+#else
   return a.is_same(IA_nt(0, CGAL_IA_MIN_DOUBLE))
       && b.is_same(IA_nt::smallest())
       && c.is_same(IA_nt(0, CGAL_IA_MIN_DOUBLE));
+#endif
 }
 
 
@@ -432,10 +498,17 @@ bool test ()
 
 int main()
 {
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  std::cout << "Stress-testing the class Interval_nt<> always rounding to nearest.\n";
+  bool ok = test<CGAL::Interval_nt<> >();
+  std::cout << "\nStress-testing the class Interval_nt_advanced always rounding to nearest.\n";
+  ok &= test<CGAL::Interval_nt_advanced>();
+#else
   std::cout << "Stress-testing the class Interval_nt<>.\n";
   bool ok = test<CGAL::Interval_nt<> >();
   std::cout << "\nStress-testing the class Interval_nt_advanced.\n";
   ok &= test<CGAL::Interval_nt_advanced>();
+#endif
 
   return !ok;
 }
diff --git a/Number_types/test/Number_types/Interval_nt_nearest.cpp b/Number_types/test/Number_types/Interval_nt_nearest.cpp
new file mode 100644
index 00000000000..fe5ce657c6e
--- /dev/null
+++ b/Number_types/test/Number_types/Interval_nt_nearest.cpp
@@ -0,0 +1,2 @@
+#define CGAL_ALWAYS_ROUND_TO_NEAREST
+#include "Interval_nt.cpp"
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
index 08a65b9ffc4..103d6eedb72 100644
--- a/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/CMakeLists.txt
@@ -77,6 +77,7 @@ create_single_source_cgal_program("triangulate_faces_hole_filling_all_search_tes
 create_single_source_cgal_program("test_pmp_remove_border_edge.cpp")
 create_single_source_cgal_program("test_pmp_distance.cpp")
 create_single_source_cgal_program("test_corefinement_and_constraints.cpp")
+create_single_source_cgal_program("test_corefinement_and_constraints_nearest.cpp")
 create_single_source_cgal_program("test_corefinement_bool_op.cpp")
 create_single_source_cgal_program("test_corefine.cpp")
 create_single_source_cgal_program("test_coref_epic_points_identity.cpp")
diff --git a/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_and_constraints_nearest.cpp b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_and_constraints_nearest.cpp
new file mode 100644
index 00000000000..b3c56b5c542
--- /dev/null
+++ b/Polygon_mesh_processing/test/Polygon_mesh_processing/test_corefinement_and_constraints_nearest.cpp
@@ -0,0 +1,2 @@
+#define CGAL_ALWAYS_ROUND_TO_NEAREST
+#include "test_corefinement_and_constraints.cpp"