diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index dc5826175a..9f011ec4ee 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -937,184 +937,137 @@ bool CpModelPresolver::PresolveIntProd(ConstraintProto* ct) { changed |= CanonicalizeLinearExpression(*ct, &exp); } - // Remove constant terms. - { - int64_t constant_factor = 1; - int new_size = 0; - for (int i = 0; i < ct->int_prod().exprs_size(); ++i) { - LinearExpressionProto expr = ct->int_prod().exprs(i); - if (context_->IsFixed(expr)) { - constant_factor = CapProd(constant_factor, context_->FixedValue(expr)); - continue; - } else { - *proto->mutable_exprs(new_size++) = expr; - } - } - - if (new_size < ct->int_prod().exprs_size()) { + // Remove constant expressions. + int64_t constant_factor = 1; + int new_size = 0; + for (int i = 0; i < ct->int_prod().exprs().size(); ++i) { + LinearExpressionProto expr = ct->int_prod().exprs(i); + if (context_->IsFixed(expr)) { context_->UpdateRuleStats("int_prod: removed constant expressions."); - proto->mutable_exprs()->erase(proto->mutable_exprs()->begin() + new_size, - proto->mutable_exprs()->end()); - } - - if (constant_factor == 0) { - context_->UpdateRuleStats("int_prod: multiplication by zero"); - if (!context_->IntersectDomainWith(ct->int_prod().target(), Domain(0))) { - return false; + changed = true; + constant_factor = CapProd(constant_factor, context_->FixedValue(expr)); + continue; + } else { + const int64_t coeff = expr.coeffs(0); + const int64_t offset = expr.offset(); + const int64_t gcd = + MathUtil::GCD64(static_cast(std::abs(coeff)), + static_cast(std::abs(offset))); + if (gcd != 1) { + constant_factor = CapProd(constant_factor, gcd); + expr.set_coeffs(0, coeff / gcd); + expr.set_offset(offset / gcd); } - return RemoveConstraint(ct); } + *proto->mutable_exprs(new_size++) = expr; + } + proto->mutable_exprs()->erase(proto->mutable_exprs()->begin() + new_size, + proto->mutable_exprs()->end()); - if (ct->int_prod().exprs().empty()) { - if (!context_->IntersectDomainWith(ct->int_prod().target(), - Domain(constant_factor))) { - return false; - } - context_->UpdateRuleStats("int_prod: constant product"); - return RemoveConstraint(ct); - } - - // We cannot have a int_max or int_min constant factor as this would have - // been prevented by the checker. - DCHECK_NE(constant_factor, std::numeric_limits::min()); - DCHECK_NE(constant_factor, std::numeric_limits::max()); - - // Replace by linear! - if (ct->int_prod().exprs().size() == 1) { - LinearConstraintProto* const lin = - context_->working_model->add_constraints()->mutable_linear(); - lin->add_domain(0); - lin->add_domain(0); - AddLinearExpressionToLinearConstraint(ct->int_prod().target(), 1, lin); - AddLinearExpressionToLinearConstraint(ct->int_prod().exprs(0), - -constant_factor, lin); - context_->UpdateNewConstraintsVariableUsage(); - context_->UpdateRuleStats("int_prod: linearize product by constant."); - return RemoveConstraint(ct); + if (ct->int_prod().exprs().empty()) { + if (!context_->IntersectDomainWith(ct->int_prod().target(), + Domain(constant_factor))) { + return false; } + context_->UpdateRuleStats("int_prod: constant product"); + return RemoveConstraint(ct); } - // Currently, the product only supports 2 terms. If one or both are fixed, - // this should have been taken care of. - CHECK_EQ(2, proto->exprs_size()); - - int64_t product_of_expr_gcds = 1; - std::vector expr_gcds; - for (const LinearExpressionProto& expr : ct->int_prod().exprs()) { - const int64_t gcd = - MathUtil::GCD64(static_cast(std::abs(expr.coeffs(0))), - static_cast(std::abs(expr.offset()))); - product_of_expr_gcds = CapProd(product_of_expr_gcds, gcd); - expr_gcds.push_back(gcd); + if (constant_factor == 0) { + context_->UpdateRuleStats("int_prod: multiplication by zero"); + if (!context_->IntersectDomainWith(ct->int_prod().target(), Domain(0))) { + return false; + } + return RemoveConstraint(ct); } // In this case, the only possible value that fit in the domains is zero. // We will check for UNSAT if zero is not achievable by the rhs below. - if (product_of_expr_gcds == std::numeric_limits::min() || - product_of_expr_gcds == std::numeric_limits::max()) { + if (constant_factor == std::numeric_limits::min() || + constant_factor == std::numeric_limits::max()) { + context_->UpdateRuleStats("int_prod: overflow if non zero"); if (!context_->IntersectDomainWith(ct->int_prod().target(), Domain(0))) { return false; } - // Make sure the target is 0. - proto->clear_target(); - proto->mutable_target()->set_offset(0); - - // Simplify all terms. Actually, we just want one of them to be zero. - for (int i = 0; i < proto->exprs_size(); ++i) { - const int64_t expr_gcd = expr_gcds[i]; - if (expr_gcd == 1) continue; - - LinearExpressionProto* expr = proto->mutable_exprs(i); - DCHECK_EQ(1, expr->coeffs_size()); - expr->set_coeffs(0, expr->coeffs(0) / expr_gcd); - expr->set_offset(expr->offset() / expr_gcd); - } - product_of_expr_gcds = 1; - context_->UpdateRuleStats("int_prod: overflow if non zero"); + constant_factor = 1; } - if (product_of_expr_gcds > 1) { + // Replace by linear! + if (ct->int_prod().exprs().size() == 1) { + LinearConstraintProto* const lin = + context_->working_model->add_constraints()->mutable_linear(); + lin->add_domain(0); + lin->add_domain(0); + AddLinearExpressionToLinearConstraint(ct->int_prod().target(), 1, lin); + AddLinearExpressionToLinearConstraint(ct->int_prod().exprs(0), + -constant_factor, lin); + context_->UpdateNewConstraintsVariableUsage(); + context_->UpdateRuleStats("int_prod: linearize product by constant."); + return RemoveConstraint(ct); + } + + if (constant_factor != 1) { + // Lets canonicalize the target by introducing a new variable if necessary. + // + // coeff * X + offset must be a multiple of constant_factor, so + // we can rewrite X so that this property is clear. const LinearExpressionProto old_target = ct->int_prod().target(); + if (!context_->IsFixed(old_target)) { + const int ref = old_target.vars(0); + const int64_t coeff = old_target.coeffs(0); + const int64_t offset = old_target.offset(); + if (!context_->CanonicalizeAffineVariable(ref, coeff, constant_factor, + -offset)) { + return false; + } + } + + // This can happen during CanonicalizeAffineVariable(). if (context_->IsFixed(old_target)) { const int64_t target_value = context_->FixedValue(old_target); - if (target_value % product_of_expr_gcds != 0) { + if (target_value % constant_factor != 0) { return context_->NotifyThatModelIsUnsat( "int_prod: constant factor does not divide constant target"); } - - // Divide target and terms by product_of_expr_gcds. proto->clear_target(); - proto->mutable_target()->set_offset(target_value / product_of_expr_gcds); - for (int i = 0; i < proto->exprs_size(); ++i) { - const int64_t expr_gcd = expr_gcds[i]; - if (expr_gcd == 1) continue; - - LinearExpressionProto* expr = proto->mutable_exprs(i); - DCHECK_EQ(1, expr->coeffs_size()); - expr->set_coeffs(0, expr->coeffs(0) / expr_gcd); - expr->set_offset(expr->offset() / expr_gcd); - } + proto->mutable_target()->set_offset(target_value / constant_factor); context_->UpdateRuleStats( "int_prod: divide product and fixed target by constant factor"); - } else { // Target is not fixed. - int64_t target_coeff = ct->int_prod().target().coeffs(0); - int64_t target_offset = ct->int_prod().target().offset(); - int64_t target_gcd = - MathUtil::GCD64(static_cast(std::abs(target_coeff)), - static_cast(std::abs(target_offset))); - if (target_gcd % product_of_expr_gcds == 0) { // Easy case. - // Divide target by product_of_expr_gcds. - proto->mutable_target()->set_coeffs( - 0, target_coeff / product_of_expr_gcds); - proto->mutable_target()->set_offset(target_offset / - product_of_expr_gcds); - // Divide all terms of the product by their respective gcds. - for (int i = 0; i < proto->exprs_size(); ++i) { - const int64_t expr_gcd = expr_gcds[i]; - if (expr_gcd == 1) continue; + } else { + // We use absl::int128 to be resistant to overflow here. + const AffineRelation::Relation r = + context_->GetAffineRelation(old_target.vars(0)); + const absl::int128 temp_coeff = + absl::int128(old_target.coeffs(0)) * absl::int128(r.coeff); + CHECK_EQ(temp_coeff % absl::int128(constant_factor), 0); + const absl::int128 temp_offset = + absl::int128(old_target.coeffs(0)) * absl::int128(r.offset) + + absl::int128(old_target.offset()); + CHECK_EQ(temp_offset % absl::int128(constant_factor), 0); + const absl::int128 new_coeff = temp_coeff / absl::int128(constant_factor); + const absl::int128 new_offset = + temp_offset / absl::int128(constant_factor); - LinearExpressionProto* expr = proto->mutable_exprs(i); - DCHECK_EQ(1, expr->coeffs_size()); - expr->set_coeffs(0, expr->coeffs(0) / expr_gcd); - expr->set_offset(expr->offset() / expr_gcd); - } - context_->UpdateRuleStats("int_prod: constant factor divides target"); - } else { - const int64_t global_gcd = MathUtil::GCD64( - static_cast(target_gcd), - static_cast(std::abs(product_of_expr_gcds))); - // target_gcd * (coeff * target_var * offset) == - // product_of_expr_gcds * PRODUCT. - // Let's divide the target by global_gcd and spread global_gcd - // across the product terms. - if (global_gcd > 1) { - proto->mutable_target()->set_coeffs(0, target_coeff / global_gcd); - proto->mutable_target()->set_offset(target_offset / global_gcd); - int64_t remainder_gcd = global_gcd; - for (int i = 0; i < proto->exprs_size(); ++i) { - const int64_t expr_gcd = expr_gcds[i]; - if (expr_gcd == 1) continue; - int64_t local_gcd = - MathUtil::GCD64(static_cast(std::abs(remainder_gcd)), - static_cast(std::abs(expr_gcd))); - if (local_gcd == 1) continue; - - remainder_gcd /= local_gcd; - - LinearExpressionProto* expr = proto->mutable_exprs(i); - DCHECK_EQ(1, expr->coeffs_size()); - expr->set_coeffs(0, expr->coeffs(0) / local_gcd); - expr->set_offset(expr->offset() / local_gcd); - } - CHECK_EQ(1, remainder_gcd); // We have consumed global_gcd. - context_->UpdateRuleStats( - "int_prod: divide product and target by gcd"); - } - // TODO(user, fdid): Instead we could "canonicalize" the target by - // creating a new variable such that - // target = product_of_expr_gcds * new_var - // and use new_var here. + // TODO(user): We try to keep coeff/offset small, if this happens, it + // probably means there is no feasible solution involving int64_t and that + // do not causes overflow while evaluating it, but it is hard to be + // exactly sure we are correct here since it depends on the evaluation + // order. Similarly, by introducing intermediate variable we might loose + // solution if this intermediate variable value do not fit on an int64_t. + if (new_coeff > absl::int128(std::numeric_limits::max()) || + new_coeff < absl::int128(std::numeric_limits::min()) || + new_offset > absl::int128(std::numeric_limits::max()) || + new_offset < absl::int128(std::numeric_limits::min())) { + return context_->NotifyThatModelIsUnsat( + "int_prod: overflow during simplification."); } + + // Rewrite the target. + proto->mutable_target()->set_coeffs(0, static_cast(new_coeff)); + proto->mutable_target()->set_vars(0, r.representative); + proto->mutable_target()->set_offset(static_cast(new_offset)); + context_->UpdateRuleStats("int_prod: divide product by constant factor"); + changed = true; } } @@ -1124,12 +1077,12 @@ bool CpModelPresolver::PresolveIntProd(ConstraintProto* ct) { implied = implied.ContinuousMultiplicationBy(context_->DomainSuperSetOf(expr)); } - bool modified = false; + bool domain_modified = false; if (!context_->IntersectDomainWith(ct->int_prod().target(), implied, - &modified)) { + &domain_modified)) { return false; } - if (modified) { + if (domain_modified) { context_->UpdateRuleStats("int_prod: reduced target domain."); } @@ -7591,6 +7544,9 @@ CpSolverStatus CpModelPresolver::Presolve() { } // Exit the loop if the reduction is not so large. + // Hack: to facilitate experiments, if the requested number of iterations + // is large, we always execute all of them. + if (context_->params().max_presolve_iterations() >= 5) continue; if (context_->num_presolve_operations - old_num_presolve_op < 0.8 * (context_->working_model->variables_size() + old_num_non_empty_constraints)) { diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index 38a108a269..a907158761 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -1001,14 +1001,14 @@ bool SquarePropagator::Propagate() { const IntegerValue min_s = integer_trail_->LowerBound(s_); const IntegerValue min_x_square(CapProd(min_x.value(), min_x.value())); if (min_x_square > min_s) { - if (!integer_trail_->Enqueue(s_.GreaterOrEqual(min_x_square), {}, - {x_.GreaterOrEqual(min_x)})) { + if (!integer_trail_->SafeEnqueue(s_.GreaterOrEqual(min_x_square), + {x_.GreaterOrEqual(min_x)})) { return false; } } else if (min_x_square < min_s) { const IntegerValue new_min = CeilSquareRoot(min_s); - if (!integer_trail_->Enqueue( - x_.GreaterOrEqual(new_min), {}, + if (!integer_trail_->SafeEnqueue( + x_.GreaterOrEqual(new_min), {s_.GreaterOrEqual((new_min - 1) * (new_min - 1) + 1)})) { return false; } @@ -1018,14 +1018,14 @@ bool SquarePropagator::Propagate() { const IntegerValue max_s = integer_trail_->UpperBound(s_); const IntegerValue max_x_square(CapProd(max_x.value(), max_x.value())); if (max_x_square < max_s) { - if (!integer_trail_->Enqueue(s_.LowerOrEqual(max_x_square), {}, - {x_.LowerOrEqual(max_x)})) { + if (!integer_trail_->SafeEnqueue(s_.LowerOrEqual(max_x_square), + {x_.LowerOrEqual(max_x)})) { return false; } } else if (max_x_square > max_s) { const IntegerValue new_max = FloorSquareRoot(max_s); - if (!integer_trail_->Enqueue( - x_.LowerOrEqual(new_max), {}, + if (!integer_trail_->SafeEnqueue( + x_.LowerOrEqual(new_max), {s_.LowerOrEqual(IntegerValue(CapProd(new_max.value() + 1, new_max.value() + 1)) - 1)})) { diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 855d3da3f1..126250f675 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -804,14 +804,30 @@ bool PresolveContext::ScaleFloatingPointObjective() { bool PresolveContext::CanonicalizeAffineVariable(int ref, int64_t coeff, int64_t mod, int64_t rhs) { - CHECK_GT(mod, 1); - CHECK_NE(coeff % mod, 0); + CHECK_NE(mod, 0); + CHECK_NE(coeff, 0); + + const int64_t gcd = std::gcd(coeff, mod); + if (gcd != 1) { + if (rhs % gcd != 0) { + return NotifyThatModelIsUnsat( + absl::StrCat("Infeasible ", coeff, " * X = ", rhs, " % ", mod)); + } + coeff /= gcd; + mod /= gcd; + rhs /= gcd; + } + + // We just abort in this case as there is no point introducing a new variable. + if (std::abs(mod) == 1) return true; + int var = ref; if (!RefIsPositive(var)) { var = NegatedRef(ref); coeff = -coeff; rhs = -rhs; } + // From var * coeff % mod = rhs // We have var = mod * X + offset. const int64_t offset = ProductWithModularInverse(coeff, mod, rhs); @@ -820,7 +836,8 @@ bool PresolveContext::CanonicalizeAffineVariable(int ref, int64_t coeff, const Domain new_domain = DomainOf(var).AdditionWith(Domain(-offset)).InverseMultiplicationBy(mod); if (new_domain.IsEmpty()) { - return NotifyThatModelIsUnsat(); + return NotifyThatModelIsUnsat( + "Empty domain in CanonicalizeAffineVariable()"); } if (new_domain.IsFixed()) { UpdateRuleStats("variables: fixed value due to affine relation"); diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 9052860c18..1f4973e8af 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -193,7 +193,7 @@ class PresolveContext { ABSL_MUST_USE_RESULT bool NotifyThatModelIsUnsat( const std::string& message = "") { // TODO(user): Report any explanation for the client in a nicer way? - VLOG(1) << "INFEASIBLE: '" << message << "'"; + SOLVER_LOG(logger_, "INFEASIBLE: '", message, "'"); DCHECK(!is_unsat_); is_unsat_ = true; return false; @@ -238,11 +238,11 @@ class PresolveContext { // Given the relation (X * coeff % mod = rhs % mod), this creates a new // variable so that X = mod * Y + cte. // - // This assumes mod > 1 and coeff % mod != 0 (CHECKed). + // This requires mod != 0 and coeff != 0. // // Note that the new variable will have a canonical domain (i.e. min == 0). - // We also do not create anything if this fixes the given variable. - // Returns false if the model is infeasible. + // We also do not create anything if this fixes the given variable or the + // relation simplifies. Returns false if the model is infeasible. bool CanonicalizeAffineVariable(int ref, int64_t coeff, int64_t mod, int64_t rhs);