[CP-SAT] add diophantine equation solver; improve CumulativeTimeTable…

… cuts

ortools/sat/BUILD.bazel

@@ -383,6 +383,7 @@ cc_library(
         ":cp_model_symmetries",
         ":cp_model_utils",
         ":diffn_util",
+        ":diophantine",
         ":inclusion",
         ":presolve_context",
         ":presolve_util",
@@ -1547,3 +1548,14 @@ cc_library(
         "@com_google_absl//absl/strings",
     ],
 )
+
+cc_library(
+    name = "diophantine",
+    srcs = ["diophantine.cc"],
+    hdrs = ["diophantine.h"],
+    deps = [
+        ":util",
+        "@com_google_absl//absl/numeric:int128",
+        "@com_google_absl//absl/types:span",
+    ],
+)

ortools/sat/cp_model.cc

@@ -16,6 +16,7 @@
 #include <cstdint>
 #include <initializer_list>
 #include <limits>
+#include <ostream>
 #include <string>
 #include <vector>
 

ortools/sat/cp_model_checker.cc

@@ -174,49 +174,6 @@ std::string ValidateIntervalsUsedInConstraint(bool after_presolve,
   return "";
 }
 
-template <class LinearExpressionProto>
-bool PossibleIntegerOverflow(const CpModelProto& model,
-                             const LinearExpressionProto& proto,
-                             int64_t offset = 0) {
-  if (offset == std::numeric_limits<int64_t>::min()) return true;
-  int64_t sum_min = -std::abs(offset);
-  int64_t sum_max = +std::abs(offset);
-  for (int i = 0; i < proto.vars_size(); ++i) {
-    const int ref = proto.vars(i);
-    const auto& var_proto = model.variables(PositiveRef(ref));
-    const int64_t min_domain = var_proto.domain(0);
-    const int64_t max_domain = var_proto.domain(var_proto.domain_size() - 1);
-    if (proto.coeffs(i) == std::numeric_limits<int64_t>::min()) return true;
-    const int64_t coeff =
-        RefIsPositive(ref) ? proto.coeffs(i) : -proto.coeffs(i);
-    const int64_t prod1 = CapProd(min_domain, coeff);
-    const int64_t prod2 = CapProd(max_domain, coeff);
-
-    // Note that we use min/max with zero to disallow "alternative" terms and
-    // be sure that we cannot have an overflow if we do the computation in a
-    // different order.
-    sum_min = CapAdd(sum_min, std::min(int64_t{0}, std::min(prod1, prod2)));
-    sum_max = CapAdd(sum_max, std::max(int64_t{0}, std::max(prod1, prod2)));
-    for (const int64_t v : {prod1, prod2, sum_min, sum_max}) {
-      if (v == std::numeric_limits<int64_t>::max() ||
-          v == std::numeric_limits<int64_t>::min())
-        return true;
-    }
-  }
-
-  // In addition to computing the min/max possible sum, we also often compare
-  // it with the constraint bounds, so we do not want max - min to overflow.
-  // We might also create an intermediate variable to represent the sum. It
-  if (sum_min < std::numeric_limits<int64_t>::min() / 2) return true;
-  if (sum_max > std::numeric_limits<int64_t>::max() / 2) return true;
-  return false;
-  //  if (sum_min < 0 && sum_min + std::numeric_limits<int64_t>::max() <
-  //  sum_max) {
-  //    return true;
-  //  }
-  //  return false;
-}
-
 int64_t MinOfRef(const CpModelProto& model, int ref) {
   const IntegerVariableProto& var_proto = model.variables(PositiveRef(ref));
   if (RefIsPositive(ref)) {
@@ -292,7 +249,8 @@ std::string ValidateLinearExpression(const CpModelProto& model,
     return absl::StrCat("coeffs_size() != vars_size() in linear expression: ",
                         ProtobufShortDebugString(expr));
   }
-  if (PossibleIntegerOverflow(model, expr, expr.offset())) {
+  if (PossibleIntegerOverflow(model, expr.vars(), expr.coeffs(),
+                              expr.offset())) {
     return absl::StrCat("Possible overflow in linear expression: ",
                         ProtobufShortDebugString(expr));
   }
@@ -319,7 +277,7 @@ std::string ValidateLinearConstraint(const CpModelProto& model,
                         ProtobufShortDebugString(ct));
   }
   const LinearConstraintProto& arg = ct.linear();
-  if (PossibleIntegerOverflow(model, arg)) {
+  if (PossibleIntegerOverflow(model, arg.vars(), arg.coeffs())) {
     return "Possible integer overflow in constraint: " +
            ProtobufDebugString(ct);
   }
@@ -421,7 +379,8 @@ std::string ValidateElementConstraint(const CpModelProto& model,
   overflow_detection.add_coeffs(-1);
   for (const int ref : element.vars()) {
     overflow_detection.set_vars(1, ref);
-    if (PossibleIntegerOverflow(model, overflow_detection)) {
+    if (PossibleIntegerOverflow(model, overflow_detection.vars(),
+                                overflow_detection.coeffs())) {
       return absl::StrCat(
           "Domain of the variables involved in element constraint may cause "
           "overflow",
@@ -612,7 +571,8 @@ std::string ValidateIntervalConstraint(const CpModelProto& model,
   RETURN_IF_NOT_EMPTY(ValidateLinearExpression(model, arg.end()));
   AppendToOverflowValidator(arg.end(), &for_overflow_validation);
 
-  if (PossibleIntegerOverflow(model, for_overflow_validation,
+  if (PossibleIntegerOverflow(model, for_overflow_validation.vars(),
+                              for_overflow_validation.coeffs(),
                               for_overflow_validation.offset())) {
     return absl::StrCat("Possible overflow in interval: ",
                         ProtobufShortDebugString(ct.interval()));
@@ -760,7 +720,7 @@ std::string ValidateObjective(const CpModelProto& model,
                           " in objective: ", ProtobufShortDebugString(obj));
     }
   }
-  if (PossibleIntegerOverflow(model, obj)) {
+  if (PossibleIntegerOverflow(model, obj.vars(), obj.coeffs())) {
     return "Possible integer overflow in objective: " +
            ProtobufDebugString(obj);
   }
@@ -887,6 +847,42 @@ std::string ValidateSolutionHint(const CpModelProto& model) {
 
 }  // namespace
 
+bool PossibleIntegerOverflow(const CpModelProto& model,
+                             absl::Span<const int> vars,
+                             absl::Span<const int64_t> coeffs, int64_t offset) {
+  if (offset == std::numeric_limits<int64_t>::min()) return true;
+  int64_t sum_min = -std::abs(offset);
+  int64_t sum_max = +std::abs(offset);
+  for (int i = 0; i < vars.size(); ++i) {
+    const int ref = vars[i];
+    const auto& var_proto = model.variables(PositiveRef(ref));
+    const int64_t min_domain = var_proto.domain(0);
+    const int64_t max_domain = var_proto.domain(var_proto.domain_size() - 1);
+    if (coeffs[i] == std::numeric_limits<int64_t>::min()) return true;
+    const int64_t coeff = RefIsPositive(ref) ? coeffs[i] : -coeffs[i];
+    const int64_t prod1 = CapProd(min_domain, coeff);
+    const int64_t prod2 = CapProd(max_domain, coeff);
+
+    // Note that we use min/max with zero to disallow "alternative" terms and
+    // be sure that we cannot have an overflow if we do the computation in a
+    // different order.
+    sum_min = CapAdd(sum_min, std::min(int64_t{0}, std::min(prod1, prod2)));
+    sum_max = CapAdd(sum_max, std::max(int64_t{0}, std::max(prod1, prod2)));
+    for (const int64_t v : {prod1, prod2, sum_min, sum_max}) {
+      if (v == std::numeric_limits<int64_t>::max() ||
+          v == std::numeric_limits<int64_t>::min())
+        return true;
+    }
+  }
+
+  // In addition to computing the min/max possible sum, we also often compare
+  // it with the constraint bounds, so we do not want max - min to overflow.
+  // We might also create an intermediate variable to represent the sum. It
+  if (sum_min < std::numeric_limits<int64_t>::min() / 2) return true;
+  if (sum_max > std::numeric_limits<int64_t>::max() / 2) return true;
+  return false;
+}
+
 std::string ValidateCpModel(const CpModelProto& model, bool after_presolve) {
   for (int v = 0; v < model.variables_size(); ++v) {
     RETURN_IF_NOT_EMPTY(ValidateIntegerVariable(model, v));

ortools/sat/cp_model_checker.h

@@ -18,6 +18,7 @@
 #include <string>
 #include <vector>
 
+#include "absl/types/span.h"
 #include "ortools/base/integral_types.h"
 #include "ortools/sat/cp_model.pb.h"
 #include "ortools/sat/sat_parameters.pb.h"
@@ -48,6 +49,13 @@ std::string ValidateCpModel(const CpModelProto& model,
 std::string ValidateInputCpModel(const SatParameters& params,
                                  const CpModelProto& model);
 
+// Check if a given linear expression can create overflow. It is exposed to test
+// new constraints created during the presolve.
+bool PossibleIntegerOverflow(const CpModelProto& model,
+                             absl::Span<const int> vars,
+                             absl::Span<const int64_t> coeffs,
+                             int64_t offset = 0);
+
 // Verifies that the given variable assignment is a feasible solution of the
 // given model. The values vector should be in one to one correspondence with
 // the model.variables() list of variables.

ortools/sat/cp_model_lns.h

@@ -19,6 +19,7 @@
 #include <functional>
 #include <string>
 #include <tuple>
+#include <utility>
 #include <vector>
 
 #include "absl/base/thread_annotations.h"

ortools/sat/cp_model_presolve.cc

@@ -18,6 +18,7 @@
 #include <cstdint>
 #include <cstdlib>
 #include <deque>
+#include <iostream>
 #include <limits>
 #include <numeric>
 #include <string>
@@ -31,7 +32,6 @@
 #include "absl/container/flat_hash_map.h"
 #include "absl/container/flat_hash_set.h"
 #include "absl/hash/hash.h"
-#include "absl/meta/type_traits.h"
 #include "absl/numeric/int128.h"
 #include "absl/strings/str_cat.h"
 #include "absl/types/span.h"
@@ -49,6 +49,7 @@
 #include "ortools/sat/cp_model_symmetries.h"
 #include "ortools/sat/cp_model_utils.h"
 #include "ortools/sat/diffn_util.h"
+#include "ortools/sat/diophantine.h"
 #include "ortools/sat/inclusion.h"
 #include "ortools/sat/integer.h"
 #include "ortools/sat/model.h"
@@ -2238,8 +2239,7 @@ bool CpModelPresolver::PresolveLinearOfSizeTwo(ConstraintProto* ct) {
       // In this case, we can solve the diophantine equation, and write
       // both x and y in term of a new affine representative z.
       //
-      // Note that PresolveLinearEqualityWithModularInverse() will have the
-      // same effect.
+      // Note that PresolveLinearEqualityWithModulo() will have the same effect.
       //
       // We can also decide to fully expand the equality if the variables
       // are fully encoded.
@@ -2325,6 +2325,147 @@ bool CpModelPresolver::PresolveSmallLinear(ConstraintProto* ct) {
   return false;
 }
 
+bool CpModelPresolver::PresolveDiophantine(ConstraintProto* ct) {
+  if (ct->constraint_case() != ConstraintProto::kLinear) return false;
+  if (ct->linear().vars().size() <= 1) return false;
+
+  if (context_->ModelIsUnsat()) return false;
+
+  const LinearConstraintProto& linear_constraint = ct->linear();
+  if (linear_constraint.domain_size() != 2) return false;
+  if (linear_constraint.domain(0) != linear_constraint.domain(1)) return false;
+
+  std::vector<int64_t> lbs(linear_constraint.vars_size());
+  std::vector<int64_t> ubs(linear_constraint.vars_size());
+  for (int i = 0; i < linear_constraint.vars_size(); ++i) {
+    lbs[i] = context_->MinOf(linear_constraint.vars(i));
+    ubs[i] = context_->MaxOf(linear_constraint.vars(i));
+  }
+  DiophantineSolution diophantine_solution = SolveDiophantine(
+      linear_constraint.coeffs(), linear_constraint.domain(0), lbs, ubs);
+
+  if (!diophantine_solution.has_solutions) {
+    context_->UpdateRuleStats("diophantine: equality has no solutions");
+    return MarkConstraintAsFalse(ct);
+  }
+  if (diophantine_solution.no_reformulation_needed) return false;
+  // Only first coefficients of kernel_basis elements and special_solution could
+  // overflow int64_t due to the reduction applied in SolveDiophantineEquation,
+  for (const std::vector<absl::int128>& b : diophantine_solution.kernel_basis) {
+    if (!IsNegatableInt64(b[0])) {
+      context_->UpdateRuleStats(
+          "diophantine: couldn't apply due to int64_t overflow");
+      return false;
+    }
+  }
+  if (!IsNegatableInt64(diophantine_solution.special_solution[0])) {
+    context_->UpdateRuleStats(
+        "diophantine: couldn't apply due to int64_t overflow");
+    return false;
+  }
+
+  const int num_replaced_variables =
+      static_cast<int>(diophantine_solution.special_solution.size());
+  const int num_new_variables =
+      static_cast<int>(diophantine_solution.kernel_vars_lbs.size());
+  DCHECK_EQ(num_new_variables + 1, num_replaced_variables);
+  for (int i = 0; i < num_new_variables; ++i) {
+    if (!IsNegatableInt64(diophantine_solution.kernel_vars_lbs[i]) ||
+        !IsNegatableInt64(diophantine_solution.kernel_vars_ubs[i])) {
+      context_->UpdateRuleStats(
+          "diophantine: couldn't apply due to int64_t overflow");
+      return false;
+    }
+  }
+  // TODO(user, demonet): Make sure the newly generated linear constraint
+  // satisfy our no-overflow precondition on the min/max activity.
+  // We should check that the model still satisfy conditions in
+
+  // Create new variables.
+  std::vector<int> new_variables(num_new_variables);
+  for (int i = 0; i < num_new_variables; ++i) {
+    new_variables[i] = context_->working_model->variables_size();
+    IntegerVariableProto* var = context_->working_model->add_variables();
+    var->add_domain(
+        static_cast<int64_t>(diophantine_solution.kernel_vars_lbs[i]));
+    var->add_domain(
+        static_cast<int64_t>(diophantine_solution.kernel_vars_ubs[i]));
+    if (!ct->name().empty()) {
+      var->set_name(absl::StrCat("u_diophantine_", ct->name(), "_", i));
+    }
+  }
+
+  // For i = 0, ..., num_replaced_variables - 1, creates
+  //  x[i] = special_solution[i]
+  //        + sum(kernel_basis[k][i]*y[k], max(1, i) <= k < vars.size - 1)
+  // where:
+  //  y[k] is the newly created variable if 0 <= k < num_new_variables
+  //  y[k] = x[index_permutation[k + 1]] otherwise.
+  for (int i = 0; i < num_replaced_variables; ++i) {
+    ConstraintProto* identity = context_->working_model->add_constraints();
+    LinearConstraintProto* lin = identity->mutable_linear();
+    if (!ct->name().empty()) {
+      identity->set_name(absl::StrCat("c_diophantine_", ct->name(), "_", i));
+    }
+    *identity->mutable_enforcement_literal() = ct->enforcement_literal();
+    lin->add_vars(
+        linear_constraint.vars(diophantine_solution.index_permutation[i]));
+    lin->add_coeffs(1);
+    lin->add_domain(
+        static_cast<int64_t>(diophantine_solution.special_solution[i]));
+    lin->add_domain(
+        static_cast<int64_t>(diophantine_solution.special_solution[i]));
+    for (int j = std::max(1, i); j < num_replaced_variables; ++j) {
+      lin->add_vars(new_variables[j - 1]);
+      lin->add_coeffs(
+          -static_cast<int64_t>(diophantine_solution.kernel_basis[j - 1][i]));
+    }
+    for (int j = num_replaced_variables; j < linear_constraint.vars_size();
+         ++j) {
+      lin->add_vars(
+          linear_constraint.vars(diophantine_solution.index_permutation[j]));
+      lin->add_coeffs(
+          -static_cast<int64_t>(diophantine_solution.kernel_basis[j - 1][i]));
+    }
+    if (PossibleIntegerOverflow(*(context_->working_model), lin->vars(),
+                                lin->coeffs())) {
+      context_->UpdateRuleStats(
+          "diophantine: couldn't apply due to overflowing activity of new "
+          "constraints");
+      // Cancel working_model changes.
+      context_->working_model->mutable_constraints()->DeleteSubrange(
+          context_->working_model->constraints_size() - i, i);
+      context_->working_model->mutable_variables()->DeleteSubrange(
+          context_->working_model->variables_size() - num_new_variables,
+          num_new_variables);
+      return false;
+    }
+  }
+  context_->InitializeNewDomains();
+
+  if (VLOG_IS_ON(2)) {
+    std::string log_eq = absl::StrCat(linear_constraint.domain(0), " = ");
+    const int terms_to_show = std::min<int>(15, linear_constraint.vars_size());
+    for (int i = 0; i < terms_to_show; ++i) {
+      if (i > 0) absl::StrAppend(&log_eq, " + ");
+      absl::StrAppend(
+          &log_eq,
+          linear_constraint.coeffs(diophantine_solution.index_permutation[i]),
+          " x",
+          linear_constraint.vars(diophantine_solution.index_permutation[i]));
+    }
+    if (terms_to_show < linear_constraint.vars_size()) {
+      absl::StrAppend(&log_eq, "+ ... (", linear_constraint.vars_size(),
+                      " terms)");
+    }
+    VLOG(2) << "[Diophantine] " << log_eq;
+  }
+
+  context_->UpdateRuleStats("diophantine: reformulated equality");
+  context_->UpdateNewConstraintsVariableUsage();
+  return RemoveConstraint(ct);
+}
+
 // This tries to decompose the constraint into coeff * part1 + part2 and show
 // that the value that part2 take is not important, thus the constraint can
 // only be transformed on a constraint on the first part.
@@ -6465,6 +6606,10 @@ bool CpModelPresolver::PresolveOneConstraint(int c) {
         }
       }
 
+      if (PresolveDiophantine(ct)) {
+        context_->UpdateConstraintVariableUsage(c);
+      }
+
       TryToReduceCoefficientsOfLinearConstraint(c, ct);
       return false;
     }