Avoid overflow in float interval medians

changelog.in

@@ -135,6 +135,17 @@ Issue: 83
 Build CMake examples only when their required variable modules and MPFR support
 are enabled.
 
+[ENTRY]
+Module: float
+What: bug
+Rank: major
+Issue: 164
+Thanks: Pierre Talbot
+[DESCRIPTION]
+Compute interval medians without overflowing at large positive or negative
+bounds. Float search over an unbounded FlatZinc variable can otherwise fail to
+make progress.
+
 [ENTRY]
 Module: flatzinc
 What:   new

gecode/float/val.hpp

@@ -80,7 +80,24 @@ namespace Gecode {
   }
   forceinline FloatNum
   FloatVal::med(void) const {
-    return gecode_boost::numeric::median(x);
+    const FloatNum l = x.lower();
+    const FloatNum u = x.upper();
+    const FloatNum inf = std::numeric_limits<FloatNum>::infinity();
+    const FloatNum max = std::numeric_limits<FloatNum>::max();
+    if (l == -inf) {
+      if (u == -inf)
+        return l;
+      // Preserve the extended endpoint for semi-infinite intervals.
+      return (u == inf) ? 0.0 : l;
+    }
+    if (u == inf)
+      return u;
+
+    Float::Rounding r;
+    const FloatNum h = max / 2.0;
+    if ((u > h) || (l < -h))
+      return std::ldexp(r.median(l / 2.0,u / 2.0),1);
+    return r.median(l,u);
   }
 
   forceinline bool

test/float/basic.cpp

@@ -35,6 +35,8 @@
 
 #include "test/float.hh"
 
+#include <cmath>
+
 namespace Test { namespace Float {
 
    /// %Tests for basic setup
@@ -73,6 +75,47 @@ namespace Test { namespace Float {
        Basic b1(3,1.5);
        Basic b2(Gecode::FloatVal(-2,10),1.5);
      }
+
+     /// Test interval medians at the numerical limits and special values
+     class Median : public Base {
+     private:
+       /// Test median inclusion and splitting of non-tight intervals
+       static bool valid(const Gecode::FloatVal& x) {
+         const Gecode::FloatNum m = x.med();
+         return (m >= x.min()) && (m <= x.max()) &&
+           (x.tight() || ((m > x.min()) && (m < x.max())));
+       }
+     public:
+       /// Create and register test
+       Median(void) : Base("Float::Basic::Median") {}
+       /// Perform test
+       virtual bool run(void) {
+         using namespace Gecode;
+         const FloatNum m = Gecode::Float::Limits::max;
+         const FloatNum i = std::numeric_limits<FloatNum>::infinity();
+         const FloatNum d = std::numeric_limits<FloatNum>::denorm_min();
+         const FloatVal negative(-m, -m / 2.0);
+         const FloatVal positive(m / 2.0, m);
+         const FloatVal mixed(-m, m);
+         // The subtraction-based formula rounds this median one ULP upward.
+         const FloatVal closest(3.07187504409887e-53,
+                                7.805365587017815e-45);
+         return valid(negative) && valid(positive) && valid(mixed) &&
+           (closest.med() == 3.902682808868283e-45) &&
+           (FloatVal(-i,-1.0).med() == -i) &&
+           (FloatVal(1.0,i).med() == i) &&
+           (FloatVal(-i,i).med() == 0.0) &&
+           (FloatVal(-i,-i).med() == -i) &&
+           (FloatVal(i,i).med() == i) &&
+           !std::signbit(FloatVal(-0.0,0.0).med()) &&
+           (FloatVal(0.0,2.0*d).med() == d) &&
+           (FloatVal(-2.0*d,0.0).med() == -d);
+       }
+     };
+
+     namespace {
+       Median median;
+     }
      //@}
 
    }