diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc
index e1b7a3f89c..d1784c70bc 100644
--- a/ortools/sat/cp_model_presolve.cc
+++ b/ortools/sat/cp_model_presolve.cc
@@ -809,14 +809,17 @@ void CpModelPresolver::DivideLinearByGcd(ConstraintProto* ct) {
 }
 
 bool CpModelPresolver::CanonicalizeLinear(ConstraintProto* ct) {
-  if (context_->ModelIsUnsat()) return false;
+  if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear ||
+      context_->ModelIsUnsat()) {
+    return false;
+  }
 
   // First regroup the terms on the same variables and sum the fixed ones.
   //
   // TODO(user): move terms in context to reuse its memory? Add a quick pass
   // to skip most of the work below if the constraint is already in canonical
   // form (strictly increasing var, no-fixed var, gcd = 1).
-  std::vector<std::pair<int, int64>> terms;
+  tmp_terms_.clear();
   const bool was_affine = gtl::ContainsKey(context_->affine_constraints, ct);
 
   int64 sum_of_fixed_terms = 0;
@@ -840,9 +843,9 @@ bool CpModelPresolver::CanonicalizeLinear(ConstraintProto* ct) {
         remapped = true;
         sum_of_fixed_terms += coeff * r.offset;
       }
-      terms.push_back({r.representative, coeff * r.coeff});
+      tmp_terms_.push_back({r.representative, coeff * r.coeff});
     } else {
-      terms.push_back({var, coeff});
+      tmp_terms_.push_back({var, coeff});
     }
   }
 
