Merge pull request #5510 from pentacular/wasm

Round-to-nearest support for intervals

Number_types/include/CGAL/FPU.h

@@ -114,10 +114,11 @@ extern "C" {
 
 // Only define CGAL_USE_SSE2 for 64 bits where malloc has a suitable
 // alignment, 32 bits is too dangerous.
-#if defined(CGAL_HAS_SSE2) && (defined(__x86_64__) || defined(_M_X64))
+#if defined CGAL_HAS_SSE2 && \
+  (defined __x86_64__ || defined _M_X64) && \
+  !defined CGAL_ALWAYS_ROUND_TO_NEAREST
 #  define CGAL_USE_SSE2 1
 #endif
-
 #ifdef CGAL_CFG_DENORMALS_COMPILE_BUG
 double& get_static_minimin(); // Defined in Interval_arithmetic_impl.h
 #endif
@@ -347,18 +348,27 @@ 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 ?
-#define CGAL_IA_ADD(a,b) CGAL_IA_FORCE_TO_DOUBLE((a)+CGAL_IA_STOP_CPROP(b))
-#define CGAL_IA_SUB(a,b) CGAL_IA_FORCE_TO_DOUBLE(CGAL_IA_STOP_CPROP(a)-(b))
-#define CGAL_IA_MUL(a,b) CGAL_IA_FORCE_TO_DOUBLE(CGAL_IA_STOP_CPROP(a)*CGAL_IA_STOP_CPROP(b))
-#define CGAL_IA_DIV(a,b) CGAL_IA_FORCE_TO_DOUBLE(CGAL_IA_STOP_CPROP(a)/CGAL_IA_STOP_CPROP(b))
+inline double IA_up(double d)
+{
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  // In round-to-nearest mode we find the successor instead.
+  // This preserves the interval invariants, but is more
+  // expensive and conservative.
+  return nextafter(d, std::numeric_limits<double>::infinity());
+#else
+  // In round-upward mode we can rely on the hardware
+  // to do the job.
+  return CGAL_IA_FORCE_TO_DOUBLE(d);
+#endif
+}
+#define CGAL_IA_ADD(a,b) IA_up((a)+CGAL_IA_STOP_CPROP(b))
+#define CGAL_IA_SUB(a,b) IA_up(CGAL_IA_STOP_CPROP(a)-(b))
+#define CGAL_IA_MUL(a,b) IA_up(CGAL_IA_STOP_CPROP(a)*CGAL_IA_STOP_CPROP(b))
+#define CGAL_IA_DIV(a,b) IA_up(CGAL_IA_STOP_CPROP(a)/CGAL_IA_STOP_CPROP(b))
 inline double CGAL_IA_SQUARE(double a){
   double b = CGAL_IA_STOP_CPROP(a); // only once
-  return CGAL_IA_FORCE_TO_DOUBLE(b*b);
+  return IA_up(b*b);
 }
-#define CGAL_IA_SQRT(a) \
-        CGAL_IA_FORCE_TO_DOUBLE(CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(a)))
-
-
 #if defined CGAL_SAFE_SSE2
 
 #define CGAL_IA_SETFPCW(CW) _MM_SET_ROUNDING_MODE(CW)
@@ -476,9 +486,15 @@ inline
 FPU_CW_t
 FPU_get_cw (void)
 {
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+    CGAL_assertion_code(FPU_CW_t cw; CGAL_IA_GETFPCW(cw);)
+    CGAL_assertion(cw == CGAL_FE_TONEAREST);
+    return CGAL_FE_TONEAREST;
+#else
     FPU_CW_t cw;
     CGAL_IA_GETFPCW(cw);
     return cw;
+#endif
 }
 
 // User interface (cont):
@@ -487,28 +503,43 @@ inline
 void
 FPU_set_cw (FPU_CW_t cw)
 {
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  CGAL_assertion(cw == CGAL_FE_TONEAREST);
+#else
   CGAL_IA_SETFPCW(cw);
+#endif
 }
 
 inline
 FPU_CW_t
 FPU_get_and_set_cw (FPU_CW_t cw)
 {
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+    CGAL_assertion(cw == CGAL_FE_TONEAREST);
+    return CGAL_FE_TONEAREST;
+#else
     FPU_CW_t old = FPU_get_cw();
     FPU_set_cw(cw);
     return old;
+#endif
 }
 
 
 // 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(FPU_CW_t r = CGAL_FE_UPWARD)
