Merge pull request #2043 from ERGO-Code/fix-2041

Treating `|x|<1` as a special case in `floor(HighsCDouble x)` and `ceil(HighsCDouble x)`

check/CMakeLists.txt

@@ -257,9 +257,10 @@ if ((NOT FAST_BUILD OR ALL_TESTS) AND NOT (BUILD_EXTRA_UNIT_ONLY))
         "standgub\; 1.25769949\;"
         )
   elseif(WIN32)
-    # on windows e226 model status is unknown, rel gap e00
+    # on windows e226 model status is unknown
+    # on windows 25fv47 model status can be unknown, with objective 5.5018458957e+03
     set(pdlpInstances
-        "25fv47\; 5.50184588\;"
+        "25fv47\; 5.5018458\;"
         "adlittle\; 2.25494963\;"
         "afiro\;-4.64753142\;"
         "avgas\;-7.749999999\;"
@@ -455,4 +456,4 @@ if (BUILD_EXTRA_UNIT_TESTS AND BUILD_EXTRA_UNIT_ONLY)
                         DEPENDS unit-test-extra-build)
     set_tests_properties(unit_tests_extra PROPERTIES TIMEOUT 10000) 
 
-endif()
+endif()

check/TestCheckSolution.cpp

@@ -1,4 +1,4 @@
-//#include <cstdio>
+// #include <cstdio>
 #include <iostream>
 
 #include "HCheckConfig.h"

check/TestHighsIntegers.cpp

@@ -1,6 +1,8 @@
 #include "HCheckConfig.h"
 #include "catch.hpp"
+#include "util/HighsCDouble.h"
 #include "util/HighsIntegers.h"
+#include "util/HighsRandom.h"
 
 const bool dev_run = false;
 
@@ -38,3 +40,107 @@ TEST_CASE("HighsIntegers", "[util]") {
 
   if (dev_run) printf("integral scalar is %g\n", integralscalar);
 }
+
+void testCeil(HighsCDouble x) {
+  double ceil_x;
+  double double_x;
+  ceil_x = double(ceil(x));
+  double_x = double(x);
+  REQUIRE(ceil_x >= double_x);
+  REQUIRE(ceil(x) >= x);
+}
+
+void testFloor(HighsCDouble x) {
+  double floor_x;
+  double double_x;
+  floor_x = double(floor(x));
+  double_x = double(x);
+  REQUIRE(floor_x <= double_x);
+  REQUIRE(floor(x) <= x);
+}
+TEST_CASE("HighsCDouble-ceil", "[util]") {
+  HighsCDouble x;
+  x = -1e-34;
+  testCeil(x);
+  x = -1e-32;
+  testCeil(x);
+  x = -1e-30;
+  testCeil(x);
+  x = -1e-23;
+  testCeil(x);
+  x = -1e-12;
+  testCeil(x);
+  x = -1e-1;
+  testCeil(x);
+  x = -0.99;
+  testCeil(x);
+
+  x = 0.99;
+  testCeil(x);
+  x = 1e-1;
+  testCeil(x);
+  x = 1e-12;
+  testCeil(x);
+  // This and rest failed in #2041
+  x = 1e-23;
+  testCeil(x);
+  x = 1e-30;
+  testCeil(x);
+  x = 1e-32;
+  testCeil(x);
+  x = 1e-34;
+  testCeil(x);
+
+  HighsRandom rand;
+  for (HighsInt k = 0; k < 1000; k++) {
+    double man = rand.fraction();
+    HighsInt power = 2 - rand.integer(5);
+    double exp = std::pow(10, power);
+    x = man * exp;
+    testCeil(x);
+  }
+}
+
+TEST_CASE("HighsCDouble-floor", "[util]") {
+  HighsCDouble x;
+
+  x = 1e-34;
+  testFloor(x);
+  x = 1e-32;
+  testFloor(x);
+  x = 1e-30;
+  testFloor(x);
+  x = 1e-23;
+  testFloor(x);
+  x = 1e-12;
+  testFloor(x);
+  x = 1e-1;
+  testFloor(x);
+  x = 0.99;
+  testFloor(x);
+
+  x = -0.99;
+  testFloor(x);
+  x = -1e-1;
+  testFloor(x);
+  x = -1e-12;
+  testFloor(x);
+  // This and rest failed in #2041
+  x = -1e-23;
+  testFloor(x);
+  x = -1e-30;
+  testFloor(x);
+  x = -1e-32;
+  testFloor(x);
+  x = -1e-34;
+  testFloor(x);
+
+  HighsRandom rand;
+  for (HighsInt k = 0; k < 1000; k++) {
+    double man = rand.fraction();
+    HighsInt power = 2 - rand.integer(5);
+    double exp = std::pow(10, power);
+    x = man * exp;
+    testFloor(x);
+  }
+}

check/TestMipSolver.cpp

@@ -679,16 +679,16 @@ TEST_CASE("IP-infeasible-unbounded", "[highs_test_mip_solver]") {
       highs.run();
       HighsModelStatus required_model_status;
       if (k == 0) {
-	// Presolve off
+        // Presolve off
         if (l == 0) {
-	  // MIP solver proves infeasiblilty
+          // MIP solver proves infeasiblilty
           required_model_status = HighsModelStatus::kInfeasible;
         } else {
-	  // Relaxation is unbounded, but origin is feasible
+          // Relaxation is unbounded, but origin is feasible
           required_model_status = HighsModelStatus::kUnbounded;
         }
       } else {
-	// Presolve on, and identifies primal infeasible or unbounded
+        // Presolve on, and identifies primal infeasible or unbounded
         required_model_status = HighsModelStatus::kUnboundedOrInfeasible;
       }
       REQUIRE(highs.getModelStatus() == required_model_status);

src/util/HighsCDouble.h

@@ -289,6 +289,12 @@ class HighsCDouble {
   }
 
   friend HighsCDouble floor(const HighsCDouble& x) {
+    // Treat |x| < 1 as special case, as per (for example)
+    // https://github.com/shibatch/tlfloat: see #2041
+    if (abs(x) < 1) {
+      if (x == 0 || x > 0) return HighsCDouble(0.0);
+      return HighsCDouble(-1.0);
+    }
     double floor_x = std::floor(double(x));
     HighsCDouble res;
 
@@ -297,6 +303,12 @@ class HighsCDouble {
   }
 
   friend HighsCDouble ceil(const HighsCDouble& x) {
+    // Treat |x| < 1 as special case, as per (for example)
+    // https://github.com/shibatch/tlfloat: see #2041
+    if (abs(x) < 1) {
+      if (x == 0 || x < 0) return HighsCDouble(0.0);
+      return HighsCDouble(1.0);
+    }
     double ceil_x = std::ceil(double(x));
     HighsCDouble res;