@@ -854,10 +857,10 @@ bool CpModelPresolver::CanonicalizeLinear(ConstraintProto* ct) {
 
   ct->mutable_linear()->clear_vars();
   ct->mutable_linear()->clear_coeffs();
-  std::sort(terms.begin(), terms.end());
+  std::sort(tmp_terms_.begin(), tmp_terms_.end());
   int current_var = 0;
   int64 current_coeff = 0;
-  for (const auto entry : terms) {
+  for (const auto entry : tmp_terms_) {
     CHECK(RefIsPositive(entry.first));
     if (entry.first == current_var) {
       current_coeff += entry.second;
@@ -889,7 +892,10 @@ bool CpModelPresolver::CanonicalizeLinear(ConstraintProto* ct) {
 }
 
 bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) {
-  if (context_->ModelIsUnsat()) return false;
+  if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear ||
+      context_->ModelIsUnsat()) {
+    return false;
+  }
 
   const bool was_affine = gtl::ContainsKey(context_->affine_constraints, ct);
   if (was_affine) return false;
@@ -959,19 +965,17 @@ bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) {
       // Special case: If the objective was a single variable, we can transfer
       // the domain of var to the objective, and just completely remove this
       // equality constraint like it is done in ExpandObjective().
-      //
-      // TODO(user): handle coeff other than one. Remove duplication with
-      // expand objective.
-      const CpObjectiveProto& obj = context_->working_model->objective();
-      if (obj.vars().size() == 1 && obj.vars(0) == var && obj.coeffs(0) == 1) {
-        if (!obj.domain().empty()) {
-          const Domain obj_domain = ReadDomainFromProto(obj);
-          if (!context_->IntersectDomainWith(var, obj_domain)) {
-            return true;
-          }
+      if (context_->ObjectiveMap().size() == 1 &&
+          context_->ObjectiveMap().contains(var)) {
+        if (!context_->IntersectDomainWith(
+                var, context_->ObjectiveDomain().InverseMultiplicationBy(
+		    gtl::FindOrDie(context_->ObjectiveMap(), var)))) {
+          return true;
         }
-        FillDomainInProto(context_->DomainOf(var),
-                          context_->working_model->mutable_objective());
+
+        // This makes sure the domain of var is propagated back to the
+        // objective.
+        context_->CanonicalizeObjective();
         context_->UpdateRuleStats("linear: singleton column define objective.");
         context_->SubstituteVariableInObjective(var, coeff, *ct);
         *(context_->mapping_model->add_constraints()) = *ct;
@@ -1010,14 +1014,16 @@ bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) {
   return true;
 }
 
-bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) {
-  if (context_->ModelIsUnsat()) return false;
-
-  Domain rhs = ReadDomainFromProto(ct->linear());
+bool CpModelPresolver::PresolveSmallLinear(ConstraintProto* ct) {
+  if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear ||
+      context_->ModelIsUnsat()) {
+    return false;
+  }
 
   // Empty constraint?
   if (ct->linear().vars().empty()) {
     context_->UpdateRuleStats("linear: empty");
+    const Domain rhs = ReadDomainFromProto(ct->linear());
     if (rhs.Contains(0)) {
       return RemoveConstraint(ct);
     } else {
@@ -1025,214 +1031,30 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) {
     }
   }
 
+  if (HasEnforcementLiteral(*ct)) return false;
+
   // Size one constraint?
-  if (ct->linear().vars().size() == 1 && !HasEnforcementLiteral(*ct)) {
+  if (ct->linear().vars().size() == 1) {
     const int64 coeff = RefIsPositive(ct->linear().vars(0))
                             ? ct->linear().coeffs(0)
                             : -ct->linear().coeffs(0);
     context_->UpdateRuleStats("linear: size one");
     const int var = PositiveRef(ct->linear().vars(0));
-    if (coeff == 1) {
-      if (!context_->IntersectDomainWith(var, rhs)) {
-        return true;
-      }
-    } else {
-      DCHECK_EQ(coeff, -1);  // Because of the GCD above.
-      if (!context_->IntersectDomainWith(var, rhs.Negation())) {
-        return true;
-      }
+    const Domain rhs = ReadDomainFromProto(ct->linear());
+    if (!context_->IntersectDomainWith(var,
+                                       rhs.InverseMultiplicationBy(coeff))) {
+      return true;
     }
     return RemoveConstraint(ct);
   }
 
-  // Compute the implied rhs bounds from the variable ones.
-  auto& term_domains = context_->tmp_term_domains;
-  auto& left_domains = context_->tmp_left_domains;
-  const int num_vars = ct->linear().vars_size();
-  term_domains.resize(num_vars + 1);
-  left_domains.resize(num_vars + 1);
-  left_domains[0] = Domain(0);
-  for (int i = 0; i < num_vars; ++i) {
-    const int var = PositiveRef(ct->linear().vars(i));
-    const int64 coeff = ct->linear().coeffs(i);
-    term_domains[i] = context_->DomainOf(var).MultiplicationBy(coeff);
-    left_domains[i + 1] =
-        left_domains[i].AdditionWith(term_domains[i]).RelaxIfTooComplex();
-  }
-  const Domain& implied_rhs = left_domains[num_vars];
-
-  // Abort if trivial.
-  if (implied_rhs.IsIncludedIn(rhs)) {
-    context_->UpdateRuleStats("linear: always true");
-    return RemoveConstraint(ct);
-  }
-
-  // Incorporate the implied rhs information.
-  rhs = rhs.SimplifyUsingImpliedDomain(implied_rhs);
-  if (rhs.IsEmpty()) {
-    context_->UpdateRuleStats("linear: infeasible");
-    return MarkConstraintAsFalse(ct);
-  }
-  if (rhs != ReadDomainFromProto(ct->linear())) {
-    context_->UpdateRuleStats("linear: simplified rhs");
-  }
-  FillDomainInProto(rhs, ct->mutable_linear());
-
-  // Propagate the variable bounds.
-  if (ct->enforcement_literal().size() <= 1) {
-    bool new_bounds = false;
-    Domain new_domain;
-    Domain right_domain(0, 0);
-    term_domains[num_vars] = rhs.Negation();
-    for (int i = num_vars - 1; i >= 0; --i) {
-      const int64 var_coeff = ct->linear().coeffs(i);
-      right_domain =
-          right_domain.AdditionWith(term_domains[i + 1]).RelaxIfTooComplex();
-      new_domain = left_domains[i]
-                       .AdditionWith(right_domain)
-                       .InverseMultiplicationBy(-var_coeff);
-
-      if (ct->enforcement_literal().size() == 1) {
-        // We cannot push the new domain, but we can add some deduction.
-        const int var = ct->linear().vars(i);
-        CHECK(RefIsPositive(var));
-        if (!context_->DomainOfVarIsIncludedIn(var, new_domain)) {
-          context_->deductions.AddDeduction(ct->enforcement_literal(0), var,
-                                            new_domain);
-        }
-        continue;
-      }
-
-      bool domain_modified = false;
-      if (!context_->IntersectDomainWith(ct->linear().vars(i), new_domain,
-                                         &domain_modified)) {
-        return true;
-      }
-      if (domain_modified) {
-        new_bounds = true;
-      }
-
-      // Can we perform some substitution?
-      //
-      // TODO(user): there is no guarante we will not miss some since we might
-      // not reprocess a constraint once other have been deleted.
-
-      // Skip affine constraint. It is more efficient to substitute them lazily
-      // when we process other constraints. Note that if we relax the fact that
-      // we substitute only equalities, we can deal with inequality of size 2
-      // here.
-      if (ct->linear().vars().size() <= 2) continue;
-
-      // TODO(user): We actually do not need a strict equality when
-      // keep_all_feasible_solutions is false, but that simplifies things as the
-      // SubstituteVariable() function cannot fail this way.
-      if (rhs.Min() != rhs.Max()) continue;
-
-      // Only consider "implied free" variables. Note that the coefficient of
-      // magnitude 1 is important otherwise we can't easily remove the
-      // constraint since the fact that the sum of the other terms must be a
-      // multiple of coeff will not be enforced anymore.
-      const int var = ct->linear().vars(i);
-      if (context_->DomainOf(var) != new_domain) continue;
-      if (std::abs(var_coeff) != 1) continue;
-      if (options_.parameters.presolve_substitution_level() <= 0) continue;
-
-      // NOTE: The mapping doesn't allow us to remove a variable if
-      // 'keep_all_feasible_solutions' is true.
-      if (context_->keep_all_feasible_solutions) continue;
-
-      bool is_in_objective = false;
-      if (context_->var_to_constraints[var].contains(-1)) {
-        is_in_objective = true;
-        DCHECK(VarFoundInObjective(
-            var, context_->working_model->mutable_objective()));
-      }
-
-      // Only consider low degree columns.
-      int col_size = context_->var_to_constraints[var].size();
-      if (is_in_objective) col_size--;
-      const int row_size = ct->linear().vars_size();
-
-      // This is actually an upper bound on the number of entries added since
-      // some of them might already be present.
-      const int num_entries_added = (row_size - 1) * (col_size - 1);
-      const int num_entries_removed = col_size + row_size - 1;
-
-      if (num_entries_added > num_entries_removed) {
-        continue;
-      }
-
-      // Check pre-conditions on all the constraints in which this variable
-      // appear. Basically they must all be linear and not already used for
-      // affine relation.
-      std::vector<int> others;
-      for (const int c : context_->var_to_constraints[var]) {
-        if (c == -1) continue;
-        if (context_->working_model->mutable_constraints(c) == ct) continue;
-        others.push_back(c);
-      }
-      bool abort = false;
-      for (const int c : others) {
-        if (context_->working_model->constraints(c).constraint_case() !=
-            ConstraintProto::ConstraintCase::kLinear) {
-          abort = true;
-          break;
-        }
-        if (context_->affine_constraints.contains(
-                &context_->working_model->constraints(c))) {
-          abort = true;
-          break;
-        }
-        for (const int ref :
-             context_->working_model->constraints(c).enforcement_literal()) {
-          if (PositiveRef(ref) == var) {
-            abort = true;
-            break;
-          }
-        }
-        if (abort) break;
-      }
-      if (abort) continue;
-
-      // Do the actual substitution.
-      for (const int c : others) {
-        // TODO(user): In some corner cases, this might create integer overflow
-        // issues. The danger is limited since the range of the linear
-        // expression used in the definition do not exceed the domain of the
-        // variable we substitute.
-        SubstituteVariable(var, var_coeff, *ct,
-                           context_->working_model->mutable_constraints(c));
-
-        // TODO(user): We should re-enqueue these constraints for presolve.
-        context_->UpdateConstraintVariableUsage(c);
-      }
-
-      // Substitute in objective.
-      if (is_in_objective) {
-        context_->SubstituteVariableInObjective(var, var_coeff, *ct);
-      }
-
-      context_->UpdateRuleStats(
-          absl::StrCat("linear: variable substitution ", others.size()));
-
-      // The variable now only appear in its definition and we can remove it
-      // because it was implied free.
-      CHECK_EQ(context_->var_to_constraints[var].size(), 1);
-      *context_->mapping_model->add_constraints() = *ct;
-      return RemoveConstraint(ct);
-    }
-    if (new_bounds) {
-      context_->UpdateRuleStats("linear: reduced variable domains");
-    }
-  }
-
   // Detect affine relation.
   //
   // TODO(user): it might be better to first add only the affine relation with
   // a coefficient of magnitude 1, and later the one with larger coeffs.
-  const bool was_affine = context_->affine_constraints.contains(ct);
-  if (!was_affine && !HasEnforcementLiteral(*ct)) {
+  if (!context_->affine_constraints.contains(ct)) {
     const LinearConstraintProto& arg = ct->linear();
+    const Domain rhs = ReadDomainFromProto(ct->linear());
     const int64 rhs_min = rhs.Min();
     const int64 rhs_max = rhs.Max();
     if (rhs_min == rhs_max && arg.vars_size() == 2) {
@@ -1255,6 +1077,192 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) {
   return false;
 }
 
+bool CpModelPresolver::PropagateDomainsInLinear(ConstraintProto* ct) {
+  if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear ||
+      context_->ModelIsUnsat()) {
+    return false;
+  }
+
+  // Compute the implied rhs bounds from the variable ones.
+  auto& term_domains = context_->tmp_term_domains;
+  auto& left_domains = context_->tmp_left_domains;
+  const int num_vars = ct->linear().vars_size();
+  term_domains.resize(num_vars + 1);
+  left_domains.resize(num_vars + 1);
+  left_domains[0] = Domain(0);
+  for (int i = 0; i < num_vars; ++i) {
+    const int var = PositiveRef(ct->linear().vars(i));
+    const int64 coeff = ct->linear().coeffs(i);
+    term_domains[i] = context_->DomainOf(var).MultiplicationBy(coeff);
+    left_domains[i + 1] =
+        left_domains[i].AdditionWith(term_domains[i]).RelaxIfTooComplex();
+  }
+  const Domain& implied_rhs = left_domains[num_vars];
+
+  // Abort if trivial.
+  const Domain old_rhs = ReadDomainFromProto(ct->linear());
+  if (implied_rhs.IsIncludedIn(old_rhs)) {
+    context_->UpdateRuleStats("linear: always true");
+    return RemoveConstraint(ct);
+  }
+
+  // Incorporate the implied rhs information.
+  const Domain rhs = old_rhs.SimplifyUsingImpliedDomain(implied_rhs);
+  if (rhs.IsEmpty()) {
+    context_->UpdateRuleStats("linear: infeasible");
+    return MarkConstraintAsFalse(ct);
+  }
+  if (rhs != old_rhs) {
+    context_->UpdateRuleStats("linear: simplified rhs");
+  }
+  FillDomainInProto(rhs, ct->mutable_linear());
+
+  // Propagate the variable bounds.
+  if (ct->enforcement_literal().size() > 1) return false;
+
+  bool new_bounds = false;
+  Domain new_domain;
+  Domain right_domain(0, 0);
+  term_domains[num_vars] = rhs.Negation();
+  for (int i = num_vars - 1; i >= 0; --i) {
+    const int64 var_coeff = ct->linear().coeffs(i);
+    right_domain =
+        right_domain.AdditionWith(term_domains[i + 1]).RelaxIfTooComplex();
+    new_domain = left_domains[i]
+                     .AdditionWith(right_domain)
+                     .InverseMultiplicationBy(-var_coeff);
+
+    if (ct->enforcement_literal().size() == 1) {
+      // We cannot push the new domain, but we can add some deduction.
+      const int var = ct->linear().vars(i);
+      CHECK(RefIsPositive(var));
+      if (!context_->DomainOfVarIsIncludedIn(var, new_domain)) {
+        context_->deductions.AddDeduction(ct->enforcement_literal(0), var,
+                                          new_domain);
+      }
+      continue;
+    }
+
+    if (!context_->IntersectDomainWith(ct->linear().vars(i), new_domain,
+                                       &new_bounds)) {
+      return true;
+    }
+
+    // Can we perform some substitution?
+    //
+    // TODO(user): there is no guarante we will not miss some since we might
+    // not reprocess a constraint once other have been deleted.
+
+    // Skip affine constraint. It is more efficient to substitute them lazily
+    // when we process other constraints. Note that if we relax the fact that
+    // we substitute only equalities, we can deal with inequality of size 2
+    // here.
+    if (ct->linear().vars().size() <= 2) continue;
+
+    // TODO(user): We actually do not need a strict equality when
+    // keep_all_feasible_solutions is false, but that simplifies things as the
+    // SubstituteVariable() function cannot fail this way.
+    if (rhs.Min() != rhs.Max()) continue;
+
+    // Only consider "implied free" variables. Note that the coefficient of
+    // magnitude 1 is important otherwise we can't easily remove the
+    // constraint since the fact that the sum of the other terms must be a
+    // multiple of coeff will not be enforced anymore.
+    const int var = ct->linear().vars(i);
+    if (context_->DomainOf(var) != new_domain) continue;
+    if (std::abs(var_coeff) != 1) continue;
+    if (options_.parameters.presolve_substitution_level() <= 0) continue;
+
+    // NOTE: The mapping doesn't allow us to remove a variable if
+    // 'keep_all_feasible_solutions' is true.
+    if (context_->keep_all_feasible_solutions) continue;
+
+    bool is_in_objective = false;
+    if (context_->var_to_constraints[var].contains(-1)) {
+      is_in_objective = true;
+      DCHECK(context_->ObjectiveMap().contains(var));
+    }
+
+    // Only consider low degree columns.
+    int col_size = context_->var_to_constraints[var].size();
+    if (is_in_objective) col_size--;
+    const int row_size = ct->linear().vars_size();
+
+    // This is actually an upper bound on the number of entries added since
+    // some of them might already be present.
+    const int num_entries_added = (row_size - 1) * (col_size - 1);
+    const int num_entries_removed = col_size + row_size - 1;
+
+    if (num_entries_added > num_entries_removed) {
+      continue;
+    }
+
+    // Check pre-conditions on all the constraints in which this variable
+    // appear. Basically they must all be linear and not already used for
+    // affine relation.
+    std::vector<int> others;
+    for (const int c : context_->var_to_constraints[var]) {
+      if (c == -1) continue;
+      if (context_->working_model->mutable_constraints(c) == ct) continue;
+      others.push_back(c);
+    }
+    bool abort = false;
+    for (const int c : others) {
+      if (context_->working_model->constraints(c).constraint_case() !=
+          ConstraintProto::ConstraintCase::kLinear) {
+        abort = true;
+        break;
+      }
+      if (context_->affine_constraints.contains(
+              &context_->working_model->constraints(c))) {
+        abort = true;
+        break;
+      }
+      for (const int ref :
+           context_->working_model->constraints(c).enforcement_literal()) {
+        if (PositiveRef(ref) == var) {
+          abort = true;
+          break;
+        }
+      }
+      if (abort) break;
+    }
+    if (abort) continue;
+
+    // Do the actual substitution.
+    for (const int c : others) {
+      // TODO(user): In some corner cases, this might create integer overflow
+      // issues. The danger is limited since the range of the linear
+      // expression used in the definition do not exceed the domain of the
+      // variable we substitute.
+      SubstituteVariable(var, var_coeff, *ct,
+                         context_->working_model->mutable_constraints(c));
+
+      // TODO(user): We should re-enqueue these constraints for presolve.
+      context_->UpdateConstraintVariableUsage(c);
+    }
+
+    // Substitute in objective.
+    if (is_in_objective) {
+      context_->SubstituteVariableInObjective(var, var_coeff, *ct);
+    }
+
+    context_->UpdateRuleStats(
+        absl::StrCat("linear: variable substitution ", others.size()));
+
+    // The variable now only appear in its definition and we can remove it
+    // because it was implied free.
+    CHECK_EQ(context_->var_to_constraints[var].size(), 1);
+    *context_->mapping_model->add_constraints() = *ct;
+    return RemoveConstraint(ct);
+  }
+  if (new_bounds) {
+    context_->UpdateRuleStats("linear: reduced variable domains");
+  }
+
+  return false;
+}
+
 // Identify Boolean variable that makes the constraint always true when set to
 // true or false. Moves such literal to the constraint enforcement literals
 // list.
@@ -1265,7 +1273,10 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) {
 // by using the same logic (i.e. real coefficient strengthening).
 void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint(
     ConstraintProto* ct) {
-  if (context_->ModelIsUnsat()) return;
+  if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear ||
+      context_->ModelIsUnsat()) {
+    return;
+  }
   if (context_->affine_constraints.contains(ct)) return;
 
   const LinearConstraintProto& arg = ct->linear();
@@ -1406,7 +1417,7 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint(
 void CpModelPresolver::ExtractAtMostOneFromLinear(ConstraintProto* ct) {
   if (context_->ModelIsUnsat()) return;
   if (HasEnforcementLiteral(*ct)) return;
-  const Domain domain = ReadDomainFromProto(ct->linear());
+  const Domain rhs = ReadDomainFromProto(ct->linear());
 
   const LinearConstraintProto& arg = ct->linear();
   const int num_vars = arg.vars_size();
@@ -1429,14 +1440,13 @@ void CpModelPresolver::ExtractAtMostOneFromLinear(ConstraintProto* ct) {
       if (context_->MaxOf(ref) != 1) continue;
 
       if (type == 0) {
-        // TODO(user): we could potentially add one more Boolean with a lower
-        // coeff as long as we have lower_coeff + min_of_other_coeff >
-        // domain.Max().
-        if (min_sum + 2 * std::abs(coeff) > domain.Max()) {
+        // TODO(user): we could add one more Boolean with a lower coeff as long
+        // as we have lower_coeff + min_of_other_coeff > rhs.Max().
+        if (min_sum + 2 * std::abs(coeff) > rhs.Max()) {
           at_most_one.push_back(coeff > 0 ? ref : NegatedRef(ref));
         }
       } else {
-        if (max_sum - 2 * std::abs(coeff) < domain.Min()) {
+        if (max_sum - 2 * std::abs(coeff) < rhs.Min()) {
           at_most_one.push_back(coeff > 0 ? NegatedRef(ref) : ref);
         }
       }
@@ -1458,7 +1468,10 @@ void CpModelPresolver::ExtractAtMostOneFromLinear(ConstraintProto* ct) {
 // Convert some linear constraint involving only Booleans to their Boolean
 // form.
 bool CpModelPresolver::PresolveLinearOnBooleans(ConstraintProto* ct) {
-  if (context_->ModelIsUnsat()) return false;
+  if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear ||
+      context_->ModelIsUnsat()) {
+    return false;
+  }
 
   // TODO(user): the alternative to mark any newly created constraints might
   // be better.
@@ -1752,7 +1765,7 @@ bool CpModelPresolver::PresolveElement(ConstraintProto* ct) {
           reduced_index_domain = true;
         } else {
           ++num_vars;
-          if (domain.Min() == domain.Max()) {
+          if (domain.IsFixed()) {
             constant_set.insert(domain.Min());
           } else {
             all_constants = false;
@@ -2140,20 +2153,6 @@ bool CpModelPresolver::IntervalsCanIntersect(
 
 namespace {
 
-void MaybeDivideByGcd(std::map<int, int64>* objective_map, int64* divisor) {
-  if (objective_map->empty()) return;
-  int64 gcd = 0;
-  for (const auto entry : *objective_map) {
-    gcd = MathUtil::GCD64(gcd, std::abs(entry.second));
-    if (gcd == 1) break;
-  }
-  if (gcd == 1) return;
-  for (auto& entry : *objective_map) {
-    entry.second /= gcd;
-  }
-  *divisor *= gcd;
-}
-
 // Returns the sorted list of literals for given bool_or or at_most_one
 // constraint.
 std::vector<int> GetLiteralsFromSetPPCConstraint(ConstraintProto* ct) {
@@ -3014,73 +3013,17 @@ void CpModelPresolver::PresolvePureSatPart() {
 void CpModelPresolver::ExpandObjective() {
   if (context_->ModelIsUnsat()) return;
 
-  // Start by simplifying the objective domain using the implied domain of the
-  // initial linear objective.
-  {
-    Domain implied_domain(0);
-    for (int i = 0; i < context_->working_model->objective().vars_size(); ++i) {
-      const int ref = context_->working_model->objective().vars(i);
-      const int64 coeff = context_->working_model->objective().coeffs(i);
-      // NOTE: Domain multiplication with zero coeff returns empty domain so we
-      // avoid it by skipping the terms with zero coeff.
-      if (coeff == 0) continue;
-      implied_domain =
-          implied_domain
-              .AdditionWith(context_->DomainOf(ref).MultiplicationBy(coeff))
-              .RelaxIfTooComplex();
-    }
-
-    CpObjectiveProto* const mutable_objective =
-        context_->working_model->mutable_objective();
-    Domain old_domain = Domain::AllValues();
-    if (!mutable_objective->domain().empty()) {
-      old_domain = ReadDomainFromProto(*mutable_objective);
-    }
-    const Domain simplified_domain =
-        old_domain.SimplifyUsingImpliedDomain(implied_domain);
-    if (simplified_domain.IsEmpty()) {
-      return (void)context_->NotifyThatModelIsUnsat();
-    }
-    FillDomainInProto(simplified_domain, mutable_objective);
-  }
-
-  // Convert the objective linear expression to a map for ease of use below.
-  // We also only use affine representative.
-  std::map<int, int64> objective_map;
-  int64 objective_offset_change = 0;
-  int64 objective_divisor = 1;
-  for (int i = 0; i < context_->working_model->objective().vars_size(); ++i) {
-    const int ref = context_->working_model->objective().vars(i);
-    int64 coeff = context_->working_model->objective().coeffs(i);
-    if (!RefIsPositive(ref)) coeff = -coeff;
-    int var = PositiveRef(ref);
-
-    // Will be added back at the end.
-    context_->var_to_constraints[var].erase(-1);
-
-    if (context_->IsFixed(var)) {
-      objective_offset_change += coeff * context_->MinOf(var);
-      continue;
-    }
-
-    const AffineRelation::Relation r = context_->GetAffineRelation(var);
-    if (r.representative != var) {
-      var = r.representative;
-      objective_offset_change += r.offset * coeff;
-      coeff *= r.coeff;
-    }
-
-    objective_map[var] += coeff;
-    if (objective_map[var] == 0) objective_map.erase(var);
-  }
-  MaybeDivideByGcd(&objective_map, &objective_divisor);
+  // The objective is already loaded in the constext, but we re-canonicalize
+  // it with the latest information.
+  context_->CanonicalizeObjective();
 
   // If the objective is a single variable, then we can usually remove this
   // variable if it is only used in one linear equality constraint and we do
   // just one expansion. This is because the domain of the variable will be
   // transferred to our objective_domain.
   int unique_expanded_constraint = -1;
-  const bool objective_was_a_single_variable = objective_map.size() == 1;
+  const bool objective_was_a_single_variable =
+      context_->ObjectiveMap().size() == 1;
 
   // To avoid a bad complexity, we need to compute the number of relevant
   // constraints for each variables.
@@ -3107,7 +3050,7 @@ void CpModelPresolver::ExpandObjective() {
   }
 
   std::set<int> var_to_process;
-  for (const auto entry : objective_map) {
+  for (const auto entry : context_->ObjectiveMap()) {
     const int var = entry.first;
     CHECK(RefIsPositive(var));
     if (var_to_num_relevant_constraints[var] != 0) {
@@ -3118,17 +3061,23 @@ void CpModelPresolver::ExpandObjective() {
   // We currently never expand a variable more than once.
   int num_expansions = 0;
   absl::flat_hash_set<int> processed_vars;
+  std::vector<int> new_vars_in_objective;
   while (!relevant_constraints.empty()) {
     // Find a not yet expanded var.
     int objective_var = -1;
     while (!var_to_process.empty()) {
       const int var = *var_to_process.begin();
-      CHECK(!processed_vars.count(var));
+      CHECK(!processed_vars.contains(var));
       if (var_to_num_relevant_constraints[var] == 0) {
         processed_vars.insert(var);
         var_to_process.erase(var);
         continue;
       }
+      if (!context_->ObjectiveMap().contains(var)) {
+        // We do not mark it as processed in case it re-appear later.
+        var_to_process.erase(var);
+        continue;
+      }
       objective_var = var;
       break;
     }
@@ -3173,7 +3122,7 @@ void CpModelPresolver::ExpandObjective() {
           objective_coeff = (ref == objective_var) ? coeff : -coeff;
         } else {
           // This is not possible since we only consider relevant constraints.
-          CHECK(!processed_vars.count(PositiveRef(ref)));
+          CHECK(!processed_vars.contains(PositiveRef(ref)));
         }
       }
       CHECK(is_present);
@@ -3197,37 +3146,17 @@ void CpModelPresolver::ExpandObjective() {
       // Update the objective map. Note that the division is possible because
       // currently we only expand with coeff with a magnitude of 1.
       CHECK_EQ(std::abs(objective_coeff_in_expanded_constraint), 1);
-      const int64 factor =
-          objective_map[objective_var] / objective_coeff_in_expanded_constraint;
-
-      objective_map.erase(objective_var);
-
       const ConstraintProto& ct =
           context_->working_model->constraints(expanded_linear_index);
-      const int num_terms = ct.linear().vars_size();
-      for (int i = 0; i < num_terms; ++i) {
-        const int ref = ct.linear().vars(i);
-        const int var = PositiveRef(ref);
-        if (var == objective_var) continue;
+      context_->SubstituteVariableInObjective(
+          objective_var, objective_coeff_in_expanded_constraint, ct,
+          &new_vars_in_objective);
 
-        int64 coeff = -ct.linear().coeffs(i) * factor;
-        if (!RefIsPositive(ref)) coeff = -coeff;
-        if (!gtl::ContainsKey(objective_map, var)) {
-          if (!gtl::ContainsKey(processed_vars, var)) {
-            var_to_process.insert(var);
-          }
-        }
-        objective_map[var] += coeff;
-        if (objective_map[var] == 0.0) {
-          objective_map.erase(var);
-          var_to_process.erase(var);
-        }
+      // Add not yet processed new variables.
+      for (const int var : new_vars_in_objective) {
+        if (!processed_vars.contains(var)) var_to_process.insert(var);
       }
 
-      objective_offset_change +=
-          ct.linear().domain(0) * factor * objective_divisor;
-      MaybeDivideByGcd(&objective_map, &objective_divisor);
-
       // If the objective variable wasn't used in other constraints and it can
       // be reconstructed whatever the value of the other variables, we can
       // remove the constraint.
@@ -3237,7 +3166,7 @@ void CpModelPresolver::ExpandObjective() {
       if (context_->var_to_constraints[objective_var].size() == 1) {
         // Compute implied domain on objective_var.
         Domain implied_domain = ReadDomainFromProto(ct.linear());
-        for (int i = 0; i < num_terms; ++i) {
+        for (int i = 0; i < size_of_expanded_constraint; ++i) {
           const int ref = ct.linear().vars(i);
           if (PositiveRef(ref) == objective_var) continue;
           implied_domain =
@@ -3279,51 +3208,9 @@ void CpModelPresolver::ExpandObjective() {
     context_->UpdateConstraintVariableUsage(unique_expanded_constraint);
   }
 
-  // Compute the implied domain from the new objective linear expression.
-  Domain implied_domain(0);
-  for (const auto& entry : objective_map) {
-    implied_domain =
-        implied_domain
-            .AdditionWith(
-                context_->DomainOf(entry.first).MultiplicationBy(entry.second))
-            .RelaxIfTooComplex();
-  }
-
-  // Re-write the objective.
-  CpObjectiveProto* const mutable_objective =
-      context_->working_model->mutable_objective();
-  mutable_objective->clear_coeffs();
-  mutable_objective->clear_vars();
-  for (const auto& entry : objective_map) {
-    context_->var_to_constraints[PositiveRef(entry.first)].insert(-1);
-    mutable_objective->add_vars(entry.first);
-    mutable_objective->add_coeffs(entry.second);
-  }
-  Domain old_domain = Domain::AllValues();
-  if (!mutable_objective->domain().empty()) {
-    // Tricky, the domain in the proto do not include the offset.
-    old_domain = ReadDomainFromProto(*mutable_objective)
-                     .AdditionWith(Domain(-objective_offset_change))
-                     .InverseMultiplicationBy(objective_divisor);
-  }
-  const Domain simplified_domain =
-      old_domain.SimplifyUsingImpliedDomain(implied_domain);
-  if (simplified_domain.IsEmpty()) {
-    return (void)context_->NotifyThatModelIsUnsat();
-  }
-  FillDomainInProto(simplified_domain, mutable_objective);
-  mutable_objective->set_offset(mutable_objective->offset() +
-                                objective_offset_change);
-  if (objective_divisor > 1) {
-    const double divisor = static_cast<double>(objective_divisor);
-    mutable_objective->set_offset(mutable_objective->offset() / divisor);
-    if (mutable_objective->scaling_factor() == 0) {
-      mutable_objective->set_scaling_factor(divisor);
-    } else {
-      mutable_objective->set_scaling_factor(
-          mutable_objective->scaling_factor() * divisor);
-    }
-  }
+  // We re-do a canonicalization with the final linear expression.
+  context_->CanonicalizeObjective();
+  context_->WriteObjectiveToProto();
 }
 
 void CpModelPresolver::MergeNoOverlapConstraints() {
@@ -3514,20 +3401,24 @@ bool CpModelPresolver::PresolveOneConstraint(int c) {
       if (CanonicalizeLinear(ct)) {
         context_->UpdateConstraintVariableUsage(c);
       }
-      if (ct->constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
-        if (RemoveSingletonInLinear(ct)) {
+      if (PresolveSmallLinear(ct)) {
+        context_->UpdateConstraintVariableUsage(c);
+      }
+      if (PropagateDomainsInLinear(ct)) {
+        context_->UpdateConstraintVariableUsage(c);
+      }
+      // We first propagate the domains before calling this presolve rule.
+      if (RemoveSingletonInLinear(ct)) {
+        context_->UpdateConstraintVariableUsage(c);
+
+        // There is no need to re-do a propagation here, but the constraint
+        // size might have been reduced.
+        if (PresolveSmallLinear(ct)) {
           context_->UpdateConstraintVariableUsage(c);
         }
       }
-      if (ct->constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
-        if (PresolveLinear(ct)) {
-          context_->UpdateConstraintVariableUsage(c);
-        }
-      }
-      if (ct->constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
-        if (PresolveLinearOnBooleans(ct)) {
-          context_->UpdateConstraintVariableUsage(c);
-        }
+      if (PresolveLinearOnBooleans(ct)) {
+        context_->UpdateConstraintVariableUsage(c);
       }
       if (ct->constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
         const int old_num_enforcement_literals = ct->enforcement_literal_size();
@@ -3537,10 +3428,9 @@ bool CpModelPresolver::PresolveOneConstraint(int c) {
           context_->UpdateConstraintVariableUsage(c);
           return true;
         }
-        if (ct->enforcement_literal_size() > old_num_enforcement_literals) {
-          if (PresolveLinear(ct)) {
-            context_->UpdateConstraintVariableUsage(c);
-          }
+        if (ct->enforcement_literal_size() > old_num_enforcement_literals &&
+            PresolveSmallLinear(ct)) {
+          context_->UpdateConstraintVariableUsage(c);
         }
       }
       return false;
@@ -3759,88 +3649,86 @@ bool CpModelPresolver::ProcessSetPPC() {
   return changed;
 }
 
-void CpModelPresolver::TryToSimplifyDomains() {
+void CpModelPresolver::TryToSimplifyDomain(int var) {
   if (context_->ModelIsUnsat()) return;
+  if (context_->IsFixed(var)) return;
 
-  const int num_vars = context_->working_model->variables_size();
-  for (int var = 0; var < num_vars; ++var) {
-    if (context_->IsFixed(var)) continue;
+  const AffineRelation::Relation r = context_->GetAffineRelation(var);
+  if (r.representative != var) return;
 
-    const AffineRelation::Relation r = context_->GetAffineRelation(var);
-    if (r.representative != var) continue;
+  // Only process discrete domain.
+  const Domain& domain = context_->DomainOf(var);
 
-    // Only process discrete domain.
-    const Domain& domain = context_->DomainOf(var);
-
-    if (domain.Size() == 2 && domain.NumIntervals() == 1 && domain.Min() != 0) {
-      // Shifted Boolean variable.
-      const int new_var_index = context_->NewBoolVar();
-      const int64 offset = domain.Min();
-
-      ConstraintProto* const ct = context_->working_model->add_constraints();
-      LinearConstraintProto* const lin = ct->mutable_linear();
-      lin->add_vars(var);
-      lin->add_coeffs(1);
-      lin->add_vars(new_var_index);
-      lin->add_coeffs(-1);
-      lin->add_domain(offset);
-      lin->add_domain(offset);
-      context_->StoreAffineRelation(*ct, var, new_var_index, 1, offset);
-      context_->UpdateRuleStats("variables: canonicalize size two domain");
-      continue;
-    }
-
-    if (domain.NumIntervals() != domain.Size()) continue;
-
-    const int64 var_min = domain.Min();
-    int64 gcd = domain[1].start - var_min;
-    for (int index = 2; index < domain.NumIntervals(); ++index) {
-      const ClosedInterval& i = domain[index];
-      CHECK_EQ(i.start, i.end);
-      const int64 shifted_value = i.start - var_min;
-      CHECK_GE(shifted_value, 0);
-
-      gcd = MathUtil::GCD64(gcd, shifted_value);
-      if (gcd == 1) break;
-    }
-    if (gcd == 1) continue;
-
-    const int new_var_index = context_->working_model->variables_size();
-    IntegerVariableProto* const var_proto =
-        context_->working_model->add_variables();
-    {
-      std::vector<int64> scaled_values;
-      for (int index = 0; index < domain.NumIntervals(); ++index) {
-        const ClosedInterval& i = domain[index];
-        CHECK_EQ(i.start, i.end);
-        const int64 shifted_value = i.start - var_min;
-        scaled_values.push_back(shifted_value / gcd);
-      }
-      const Domain scaled_domain = Domain::FromValues(scaled_values);
-      FillDomainInProto(scaled_domain, var_proto);
-    }
-    context_->InitializeNewDomains();
-    if (context_->ModelIsUnsat()) return;
+  if (domain.Size() == 2 && domain.NumIntervals() == 1 && domain.Min() != 0) {
+    // Shifted Boolean variable.
+    const int new_var_index = context_->NewBoolVar();
+    const int64 offset = domain.Min();
 
     ConstraintProto* const ct = context_->working_model->add_constraints();
     LinearConstraintProto* const lin = ct->mutable_linear();
     lin->add_vars(var);
     lin->add_coeffs(1);
     lin->add_vars(new_var_index);
-    lin->add_coeffs(-gcd);
-    lin->add_domain(var_min);
-    lin->add_domain(var_min);
-
-    context_->StoreAffineRelation(*ct, var, new_var_index, gcd, var_min);
-    context_->UpdateRuleStats("variables: canonicalize affine domain");
+    lin->add_coeffs(-1);
+    lin->add_domain(offset);
+    lin->add_domain(offset);
+    context_->StoreAffineRelation(*ct, var, new_var_index, 1, offset);
+    context_->UpdateRuleStats("variables: canonicalize size two domain");
+    context_->UpdateNewConstraintsVariableUsage();
+    return;
   }
 
+  if (domain.NumIntervals() != domain.Size()) return;
+
+  const int64 var_min = domain.Min();
+  int64 gcd = domain[1].start - var_min;
+  for (int index = 2; index < domain.NumIntervals(); ++index) {
+    const ClosedInterval& i = domain[index];
+    CHECK_EQ(i.start, i.end);
+    const int64 shifted_value = i.start - var_min;
+    CHECK_GE(shifted_value, 0);
+
+    gcd = MathUtil::GCD64(gcd, shifted_value);
+    if (gcd == 1) break;
+  }
+  if (gcd == 1) return;
+
+  int new_var_index;
+  {
+    std::vector<int64> scaled_values;
+    for (int index = 0; index < domain.NumIntervals(); ++index) {
+      const ClosedInterval& i = domain[index];
+      CHECK_EQ(i.start, i.end);
+      const int64 shifted_value = i.start - var_min;
+      scaled_values.push_back(shifted_value / gcd);
+    }
+    new_var_index = context_->NewIntVar(Domain::FromValues(scaled_values));
+  }
+  if (context_->ModelIsUnsat()) return;
+
+  ConstraintProto* const ct = context_->working_model->add_constraints();
+  LinearConstraintProto* const lin = ct->mutable_linear();
+  lin->add_vars(var);
+  lin->add_coeffs(1);
+  lin->add_vars(new_var_index);
+  lin->add_coeffs(-gcd);
+  lin->add_domain(var_min);
+  lin->add_domain(var_min);
+
+  context_->StoreAffineRelation(*ct, var, new_var_index, gcd, var_min);
+  context_->UpdateRuleStats("variables: canonicalize affine domain");
   context_->UpdateNewConstraintsVariableUsage();
 }
 
 void CpModelPresolver::PresolveToFixPoint() {
   if (context_->ModelIsUnsat()) return;
 
+  // Do one pass to try to simplify the domain of variables.
+  const int num_vars = context_->working_model->variables_size();
+  for (int var = 0; var < num_vars; ++var) {
+    TryToSimplifyDomain(var);
+  }
+
   // This is used for constraint having unique variables in them (i.e. not
   // appearing anywhere else) to not call the presolve more than once for this
   // reason.
@@ -3896,19 +3784,25 @@ void CpModelPresolver::PresolveToFixPoint() {
       }
     }
 
-    // Look at variables to see if we can canonicalize the domain.
-    // Note that all the new constraint are just affine constraint and do
-    // not need to be presolved at this time.
-    TryToSimplifyDomains();
-    in_queue.resize(context_->working_model->constraints_size(), false);
-
     // Re-add to the queue the constraints that touch a variable that changed.
     //
     // TODO(user): Avoid reprocessing the constraints that changed the variables
     // with the use of timestamp.
     if (context_->ModelIsUnsat()) return;
     for (const int v : context_->modified_domains.PositionsSetAtLeastOnce()) {
-      if (context_->IsFixed(v)) context_->ExploitFixedDomain(v);
+      if (context_->IsFixed(v)) {
+        context_->ExploitFixedDomain(v);
+      } else {
+        // The domain changed, maybe we can canonicalize it.
+        //
+        // Important: This code is a bit brittle, because it assumes
+        // PositionsSetAtLeastOnce() will not change behind our back. That
+        // should however be the case because TryToSimplifyDomain() will only
+        // mark as modified via AddAffineRelation a variable that is already
+        // present in the modified set.
+        TryToSimplifyDomain(v);
+        in_queue.resize(context_->working_model->constraints_size(), false);
+      }
       for (const int c : context_->var_to_constraints[v]) {
         if (c >= 0 && !in_queue[c]) {
           in_queue[c] = true;
@@ -4119,13 +4013,9 @@ CpModelPresolver::CpModelPresolver(const PresolveOptions& options,
       context_->working_model->variables_size());
   context_->UpdateNewConstraintsVariableUsage();
 
-  // Hack for the objective variable(s) so that they are never considered to
-  // appear in only one constraint.
-  if (context_->working_model->has_objective()) {
-    for (const int obj_var : context_->working_model->objective().vars()) {
-      context_->var_to_constraints[PositiveRef(obj_var)].insert(-1);
-    }
-  }
+  // Initialize the objective.
+  context_->ReadObjectiveFromProto();
+  context_->CanonicalizeObjective();
 }
 
 // The presolve works as follow:
diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h
index ed6ae00794..b9c82229a8 100644
--- a/ortools/sat/cp_model_presolve.h
+++ b/ortools/sat/cp_model_presolve.h
@@ -102,10 +102,6 @@ class CpModelPresolver {
   bool PresolveTable(ConstraintProto* ct);
   bool PresolveElement(ConstraintProto* ct);
   bool PresolveInterval(int c, ConstraintProto* ct);
-  bool PresolveLinear(ConstraintProto* ct);
-  bool PresolveLinearOnBooleans(ConstraintProto* ct);
-  bool CanonicalizeLinear(ConstraintProto* ct);
-  bool RemoveSingletonInLinear(ConstraintProto* ct);
   bool PresolveIntDiv(ConstraintProto* ct);
   bool PresolveIntProd(ConstraintProto* ct);
   bool PresolveIntMin(ConstraintProto* ct);
@@ -116,6 +112,13 @@ class CpModelPresolver {
   bool PresolveBoolOr(ConstraintProto* ct);
   bool PresolveEnforcementLiteral(ConstraintProto* ct);
 
+  // For the linear constraints, we have more than one function.
+  bool CanonicalizeLinear(ConstraintProto* ct);
+  bool PropagateDomainsInLinear(ConstraintProto* ct);
+  bool RemoveSingletonInLinear(ConstraintProto* ct);
+  bool PresolveSmallLinear(ConstraintProto* ct);
+  bool PresolveLinearOnBooleans(ConstraintProto* ct);
+
   // SetPPC is short for set packing, partitioning and covering constraints.
   // These are sum of booleans <=, = and >= 1 respectively.
   bool ProcessSetPPC();
@@ -144,7 +147,7 @@ class CpModelPresolver {
 
   void ExpandObjective();
 
-  void TryToSimplifyDomains();
+  void TryToSimplifyDomain(int var);
 
   void MergeNoOverlapConstraints();
 
@@ -161,6 +164,9 @@ class CpModelPresolver {
   const PresolveOptions& options_;
   std::vector<int>* postsolve_mapping_;
   PresolveContext* context_;
+
+  // Used by CanonicalizeLinear().
+  std::vector<std::pair<int, int64>> tmp_terms_;
 };
 
 // Convenient wrapper to call the full presolve.
diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc
index 6f8c126d5b..77ed1073cf 100644
--- a/ortools/sat/linear_constraint_manager.cc
+++ b/ortools/sat/linear_constraint_manager.cc
@@ -163,8 +163,17 @@ void LinearConstraintManager::ComputeObjectiveParallelism(
     constraint_infos_[ct_index].objective_parallelism = 0.0;
     return;
   }
+
+  const LinearConstraint& lc = constraint_infos_[ct_index].constraint;
+  double unscaled_objective_parallelism = 0.0;
+  for (int i = 0; i < lc.vars.size(); ++i) {
+    if (lc.vars[i] < dense_objective_coeffs_.size()) {
+      unscaled_objective_parallelism +=
+          ToDouble(lc.coeffs[i]) * dense_objective_coeffs_[lc.vars[i]];
+    }
+  }
   const double objective_parallelism =
-      ScalarProduct(constraint_infos_[ct_index].constraint, objective_) /
+      unscaled_objective_parallelism /
       (constraint_infos_[ct_index].l2_norm * objective_l2_norm_);
   constraint_infos_[ct_index].objective_parallelism =
       std::abs(objective_parallelism);
@@ -335,6 +344,8 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
 
 bool LinearConstraintManager::ChangeLp(
     const gtl::ITIVector<IntegerVariable, double>& lp_solution) {
+  VLOG(3) << "Enter ChangeLP, scan " << constraint_infos_.size()
+          << " constraints";
   std::vector<ConstraintIndex> new_constraints;
   std::vector<double> new_constraints_efficacies;
   std::vector<double> new_constraints_orthogonalities;
@@ -346,6 +357,7 @@ bool LinearConstraintManager::ChangeLp(
   // We keep any constraints that is already present, and otherwise, we add the
   // ones that are currently not satisfied by at least "tolerance".
   const double tolerance = 1e-6;
+  FillDenseObjectiveCoeffs();
   for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
     if (constraint_infos_[i].permanently_removed) continue;
 
@@ -387,6 +399,17 @@ bool LinearConstraintManager::ChangeLp(
       new_constraints_efficacies.push_back(violation /
                                            constraint_infos_[i].l2_norm);
       new_constraints_orthogonalities.push_back(1.0);
+
+      if (objective_is_defined_ &&
+          !constraint_infos_[i].objective_parallelism_computed) {
+        ComputeObjectiveParallelism(i);
+      } else if (!objective_is_defined_) {
+        constraint_infos_[i].objective_parallelism = 0.0;
+      }
+
+      constraint_infos_[i].current_score =
+          new_constraints_efficacies.back() +
+          constraint_infos_[i].objective_parallelism;
     }
   }
 
@@ -396,18 +419,36 @@ bool LinearConstraintManager::ChangeLp(
     RemoveMarkedConstraints();
   }
 
-  // Note that the algo below is in O(limit * new_constraint), so this limit
-  // should stay low.
+  // Note that the algo below is in O(limit * new_constraint). In order to
+  // limit spending too much time on this, we first sort all the constraints
+  // with an imprecise score (no orthogonality), then limit the size of the
+  // vector of constraints to precisely score, then we do the actual scoring.
+  //
+  // On problem crossword_opt_grid-19.05_dict-80_sat with linearization_level=2,
+  // new_constraint.size() > 1.5M.
+  //
+  // TODO(user): This blowup factor could be adaptative w.r.t. the constraint
+  // limit.
+  const int kBlowupFactor = 4;
   int constraint_limit = std::min(50, static_cast<int>(new_constraints.size()));
   if (lp_constraints_.empty()) {
     constraint_limit = std::min(1000, static_cast<int>(new_constraints.size()));
   }
+  VLOG(3) << "   - size = " << new_constraints.size()
+          << ", limit = " << constraint_limit;
 
-  // TODO(user,user): on problem crossword_opt_grid-19.05_dict-80_sat with
-  //     linearization_level=2, new_constraint.size() > 1.5M. Improve
-  //     complexity of the following loop.
+  std::stable_sort(new_constraints.begin(), new_constraints.end(),
+                   [&](ConstraintIndex a, ConstraintIndex b) {
+                     return constraint_infos_[a].current_score >
+                            constraint_infos_[b].current_score;
+                   });
+  if (new_constraints.size() > kBlowupFactor * constraint_limit) {
+    VLOG(3) << "Resize candidate constraints from " << new_constraints.size()
+            << " down to " << kBlowupFactor * constraint_limit;
+    new_constraints.resize(kBlowupFactor * constraint_limit);
+  }
 
-  int skipped_checks = 0;
+  int num_skipped_checks = 0;
   const int kCheckFrequency = 100;
   ConstraintIndex last_added_candidate = kInvalidConstraintIndex;
   for (int i = 0; i < constraint_limit; ++i) {
@@ -416,10 +457,10 @@ bool LinearConstraintManager::ChangeLp(
     double best_score = 0.0;
     ConstraintIndex best_candidate = kInvalidConstraintIndex;
     for (int j = 0; j < new_constraints.size(); ++j) {
-      // Checks the time limit, and returns not_changed if crossed.
-      if (++skipped_checks >= kCheckFrequency) {
-        if (time_limit_->LimitReached()) return false;
-        skipped_checks = 0;
+      // Checks the time limit, and returns if the lp has changed.
+      if (++num_skipped_checks >= kCheckFrequency) {
+        if (time_limit_->LimitReached()) return current_lp_is_changed_;
+        num_skipped_checks = 0;
       }
 
       const ConstraintIndex new_constraint = new_constraints[j];
@@ -449,16 +490,10 @@ bool LinearConstraintManager::ChangeLp(
         continue;
       }
 
-      if (objective_is_defined_ &&
-          !constraint_infos_[new_constraint].objective_parallelism_computed) {
-        ComputeObjectiveParallelism(new_constraint);
-      }
-
       // TODO(user): Experiment with different weights or different
       // functions for computing score.
-      const double score =
-          new_constraints_orthogonalities[j] + new_constraints_efficacies[j] +
-          constraint_infos_[new_constraint].objective_parallelism;
+      const double score = new_constraints_orthogonalities[j] +
+                           constraint_infos_[new_constraint].current_score;
       CHECK_GE(score, 0.0);
       if (score > best_score || best_candidate == kInvalidConstraintIndex) {
         best_score = score;
@@ -478,6 +513,7 @@ bool LinearConstraintManager::ChangeLp(
     }
   }
 
+  VLOG(3) << "   - Exit ChangeLP";
   // The LP changed only if we added new constraints or if some constraints
   // already inside changed (simplification or tighter bounds).
   if (current_lp_is_changed_) {
@@ -495,5 +531,16 @@ void LinearConstraintManager::AddAllConstraintsToLp() {
   }
 }
 
+void LinearConstraintManager::FillDenseObjectiveCoeffs() {
+  if (objective_.vars.empty()) return;
+  DCHECK(std::is_sorted(objective_.vars.begin(), objective_.vars.end()));
+  const IntegerVariable last_var = objective_.vars.back();
+  dense_objective_coeffs_.assign(last_var.value() + 1, 0.0);
+  for (int i = 0; i < objective_.vars.size(); ++i) {
+    dense_objective_coeffs_[objective_.vars[i]] =
+        ToDouble(objective_.coeffs[i]);
+  }
+}
+
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h
index c126dbbe9b..0724116eee 100644
--- a/ortools/sat/linear_constraint_manager.h
+++ b/ortools/sat/linear_constraint_manager.h
@@ -53,6 +53,7 @@ class LinearConstraintManager {
     // parallel to one of the existing constraints in the LP.
     bool permanently_removed;
     size_t hash;
+    double current_score;
   };
 
   explicit LinearConstraintManager(Model* model)
@@ -114,6 +115,9 @@ class LinearConstraintManager {
   // Returns true if the terms of the constraint changed.
   bool SimplifyConstraint(LinearConstraint* ct);
 
+  // Helper method to fill in the objective_values_ vector.
+  void FillDenseObjectiveCoeffs();
+
   // Helper method to compute objective parallelism for a given constraint. This
   // also lazily computes objective norm.
   void ComputeObjectiveParallelism(const ConstraintIndex ct_index);
@@ -155,6 +159,9 @@ class LinearConstraintManager {
   bool objective_norm_computed_ = false;
   LinearConstraint objective_;
   double objective_l2_norm_ = 0.0;
+  // Dense representation of the objective coeffs indexed by variables indices.
+  // It contains 0.0 where the variables does not appear in the objective.
+  gtl::ITIVector<IntegerVariable, double> dense_objective_coeffs_;
 
   TimeLimit* time_limit_;
 };
diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc
index 7363dac559..a8d53a79a5 100644
--- a/ortools/sat/optimization.cc
+++ b/ortools/sat/optimization.cc
@@ -1243,7 +1243,10 @@ SatSolver::Status FindCores(std::vector<Literal> assumptions,
                             std::vector<std::vector<Literal>>* cores) {
   cores->clear();
   SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
+  TimeLimit* limit = model->GetOrCreate<TimeLimit>();
   do {
+    if (limit->LimitReached()) return SatSolver::LIMIT_REACHED;
+
     const SatSolver::Status result =
         ResetAndSolveIntegerProblem(assumptions, model);
     if (result != SatSolver::ASSUMPTIONS_UNSAT) return result;
diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc
index 11960c42cd..38d3a2787d 100644
--- a/ortools/sat/presolve_context.cc
+++ b/ortools/sat/presolve_context.cc
@@ -14,6 +14,7 @@
 #include "ortools/sat/presolve_context.h"
 
 #include "ortools/base/map_util.h"
+#include "ortools/base/mathutil.h"
 #include "ortools/port/proto_utils.h"
 
 namespace operations_research {
@@ -66,8 +67,8 @@ bool PresolveContext::DomainIsEmpty(int ref) const {
 }
 
 bool PresolveContext::IsFixed(int ref) const {
-  CHECK(!DomainIsEmpty(ref));
-  return domains[PositiveRef(ref)].Min() == domains[PositiveRef(ref)].Max();
+  DCHECK(!DomainIsEmpty(ref));
+  return domains[PositiveRef(ref)].IsFixed();
 }
 
 bool PresolveContext::LiteralIsTrue(int lit) const {
@@ -263,10 +264,13 @@ void PresolveContext::StoreAffineRelation(const ConstraintProto& ct, int ref_x,
     added |= var_equiv_relations.TryAdd(x, y, c, o, allow_rep_x, allow_rep_y);
   }
   if (added) {
-    // The domain didn't change, but this notification allows to re-process
-    // any constraint containing these variables.
-    modified_domains.Set(x);
-    modified_domains.Set(y);
+    // The domain didn't change, but this notification allows to re-process any
+    // constraint containing these variables. Note that we do not need to
+    // retrigger a propagation of the constraint containing a variable whose
+    // representative didn't change.
+    const int new_rep = affine_relations.Get(rep_x).representative;
+    if (new_rep != rep_x) modified_domains.Set(x);
+    if (new_rep != rep_y) modified_domains.Set(y);
     affine_constraints.insert(&ct);
   }
 }
@@ -277,20 +281,6 @@ void PresolveContext::StoreBooleanEqualityRelation(int ref_a, int ref_b) {
     is_unsat = true;
     return;
   }
-  bool added = false;
-  if (RefIsPositive(ref_a) == RefIsPositive(ref_b)) {
-    added |=
-        affine_relations.TryAdd(PositiveRef(ref_a), PositiveRef(ref_b), 1, 0);
-    added |= var_equiv_relations.TryAdd(PositiveRef(ref_a), PositiveRef(ref_b),
-                                        1, 0);
-  } else {
-    added |=
-        affine_relations.TryAdd(PositiveRef(ref_a), PositiveRef(ref_b), -1, 1);
-  }
-  if (!added) return;
-
-  modified_domains.Set(PositiveRef(ref_a));
-  modified_domains.Set(PositiveRef(ref_b));
 
   // For now, we do need to add the relation ref_a == ref_b so we have a
   // proper variable usage count and propagation between ref_a and ref_b.
@@ -308,14 +298,17 @@ void PresolveContext::StoreBooleanEqualityRelation(int ref_a, int ref_b) {
     arg->add_coeffs(-1);
     arg->add_domain(0);
     arg->add_domain(0);
+    StoreAffineRelation(*ct, PositiveRef(ref_a), PositiveRef(ref_b), 1,
+                        /*offset=*/0);
   } else {
     // a = 1 - b
     arg->add_coeffs(1);
     arg->add_coeffs(1);
     arg->add_domain(1);
     arg->add_domain(1);
+    StoreAffineRelation(*ct, PositiveRef(ref_a), PositiveRef(ref_b), -1,
+                        /*offset=*/1);
   }
-  affine_constraints.insert(ct);
   UpdateNewConstraintsVariableUsage();
 }
 
@@ -454,22 +447,194 @@ int PresolveContext::GetOrCreateVarValueEncoding(int ref, int64 value) {
   return literal;
 }
 
-void PresolveContext::SubstituteVariableInObjective(
-    int var, int64 coeff, const ConstraintProto& equality) {
-  // Remove objective entry from var_to_constraint lists. The objective
-  // entries will be added back after the substitution.
-  const CpObjectiveProto& objective = working_model->objective();
-  for (const int ref : objective.vars()) {
-    var_to_constraints[PositiveRef(ref)].erase(-1);
+void PresolveContext::ReadObjectiveFromProto() {
+  const CpObjectiveProto& obj = working_model->objective();
+
+  objective_offset = obj.offset();
+  objective_scaling_factor = obj.scaling_factor();
+  if (objective_scaling_factor == 0.0) {
+    objective_scaling_factor = 1.0;
+  }
+  if (!obj.domain().empty()) {
+    objective_domain = ReadDomainFromProto(obj);
+  } else {
+    objective_domain = Domain::AllValues();
   }
 
-  ::operations_research::sat::SubstituteVariableInObjective(
-      var, coeff, equality, working_model->mutable_objective());
+  objective_map.clear();
+  for (int i = 0; i < obj.vars_size(); ++i) {
+    const int ref = obj.vars(i);
+    int64 coeff = obj.coeffs(i);
+    if (!RefIsPositive(ref)) coeff = -coeff;
+    int var = PositiveRef(ref);
 
-  // Add back the objective entry for the variables in objective in
-  // var_to_constraint lists.
-  for (const int ref : objective.vars()) {
-    var_to_constraints[PositiveRef(ref)].insert(-1);
+    objective_map[var] += coeff;
+    if (objective_map[var] == 0) {
+      objective_map.erase(var);
+      var_to_constraints[var].erase(-1);
+    } else {
+      var_to_constraints[var].insert(-1);
+    }
+  }
+}
+
+void PresolveContext::CanonicalizeObjective() {
+  int64 offset_change = 0;
+
+  // We replace each entry by its affine representative.
+  // Note that the non-deterministic loop is fine.
+  //
+  // TODO(user): This is a bit dupplicated with the presolve linear code.
+  // We also do not propagate back any domain restriction from the objective to
+  // the variables if any.
+  for (const auto& entry : objective_map) {
+    const int var = entry.first;
+    const int64 coeff = entry.second;
+
+    if (IsFixed(var)) {
+      offset_change += coeff * MinOf(var);
+      var_to_constraints[var].erase(-1);
+      objective_map.erase(var);
+      continue;
+    }
+
+    const AffineRelation::Relation r = GetAffineRelation(var);
+    if (r.representative == var) continue;
+
+    objective_map.erase(var);
+    var_to_constraints[var].erase(-1);
+
+    // Do the substitution.
+    offset_change += coeff * r.offset;
+    const int64 new_coeff = objective_map[r.representative] += coeff * r.coeff;
+
+    // Process new term.
+    if (new_coeff == 0) {
+      objective_map.erase(r.representative);
+      var_to_constraints[r.representative].erase(-1);
+    } else {
+      var_to_constraints[r.representative].insert(-1);
+      if (IsFixed(r.representative)) {
+        offset_change += new_coeff * MinOf(r.representative);
+        var_to_constraints[r.representative].erase(-1);
+        objective_map.erase(r.representative);
+      }
+    }
+  }
+
+  Domain implied_domain(0);
+  int64 gcd(0);
+
+  // We need to sort the entries to be deterministic.
+  {
+    std::vector<std::pair<int, int64>> entries;
+    for (const auto& entry : objective_map) {
+      entries.push_back(entry);
+    }
+    std::sort(entries.begin(), entries.end());
+    for (const auto& entry : entries) {
+      const int var = entry.first;
+      const int64 coeff = entry.second;
+      gcd = MathUtil::GCD64(gcd, std::abs(coeff));
+      implied_domain =
+          implied_domain.AdditionWith(DomainOf(var).MultiplicationBy(coeff))
+              .RelaxIfTooComplex();
+    }
+  }
+
+  // This is the new domain.
+  // Note that the domain never include the offset.
+  objective_domain = objective_domain.AdditionWith(Domain(-offset_change))
+                         .IntersectionWith(implied_domain);
+  objective_domain =
+      objective_domain.SimplifyUsingImpliedDomain(implied_domain);
+
+  // Updat the offset.
+  objective_offset += offset_change;
+
+  // Maybe divide by GCD.
+  if (gcd > 1) {
+    for (auto& entry : objective_map) {
+      entry.second /= gcd;
+    }
+    objective_domain = objective_domain.InverseMultiplicationBy(gcd);
+    objective_offset /= static_cast<double>(gcd);
+    objective_scaling_factor *= static_cast<double>(gcd);
+  }
+}
+
+void PresolveContext::SubstituteVariableInObjective(
+    int var_in_equality, int64 coeff_in_equality,
+    const ConstraintProto& equality, std::vector<int>* new_vars_in_objective) {
+  CHECK(equality.enforcement_literal().empty());
+  CHECK(RefIsPositive(var_in_equality));
+  CHECK_EQ(std::abs(coeff_in_equality), 1);
+
+  if (new_vars_in_objective != nullptr) new_vars_in_objective->clear();
+
+  const int64 coeff_in_objective = gtl::FindOrDie(objective_map, var_in_equality);
+  for (int i = 0; i < equality.linear().vars().size(); ++i) {
+    int var = equality.linear().vars(i);
+    int64 coeff = equality.linear().coeffs(i);
+    if (!RefIsPositive(var)) {
+      var = NegatedRef(var);
+      coeff = -coeff;
+    }
+    if (var == var_in_equality) continue;
+
+    // We should divided by coeff_in_equality, but since its magnitude is one,
+    // multiplying is the same.
+    int64& map_ref = objective_map[var];
+    if (map_ref == 0 && new_vars_in_objective != nullptr) {
+      new_vars_in_objective->push_back(var);
+    }
+    map_ref -= coeff_in_objective * coeff * coeff_in_equality;
+
+    if (map_ref == 0) {
+      objective_map.erase(var);
+      var_to_constraints[var].erase(-1);
+    } else {
+      var_to_constraints[var].insert(-1);
+    }
+  }
+
+  objective_map.erase(var_in_equality);
+  var_to_constraints[var_in_equality].erase(-1);
+
+  // Deal with the offset.
+  Domain offset = ReadDomainFromProto(equality.linear());
+  DCHECK_EQ(offset.Min(), offset.Max());
+  bool exact = true;
+  offset =
+      offset.MultiplicationBy(coeff_in_objective * coeff_in_equality, &exact);
+  CHECK(exact);
+
+  // Tricky: The objective domain is without the offset, so we need to shift it
+  objective_offset += static_cast<double>(offset.Min());
+  objective_domain = objective_domain.AdditionWith(Domain(-offset.Min()));
+}
+
+void PresolveContext::WriteObjectiveToProto() {
+  if (objective_domain.IsEmpty()) {
+    return (void)NotifyThatModelIsUnsat();
+  }
+
+  // We need to sort the entries to be deterministic.
+  std::vector<std::pair<int, int64>> entries;
+  for (const auto& entry : objective_map) {
+    entries.push_back(entry);
+  }
+  std::sort(entries.begin(), entries.end());
+
+  CpObjectiveProto* mutable_obj = working_model->mutable_objective();
+  mutable_obj->set_offset(objective_offset);
+  mutable_obj->set_scaling_factor(objective_scaling_factor);
+  FillDomainInProto(objective_domain, mutable_obj);
+  mutable_obj->clear_vars();
+  mutable_obj->clear_coeffs();
+  for (const auto& entry : entries) {
+    mutable_obj->add_vars(entry.first);
+    mutable_obj->add_coeffs(entry.second);
   }
 }
 
diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h
index 459634fb3c..c2bb5946f7 100644
--- a/ortools/sat/presolve_context.h
+++ b/ortools/sat/presolve_context.h
@@ -133,19 +133,33 @@ struct PresolveContext {
   // create it, add the corresponding constraints and returns it.
   int GetOrCreateVarValueEncoding(int ref, int64 value);
 
+  // Objective handling functions. We load it at the beginning so that during
+  // presolve we can work on the more efficient hash_map representation.
+  //
+  // Note that ReadObjectiveFromProto() makes sure that var_to_constraints of
+  // all the variable that appear in the objective contains -1. This is later
+  // enforced by all the functions modifying the objective.
+  void ReadObjectiveFromProto();
+  void CanonicalizeObjective();
+  void WriteObjectiveToProto();
+
   // Given a variable defined by the given inequality that also appear in the
   // objective, remove it from the objective by transferring its cost to other
   // variables in the equality.
   //
-  // TODO(user): Each time this is called, we do a linear scan of the objective
-  // that can sometimes be really large (millions of variables). It will be more
-  // efficient to just store the objective in a map at the beginning of the
-  // presolve so this can be done in O(equality_size) and re-encode the
-  // objective in the proto at the end. This code already exist in
-  // ExpandObjective() and can be factored out. Look at the "ivu*.mps" problems
-  // where the presolve is slow because of that.
-  void SubstituteVariableInObjective(int var, int64 coeff,
-                                     const ConstraintProto& equality);
+  // If new_vars_in_objective is not nullptr, it will be filled with "new"
+  // variables that where not in the objective before and are after
+  // substitution.
+  void SubstituteVariableInObjective(
+      int var_in_equality, int64 coeff_in_equality,
+      const ConstraintProto& equality,
+      std::vector<int>* new_vars_in_objective = nullptr);
+
+  // Objective getters.
+  const Domain& ObjectiveDomain() const { return objective_domain; }
+  const absl::flat_hash_map<int, int64>& ObjectiveMap() const {
+    return objective_map;
+  }
 
   // This regroups all the affine relations between variables. Note that the
   // constraints used to detect such relations will not be removed from the
@@ -228,6 +242,15 @@ struct PresolveContext {
 
   // The current domain of each variables.
   std::vector<Domain> domains;
+
+  // Internal representation of the objective. During presolve, we first load
+  // the objective in this format in order to have more efficient substitution
+  // on large problems (also because the objective is often dense). At the end
+  // we re-convert it to its proto form.
+  absl::flat_hash_map<int, int64> objective_map;
+  Domain objective_domain;
+  double objective_offset;
+  double objective_scaling_factor;
 };
 
 }  // namespace sat
diff --git a/ortools/sat/presolve_util.cc b/ortools/sat/presolve_util.cc
index 3ae46dd45b..4c1069036f 100644
--- a/ortools/sat/presolve_util.cc
+++ b/ortools/sat/presolve_util.cc
@@ -206,46 +206,5 @@ void SubstituteVariable(int var, int64 var_coeff_in_definition,
   SortAndMergeTerms(&terms, ct->mutable_linear());
 }
 
-bool VarFoundInObjective(int var, CpObjectiveProto* obj) {
-  const int size = obj->vars_size();
-  for (int i = 0; i < size; ++i) {
-    if (PositiveRef(obj->vars(i)) == PositiveRef(var)) {
-      return true;
-    }
-  }
-  return false;
-}
-
-void SubstituteVariableInObjective(int var, int64 var_coeff_in_definition,
-                                   const ConstraintProto& definition,
-                                   CpObjectiveProto* obj) {
-  CHECK(RefIsPositive(var));
-  CHECK_EQ(std::abs(var_coeff_in_definition), 1);
-
-  // Copy all the terms (except the one refering to var).
-  std::vector<std::pair<int, int64>> terms;
-  int64 var_coeff_in_obj = GetVarCoeffAndCopyOtherTerms(var, *obj, &terms);
-
-  if (var_coeff_in_definition < 0) var_coeff_in_obj *= -1;
-
-  AddTermsFromVarDefinition(var, var_coeff_in_obj, definition, &terms);
-
-  bool exact = false;
-  Domain offset = ReadDomainFromProto(definition.linear());
-  DCHECK_EQ(offset.Min(), offset.Max());
-  offset = offset.MultiplicationBy(var_coeff_in_obj, &exact);
-  CHECK(exact);
-
-  // Tricky: The objective domain is without the offset, so we need to shift it.
-  obj->set_offset(offset.Min() + obj->offset());
-  if (!obj->domain().empty()) {
-    Domain old_domain = ReadDomainFromProto(*obj);
-    FillDomainInProto(old_domain.AdditionWith(Domain(-offset.Min())), obj);
-  }
-
-  // Sort and merge terms refering to the same variable.
-  SortAndMergeTerms(&terms, obj);
-}
-
 }  // namespace sat
 }  // namespace operations_research
diff --git a/ortools/sat/presolve_util.h b/ortools/sat/presolve_util.h
index 5d6f80263e..6f9d139a5a 100644
--- a/ortools/sat/presolve_util.h
+++ b/ortools/sat/presolve_util.h
@@ -88,15 +88,6 @@ class DomainDeductions {
 void SubstituteVariable(int var, int64 var_coeff_in_definition,
                         const ConstraintProto& definition, ConstraintProto* ct);
 
-// Returns true if the variable is found in objective linear terms.
-bool VarFoundInObjective(int var, CpObjectiveProto* obj);
-
-// Replaces the variable var in objective using the definition constraint.
-// Currently the coefficient in the definition must be 1 or -1.
-void SubstituteVariableInObjective(int var, int64 var_coeff_in_definition,
-                                   const ConstraintProto& definition,
-                                   CpObjectiveProto* obj);
-
 }  // namespace sat
 }  // namespace operations_research
 
diff --git a/ortools/util/sorted_interval_list.cc b/ortools/util/sorted_interval_list.cc
index 1e797256f3..54b063a703 100644
--- a/ortools/util/sorted_interval_list.cc
+++ b/ortools/util/sorted_interval_list.cc
@@ -167,6 +167,8 @@ Domain Domain::FromVectorIntervals(
 
 bool Domain::IsEmpty() const { return intervals_.empty(); }
 
+bool Domain::IsFixed() const { return Min() == Max(); }
+
 int64 Domain::Size() const {
   int64 size = 0;
   for (const ClosedInterval interval : intervals_) {
@@ -180,12 +182,12 @@ int64 Domain::Size() const {
 }
 
 int64 Domain::Min() const {
-  CHECK(!IsEmpty());
+  DCHECK(!IsEmpty());
   return intervals_.front().start;
 }
 
 int64 Domain::Max() const {
-  CHECK(!IsEmpty());
+  DCHECK(!IsEmpty());
   return intervals_.back().end;
 }
 
diff --git a/ortools/util/sorted_interval_list.h b/ortools/util/sorted_interval_list.h
index 8d2aee2036..3ac436e887 100644
--- a/ortools/util/sorted_interval_list.h
+++ b/ortools/util/sorted_interval_list.h
@@ -143,18 +143,22 @@ class Domain {
 
   /**
    * Returns the min value of the domain.
-   *
-   * This checks that the domain is not empty.
+   * The domain must not be empty.
    */
   int64 Min() const;
 
   /**
    * Returns the max value of the domain.
-   *
-   * This checks that the domain is not empty.
+   * The domain must not be empty.
    */
   int64 Max() const;
 
+  /**
+   * Returns true iff the domain is reduced to a single value.
+   * The domain must not be empty.
+   */
+  bool IsFixed() const;
+
   /**
    * Returns true iff value is in Domain.
    */