+  Protect_FPU_rounding(FPU_CW_t r = CGAL_FE_PROTECTED)
     : backup( FPU_get_and_set_cw(r) ) {}
 
   ~Protect_FPU_rounding()
@@ -524,7 +555,7 @@ template <>
 struct Protect_FPU_rounding<false>
 {
   Protect_FPU_rounding() {}
-  Protect_FPU_rounding(FPU_CW_t /*= CGAL_FE_UPWARD*/) {}
+  Protect_FPU_rounding(FPU_CW_t /*= CGAL_FE_PROTECTED */) {}
 };
 
 
@@ -538,13 +569,13 @@ struct Checked_protect_FPU_rounding
 {
   Checked_protect_FPU_rounding()
   {
-    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_UPWARD);
+    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_PROTECTED);
   }
 
   Checked_protect_FPU_rounding(FPU_CW_t r)
     : Protect_FPU_rounding<Protected>(r)
   {
-    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_UPWARD);
+    CGAL_expensive_assertion(FPU_get_cw() == CGAL_FE_PROTECTED);
   }
 };
 
@@ -555,6 +586,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
@@ -579,6 +612,21 @@ inline void force_ieee_double_precision()
 #endif
 }
 
+inline double IA_sqrt_up(double a) {
+  return IA_up(CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(a)));
+}
+
+inline double IA_sqrt_toward_zero(double d) {
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  return (d > 0.0) ? nextafter(std::sqrt(d), 0.) : 0.0;
+#else
+  FPU_set_cw(CGAL_FE_DOWNWARD);
+  double i = (d > 0.0) ? CGAL_IA_FORCE_TO_DOUBLE(CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(d))) : 0.0;
+  FPU_set_cw(CGAL_FE_UPWARD);
+  return i;
+#endif
+}
+
 } //namespace CGAL
 
 #ifdef CGAL_HEADER_ONLY

Number_types/include/CGAL/Interval_nt.h

@@ -1141,17 +1141,15 @@ namespace INTERN_INTERVAL_NT {
       // it helps significantly, it might even hurt by introducing a
       // dependency.
     }
-#else
+#else // no __AVX512F__
     // TODO: Alternative for computing CGAL_IA_SQRT_DOWN(d.inf()) exactly
     // without changing the rounding mode:
     // - compute x = CGAL_IA_SQRT(d.inf())
     // - compute y = CGAL_IA_SQUARE(x)
     // - if y==d.inf() use x, else use -CGAL_IA_SUB(CGAL_IA_MIN_DOUBLE,x)
-    FPU_set_cw(CGAL_FE_DOWNWARD);
-    double i = (d.inf() > 0.0) ? CGAL_IA_SQRT(d.inf()) : 0.0;
-    FPU_set_cw(CGAL_FE_UPWARD);
-#endif
-    return Interval_nt<Protected>(i, CGAL_IA_SQRT(d.sup()));
+    double i = IA_sqrt_toward_zero(d.inf());
+#endif // no __AVX512F__
+    return Interval_nt<Protected>(i, IA_sqrt_up(d.sup()));
   }
 
   template <bool Protected>

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")

Number_types/test/Number_types/Interval_nt.cpp

@@ -31,12 +31,18 @@ bool spiral_test()
       break;
   };
 
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+  return i == 365;
+#else
   return i == 396;
+#endif
 }
 
 // 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()
 {
@@ -92,12 +98,31 @@ 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;
+  ) // DEBUG
+  if (i != 54) {
+    return false;
+  }
+  // When we round to nearest it doesn't quite converge.
+  if (a.sup() > 3/(double(1<<30)*(1<<22))) {
+    return false;
+  }
+  if (-3/(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
 }
 
 
@@ -164,9 +189,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(-CGAL_IA_MIN_DOUBLE, CGAL_IA_MIN_DOUBLE))
+      && b.is_same(IA_nt::smallest())
+      && c.is_same(IA_nt(0, CGAL_IA_MIN_DOUBLE));
+#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 +463,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;
 }

Number_types/test/Number_types/Interval_nt_nearest.cpp

@@ -0,0 +1,4 @@
+#define CGAL_ALWAYS_ROUND_TO_NEAREST
+#define CGAL_DEBUG
+#define CGAL_CHECK_EXPENSIVE
+#include "Interval_nt.cpp"

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")

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"