// Copyright 2010-2025 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "ortools/sat/util.h"

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <deque>
#include <limits>
#include <numeric>
#include <string>
#include <utility>
#include <vector>

#include "absl/algorithm/container.h"
#include "absl/container/btree_set.h"
#include "absl/log/check.h"
#include "absl/numeric/bits.h"
#include "absl/numeric/int128.h"
#include "absl/random/bit_gen_ref.h"
#include "absl/random/distributions.h"
#include "absl/strings/str_cat.h"
#include "absl/types/span.h"
#include "google/protobuf/descriptor.h"
#include "ortools/base/logging.h"
#include "ortools/base/mathutil.h"
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/util/saturated_arithmetic.h"
#include "ortools/util/sorted_interval_list.h"

namespace operations_research {
namespace sat {

std::string FormatCounter(int64_t num) {
  std::string s = absl::StrCat(num);
  std::string out;
  const int size = s.size();
  for (int i = 0; i < size; ++i) {
    if (i > 0 && (size - i) % 3 == 0) {
      out.push_back('\'');
    }
    out.push_back(s[i]);
  }
  return out;
}

namespace {

inline std::string LeftAlign(std::string s, int size = 16) {
  if (s.size() >= size) return s;
  s.resize(size, ' ');
  return s;
}

inline std::string RightAlign(std::string s, int size = 16) {
  if (s.size() >= size) return s;
  return absl::StrCat(std::string(size - s.size(), ' '), s);
}

}  // namespace

std::string FormatTable(std::vector<std::vector<std::string>>& table,
                        int spacing) {
  if (table.size() > 1) {
    // We order by name.
    std::sort(table.begin() + 1, table.end());
  }

  std::vector<int> widths;
  for (const std::vector<std::string>& line : table) {
    if (line.size() > widths.size()) widths.resize(line.size(), spacing);
    for (int j = 0; j < line.size(); ++j) {
      widths[j] = std::max<int>(widths[j], line[j].size() + spacing);
    }
  }
  std::string output;
  for (int i = 0; i < table.size(); ++i) {
    for (int j = 0; j < table[i].size(); ++j) {
      if (i == 0 && j == 0) {
        // We currently only left align the table name.
        absl::StrAppend(&output, LeftAlign(table[i][j], widths[j]));
      } else {
        absl::StrAppend(&output, RightAlign(table[i][j], widths[j]));
      }
    }
    absl::StrAppend(&output, "\n");
  }
  return output;
}

void RandomizeDecisionHeuristic(absl::BitGenRef random,
                                SatParameters* parameters) {
#if !defined(__PORTABLE_PLATFORM__)
  // Random preferred variable order.
  const google::protobuf::EnumDescriptor* order_d =
      SatParameters::VariableOrder_descriptor();
  parameters->set_preferred_variable_order(
      static_cast<SatParameters::VariableOrder>(
          order_d->value(absl::Uniform(random, 0, order_d->value_count()))
              ->number()));

  // Random polarity initial value.
  const google::protobuf::EnumDescriptor* polarity_d =
      SatParameters::Polarity_descriptor();
  parameters->set_initial_polarity(static_cast<SatParameters::Polarity>(
      polarity_d->value(absl::Uniform(random, 0, polarity_d->value_count()))
          ->number()));
#endif  // __PORTABLE_PLATFORM__
  // Other random parameters.
  parameters->set_use_phase_saving(absl::Bernoulli(random, 0.5));
  parameters->set_random_polarity_ratio(absl::Bernoulli(random, 0.5) ? 0.01
                                                                     : 0.0);
  parameters->set_random_branches_ratio(absl::Bernoulli(random, 0.5) ? 0.01
                                                                     : 0.0);
}

namespace {

// This will be optimized into one division. I tested that in
// `BM_DivisionAndRemainderAlternative` (sat/integer_test.cc)
//
// Note that I am not 100% sure we need the indirection for the optimization
// to kick in though, but this seemed safer given our weird r[i ^ 1] inputs.
void QuotientAndRemainder(int64_t a, int64_t b, int64_t& q, int64_t& r) {
  q = a / b;
  r = a % b;
}

}  // namespace

// Using the extended Euclidean algo, we find a and b such that
//     a x + b m = gcd(x, m)
// https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
int64_t ModularInverse(int64_t x, int64_t m) {
  DCHECK_GE(x, 0);
  DCHECK_LT(x, m);

  int64_t r[2] = {m, x};
  int64_t t[2] = {0, 1};
  int64_t q;

  // We only keep the last two terms of the sequences with the "^1" trick:
  //
  // q = r[i-2] / r[i-1]
  // r[i] = r[i-2] % r[i-1]
  // t[i] = t[i-2] - t[i-1] * q
  //
  // We always have:
  // - gcd(r[i], r[i - 1]) = gcd(r[i - 1], r[i - 2])
  // - x * t[i] + m * t[i - 1] = r[i]
  int i = 0;
  for (; r[i ^ 1] != 0; i ^= 1) {
    QuotientAndRemainder(r[i], r[i ^ 1], q, r[i]);
    t[i] -= t[i ^ 1] * q;
  }

  // If the gcd is not one, there is no inverse, we returns 0.
  if (r[i] != 1) return 0;

  // Correct the result so that it is in [0, m). Note that abs(t[i]) is known to
  // be less than or equal to x / 2, and we have thorough unit-tests.
  if (t[i] < 0) t[i] += m;

  return t[i];
}

int64_t PositiveMod(int64_t x, int64_t m) {
  const int64_t r = x % m;
  return r < 0 ? r + m : r;
}

int64_t ProductWithModularInverse(int64_t coeff, int64_t mod, int64_t rhs) {
  DCHECK_NE(coeff, 0);
  DCHECK_NE(mod, 0);

  mod = std::abs(mod);
  if (rhs == 0 || mod == 1) return 0;
  DCHECK_EQ(std::gcd(std::abs(coeff), mod), 1);

  // Make both in [0, mod).
  coeff = PositiveMod(coeff, mod);
  rhs = PositiveMod(rhs, mod);

  // From X * coeff % mod = rhs
  // We deduce that X % mod = rhs * inverse % mod
  const int64_t inverse = ModularInverse(coeff, mod);
  CHECK_NE(inverse, 0);

  // We make the operation in 128 bits to be sure not to have any overflow here.
  const absl::int128 p = absl::int128{inverse} * absl::int128{rhs};
  return static_cast<int64_t>(p % absl::int128{mod});
}

bool SolveDiophantineEquationOfSizeTwo(int64_t& a, int64_t& b, int64_t& cte,
                                       int64_t& x0, int64_t& y0) {
  CHECK_NE(a, 0);
  CHECK_NE(b, 0);
  CHECK_NE(a, std::numeric_limits<int64_t>::min());
  CHECK_NE(b, std::numeric_limits<int64_t>::min());

  const int64_t gcd = std::gcd(std::abs(a), std::abs(b));
  if (cte % gcd != 0) return false;
  a /= gcd;
  b /= gcd;
  cte /= gcd;

  // The simple case where (0, 0) is a solution.
  if (cte == 0) {
    x0 = y0 = 0;
    return true;
  }

  // We solve a * X + b * Y = cte
  // We take a valid x0 in [0, b) by considering the equation mod b.
  x0 = ProductWithModularInverse(a, b, cte);

  // We choose x0 of the same sign as cte.
  if (cte < 0 && x0 != 0) x0 -= std::abs(b);

  // By plugging X = x0 + b * Z
  // We have a * (x0 + b * Z) + b * Y = cte
  // so a * b * Z + b * Y = cte - a * x0;
  // and y0 = (cte - a * x0) / b (with an exact division by construction).
  const absl::int128 t = absl::int128{cte} - absl::int128{a} * absl::int128{x0};
  DCHECK_EQ(t % absl::int128{b}, absl::int128{0});

  // Overflow-wise, there is two cases for cte > 0:
  // - a * x0 <= cte, in this case y0 will not overflow (<= cte).
  // - a * x0 > cte, in this case y0 will be in (-a, 0].
  const absl::int128 r = t / absl::int128{b};
  DCHECK_LE(r, absl::int128{std::numeric_limits<int64_t>::max()});
  DCHECK_GE(r, absl::int128{std::numeric_limits<int64_t>::min()});

  y0 = static_cast<int64_t>(r);
  return true;
}

bool DiophantineEquationOfSizeTwoHasSolutionInDomain(const Domain& x, int64_t a,
                                                     const Domain& y, int64_t b,
                                                     int64_t cte) {
  if (x.IsEmpty() || y.IsEmpty()) return false;
  if (a == 0 && b == 0) {
    return cte == 0;
  }
  if (a == 0) {
    const int64_t div = cte / b;
    if (b * div != cte) return false;
    return y.Contains(div);
  }
  if (b == 0) {
    const int64_t div = cte / a;
    if (a * div != cte) return false;
    return x.Contains(div);
  }

  if (a == b) {
    const int64_t div = cte / a;
    if (a * div != cte) {
      return false;
    }
    return !Domain(x)
                .Negation()
                .AdditionWith(Domain(div))
                .IntersectionWith(y)
                .IsEmpty();
  }

  int64_t x0, y0;
  if (!SolveDiophantineEquationOfSizeTwo(a, b, cte, x0, y0)) {
    return false;
  }
  const Domain z_domain =
      x.AdditionWith(Domain(-x0))
          .InverseMultiplicationBy(b)
          .IntersectionWith(
              y.AdditionWith(Domain(-y0)).InverseMultiplicationBy(-a));

  const Domain z_restricted_d1 =
      x.AdditionWith(Domain(-x0)).InverseMultiplicationBy(b);
  const Domain z_restricted_d2 =
      y.AdditionWith(Domain(-y0)).InverseMultiplicationBy(-a);
  const Domain z_restricted_domain =
      z_restricted_d1.IntersectionWith(z_restricted_d2);
  return !z_restricted_domain.IsEmpty();
}

// TODO(user): Find better implementation? In practice passing via double is
// almost always correct, but the CapProd() might be a bit slow. However this
// is only called when we do propagate something.
int64_t FloorSquareRoot(int64_t a) {
  int64_t result =
      static_cast<int64_t>(std::floor(std::sqrt(static_cast<double>(a))));
  while (CapProd(result, result) > a) --result;
  while (CapProd(result + 1, result + 1) <= a) ++result;
  return result;
}

// TODO(user): Find better implementation?
int64_t CeilSquareRoot(int64_t a) {
  int64_t result =
      static_cast<int64_t>(std::ceil(std::sqrt(static_cast<double>(a))));
  while (CapProd(result, result) < a) ++result;
  while ((result - 1) * (result - 1) >= a) --result;
  return result;
}

int64_t ClosestMultiple(int64_t value, int64_t base) {
  if (value < 0) return -ClosestMultiple(-value, base);
  int64_t result = value / base * base;
  if (value - result > base / 2) result += base;
  return result;
}

bool LinearInequalityCanBeReducedWithClosestMultiple(
    int64_t base, absl::Span<const int64_t> coeffs,
    absl::Span<const int64_t> lbs, absl::Span<const int64_t> ubs, int64_t rhs,
    int64_t* new_rhs) {
  // Precompute some bounds for the equation base * X + error <= rhs.
  int64_t max_activity = 0;
  int64_t max_x = 0;
  int64_t min_error = 0;
  const int num_terms = coeffs.size();
  if (num_terms == 0) return false;
  for (int i = 0; i < num_terms; ++i) {
    const int64_t coeff = coeffs[i];
    CHECK_GT(coeff, 0);
    const int64_t closest = ClosestMultiple(coeff, base);
    max_activity += coeff * ubs[i];
    max_x += closest / base * ubs[i];

    const int64_t error = coeff - closest;
    if (error >= 0) {
      min_error += error * lbs[i];
    } else {
      min_error += error * ubs[i];
    }
  }

  if (max_activity <= rhs) {
    // The constraint is trivially true.
    *new_rhs = max_x;
    return true;
  }

  // This is the max error assuming that activity > rhs.
  int64_t max_error_if_invalid = 0;
  const int64_t slack = max_activity - rhs - 1;
  for (int i = 0; i < num_terms; ++i) {
    const int64_t coeff = coeffs[i];
    const int64_t closest = ClosestMultiple(coeff, base);
    const int64_t error = coeff - closest;
    if (error >= 0) {
      max_error_if_invalid += error * ubs[i];
    } else {
      const int64_t lb = std::max(lbs[i], ubs[i] - slack / coeff);
      max_error_if_invalid += error * lb;
    }
  }

  // We have old solution valid =>
  //     base * X + error <= rhs
  //     base * X <= rhs - error
  //     base * X <= rhs - min_error
  //     X <= new_rhs
  *new_rhs = std::min(max_x, MathUtil::FloorOfRatio(rhs - min_error, base));

  // And we have old solution invalid =>
  //     base * X + error >= rhs + 1
  //     base * X >= rhs + 1 - max_error_if_invalid
  //     X >= infeasibility_bound
  const int64_t infeasibility_bound =
      MathUtil::CeilOfRatio(rhs + 1 - max_error_if_invalid, base);

  // If the two bounds can be separated, we have an equivalence !
  return *new_rhs < infeasibility_bound;
}

int MoveOneUnprocessedLiteralLast(
    const absl::btree_set<LiteralIndex>& processed, int relevant_prefix_size,
    std::vector<Literal>* literals) {
  if (literals->empty()) return -1;
  if (!processed.contains(literals->back().Index())) {
    return std::min<int>(relevant_prefix_size, literals->size());
  }

  // To get O(n log n) size of suffixes, we will first process the last n/2
  // literals, we then move all of them first and process the n/2 literals left.
  // We use the same algorithm recursively. The sum of the suffixes' size S(n)
  // is thus S(n/2) + n + S(n/2). That gives us the correct complexity. The code
  // below simulates one step of this algorithm and is made to be "robust" when
  // from one call to the next, some literals have been removed (but the order
  // of literals is preserved).
  int num_processed = 0;
  int num_not_processed = 0;
  int target_prefix_size = literals->size() - 1;
  for (int i = literals->size() - 1; i >= 0; i--) {
    if (processed.contains((*literals)[i].Index())) {
      ++num_processed;
    } else {
      ++num_not_processed;
      target_prefix_size = i;
    }
    if (num_not_processed >= num_processed) break;
  }
  if (num_not_processed == 0) return -1;
  target_prefix_size = std::min(target_prefix_size, relevant_prefix_size);

  // Once a prefix size has been decided, it is always better to
  // enqueue the literal already processed first.
  std::stable_partition(
      literals->begin() + target_prefix_size, literals->end(),
      [&processed](Literal l) { return processed.contains(l.Index()); });
  return target_prefix_size;
}

int WeightedPick(absl::Span<const double> input, absl::BitGenRef random) {
  DCHECK(!input.empty());

  double total_weight = 0;
  for (const double w : input) {
    total_weight += w;
  }

  const double weight_point = absl::Uniform(random, 0.0f, total_weight);
  double total_weight2 = 0;
  for (int i = 0; i < input.size(); ++i) {
    total_weight2 += input[i];
    if (total_weight2 > weight_point) {
      return i;
    }
  }
  return input.size() - 1;
}

void IncrementalAverage::Reset(double reset_value) {
  num_records_ = 0;
  average_ = reset_value;
}

void IncrementalAverage::AddData(double new_record) {
  num_records_++;
  average_ += (new_record - average_) / num_records_;
}

void ExponentialMovingAverage::AddData(double new_record) {
  num_records_++;
  average_ = (num_records_ == 1)
                 ? new_record
                 : (new_record + decaying_factor_ * (average_ - new_record));
}

void Percentile::AddRecord(double record) {
  records_.push_front(record);
  if (records_.size() > record_limit_) {
    records_.pop_back();
  }
}

double Percentile::GetPercentile(double percent) {
  CHECK_GT(records_.size(), 0);
  CHECK_LE(percent, 100.0);
  CHECK_GE(percent, 0.0);

  const int num_records = records_.size();
  const double percentile_rank =
      static_cast<double>(num_records) * percent / 100.0 - 0.5;
  if (percentile_rank <= 0) {
    return *absl::c_min_element(records_);
  } else if (percentile_rank >= num_records - 1) {
    return *absl::c_max_element(records_);
  }
  std::vector<double> sorted_records;
  sorted_records.assign(records_.begin(), records_.end());
  // Interpolate.
  DCHECK_GE(num_records, 2);
  DCHECK_LT(percentile_rank, num_records - 1);
  const int lower_rank = static_cast<int>(std::floor(percentile_rank));
  DCHECK_LT(lower_rank, num_records - 1);
  auto upper_it = sorted_records.begin() + lower_rank + 1;
  // Ensure that sorted_records[lower_rank + 1] is in the correct place as if
  // records were actually sorted, the next closest is then the largest element
  // to the left of this element.
  absl::c_nth_element(sorted_records, upper_it);
  auto lower_it = std::max_element(sorted_records.begin(), upper_it);
  return *lower_it + (percentile_rank - lower_rank) * (*upper_it - *lower_it);
}

void MaxBoundedSubsetSum::Reset(int64_t bound) {
  DCHECK_GE(bound, 0);
  gcd_ = 0;
  sums_ = {0};
  expanded_sums_.clear();
  current_max_ = 0;
  bound_ = bound;
}

void MaxBoundedSubsetSum::Add(int64_t value) {
  if (value == 0) return;
  if (value > bound_) return;
  gcd_ = std::gcd(gcd_, value);
  AddChoicesInternal({value});
}

void MaxBoundedSubsetSum::AddChoices(absl::Span<const int64_t> choices) {
  if (DEBUG_MODE) {
    for (const int64_t c : choices) {
      DCHECK_GE(c, 0);
    }
  }

  // The max is already reachable or we aborted.
  if (current_max_ == bound_) return;

  // Filter out zero and values greater than bound_.
  filtered_values_.clear();
  for (const int64_t c : choices) {
    if (c == 0 || c > bound_) continue;
    filtered_values_.push_back(c);
    gcd_ = std::gcd(gcd_, c);
  }
  if (filtered_values_.empty()) return;

  // So we can abort early in the AddChoicesInternal() inner loops.
  std::sort(filtered_values_.begin(), filtered_values_.end());
  AddChoicesInternal(filtered_values_);
}

void MaxBoundedSubsetSum::AddMultiples(int64_t coeff, int64_t max_value) {
  DCHECK_GE(coeff, 0);
  DCHECK_GE(max_value, 0);

  if (coeff == 0 || max_value == 0) return;
  if (coeff > bound_) return;
  if (current_max_ == bound_) return;
  gcd_ = std::gcd(gcd_, coeff);

  const int64_t num_values =
      std::min(max_value, MathUtil::FloorOfRatio(bound_, coeff));
  if (num_values > 10) {
    // We only keep GCD in this case.
    sums_.clear();
    expanded_sums_.clear();
    current_max_ = MathUtil::FloorOfRatio(bound_, gcd_) * gcd_;
    return;
  }

  filtered_values_.clear();
  for (int multiple = 1; multiple <= num_values; ++multiple) {
    const int64_t v = multiple * coeff;
    if (v == bound_) {
      current_max_ = bound_;
      return;
    }
    filtered_values_.push_back(v);
  }
  AddChoicesInternal(filtered_values_);
}

void MaxBoundedSubsetSum::AddChoicesInternal(absl::Span<const int64_t> values) {
  // Mode 1: vector of all possible sums (with duplicates).
  if (!sums_.empty() && sums_.size() <= max_complexity_per_add_) {
    const int old_size = sums_.size();
    for (int i = 0; i < old_size; ++i) {
      for (const int64_t value : values) {
        const int64_t s = sums_[i] + value;
        if (s > bound_) break;

        sums_.push_back(s);
        current_max_ = std::max(current_max_, s);
        if (current_max_ == bound_) return;  // Abort
      }
    }
    return;
  }

  // Mode 2: bitset of all possible sums.
  if (bound_ <= max_complexity_per_add_) {
    if (!sums_.empty()) {
      expanded_sums_.assign(bound_ + 1, false);
      for (const int64_t s : sums_) {
        expanded_sums_[s] = true;
      }
      sums_.clear();
    }

    // The reverse order is important to not add the current value twice.
    if (!expanded_sums_.empty()) {
      for (int64_t i = bound_ - 1; i >= 0; --i) {
        if (!expanded_sums_[i]) continue;
        for (const int64_t value : values) {
          if (i + value > bound_) break;

          expanded_sums_[i + value] = true;
          current_max_ = std::max(current_max_, i + value);
          if (current_max_ == bound_) return;  // Abort
        }
      }
      return;
    }
  }

  // Fall back to gcd_.
  DCHECK_NE(gcd_, 0);
  if (gcd_ == 1) {
    current_max_ = bound_;
  } else {
    current_max_ = MathUtil::FloorOfRatio(bound_, gcd_) * gcd_;
  }
}

int64_t MaxBoundedSubsetSum::MaxIfAdded(int64_t candidate) const {
  if (candidate > bound_ || current_max_ == bound_) return current_max_;

  int64_t current_max = current_max_;
  // Mode 1: vector of all possible sums (with duplicates).
  if (!sums_.empty()) {
    for (const int64_t v : sums_) {
      if (v + candidate > bound_) continue;
      if (v + candidate > current_max) {
        current_max = v + candidate;
        if (current_max == bound_) return current_max;
      }
    }
    return current_max;
  }

  // Mode 2: bitset of all possible sums.
  if (!expanded_sums_.empty()) {
    const int64_t min_useful = std::max<int64_t>(0, current_max_ - candidate);
    const int64_t max_useful = bound_ - candidate;
    for (int64_t v = max_useful; v >= min_useful; --v) {
      if (expanded_sums_[v]) return v + candidate;
    }
  }
  return current_max_;
}

BasicKnapsackSolver::Result BasicKnapsackSolver::Solve(
    absl::Span<const Domain> domains, absl::Span<const int64_t> coeffs,
    absl::Span<const int64_t> costs, const Domain& rhs) {
  const int num_vars = domains.size();
  if (num_vars == 0) return {};

  int64_t min_activity = 0;
  int64_t max_domain_size = 0;
  for (int i = 0; i < num_vars; ++i) {
    max_domain_size = std::max(max_domain_size, domains[i].Size());
    if (coeffs[i] > 0) {
      min_activity += coeffs[i] * domains[i].Min();
    } else {
      min_activity += coeffs[i] * domains[i].Max();
    }
  }

  // The complexity of our DP will depends on the number of "activity" values
  // that need to be considered.
  //
  // TODO(user): We can also solve efficiently if max_activity - rhs.Min() is
  // small. Implement.
  const int64_t num_values = rhs.Max() - min_activity + 1;
  if (num_values < 0) {
    // Problem is clearly infeasible, we can report the result right away.
    Result result;
    result.solved = true;
    result.infeasible = true;
    return result;
  }

  // Abort if complexity too large.
  const int64_t max_work_per_variable = std::min(num_values, max_domain_size);
  if (rhs.Max() - min_activity > 1e6) return {};
  if (num_vars * num_values * max_work_per_variable > 1e8) return {};

  // Canonicalize to positive coeffs and non-negative variables.
  domains_.clear();
  coeffs_.clear();
  costs_.clear();
  for (int i = 0; i < num_vars; ++i) {
    if (coeffs[i] > 0) {
      domains_.push_back(domains[i].AdditionWith(Domain(-domains[i].Min())));
      coeffs_.push_back(coeffs[i]);
      costs_.push_back(costs[i]);
    } else {
      domains_.push_back(
          domains[i].Negation().AdditionWith(Domain(domains[i].Max())));
      coeffs_.push_back(-coeffs[i]);
      costs_.push_back(-costs[i]);
    }
  }

  Result result =
      InternalSolve(num_values, rhs.AdditionWith(Domain(-min_activity)));
  if (result.solved && !result.infeasible) {
    // Transform solution back.
    for (int i = 0; i < num_vars; ++i) {
      if (coeffs[i] > 0) {
        result.solution[i] += domains[i].Min();
      } else {
        result.solution[i] = domains[i].Max() - result.solution[i];
      }
    }
  }
  return result;
}

BasicKnapsackSolver::Result BasicKnapsackSolver::InternalSolve(
    int64_t num_values, const Domain& rhs) {
  const int num_vars = domains_.size();

  // The set of DP states that we will fill.
  var_activity_states_.assign(num_vars, std::vector<State>(num_values));

  // Initialize with first variable.
  for (const int64_t v : domains_[0].Values()) {
    const int64_t value = v * coeffs_[0];
    CHECK_GE(value, 0);
    if (value >= num_values) break;
    var_activity_states_[0][value].cost = v * costs_[0];
    var_activity_states_[0][value].value = v;
  }

  // Fill rest of the DP states.
  for (int i = 1; i < num_vars; ++i) {
    const std::vector<State>& prev = var_activity_states_[i - 1];
    std::vector<State>& current = var_activity_states_[i];
    for (int prev_value = 0; prev_value < num_values; ++prev_value) {
      if (prev[prev_value].cost == std::numeric_limits<int64_t>::max()) {
        continue;
      }
      for (const int64_t v : domains_[i].Values()) {
        const int64_t value = prev_value + v * coeffs_[i];
        CHECK_GE(value, 0);
        if (value >= num_values) break;
        const int64_t new_cost = prev[prev_value].cost + v * costs_[i];
        if (new_cost < current[value].cost) {
          current[value].cost = new_cost;
          current[value].value = v;
        }
      }
    }
  }

  Result result;
  result.solved = true;

  int64_t best_cost = std::numeric_limits<int64_t>::max();
  int64_t best_activity;
  for (int v = 0; v < num_values; ++v) {
    // TODO(user): optimize this?
    if (!rhs.Contains(v)) continue;
    if (var_activity_states_.back()[v].cost < best_cost) {
      best_cost = var_activity_states_.back()[v].cost;
      best_activity = v;
    }
  }

  if (best_cost == std::numeric_limits<int64_t>::max()) {
    result.infeasible = true;
    return result;
  }

  // Recover the values.
  result.solution.resize(num_vars);
  int64_t current_activity = best_activity;
  for (int i = num_vars - 1; i >= 0; --i) {
    const int64_t var_value = var_activity_states_[i][current_activity].value;
    result.solution[i] = var_value;
    current_activity -= coeffs_[i] * var_value;
  }

  return result;
}

namespace {

class CliqueDecomposition {
 public:
  CliqueDecomposition(const std::vector<std::vector<int>>& graph,
                      absl::BitGenRef random, std::vector<int>* buffer)
      : graph_(graph), random_(random), buffer_(buffer) {
    const int n = graph.size();
    permutation_.resize(n);
    buffer_->resize(n);

    // TODO(user): Start by sorting by decreasing size of adjacency list?
    for (int i = 0; i < n; ++i) permutation_[i] = i;
  }

  // This works in O(m). All possible decomposition are reachable, depending on
  // the initial permutation.
  //
  // TODO(user): It can be made faster within the same complexity though.
  void DecomposeGreedily() {
    decomposition_.clear();
    candidates_.clear();

    const int n = graph_.size();
    taken_.assign(n, false);
    temp_.assign(n, false);

    int buffer_index = 0;
    for (const int i : permutation_) {
      if (taken_[i]) continue;

      // Save start of amo and output i.
      const int start = buffer_index;
      taken_[i] = true;
      (*buffer_)[buffer_index++] = i;

      candidates_.clear();
      for (const int c : graph_[i]) {
        if (!taken_[c]) candidates_.push_back(c);
      }
      while (!candidates_.empty()) {
        // Extend at most one with lowest possible permutation index.
        int next = candidates_.front();
        for (const int n : candidates_) {
          if (permutation_[n] < permutation_[next]) next = n;
        }

        // Add to current amo.
        taken_[next] = true;
        (*buffer_)[buffer_index++] = next;

        // Filter candidates.
        // TODO(user): With a sorted graph_, same complexity, but likely faster.
        for (const int head : graph_[next]) temp_[head] = true;
        int new_size = 0;
        for (const int c : candidates_) {
          if (taken_[c]) continue;
          if (!temp_[c]) continue;
          candidates_[new_size++] = c;
        }
        candidates_.resize(new_size);
        for (const int head : graph_[next]) temp_[head] = false;
      }

      // Output amo.
      decomposition_.push_back(
          absl::MakeSpan(buffer_->data() + start, buffer_index - start));
    }
    DCHECK_EQ(buffer_index, n);
  }

  // We follow the heuristic in the paper: "Partitioning Networks into Cliques:
  // A Randomized Heuristic Approach", David Chaluppa, 2014.
  //
  // Note that by keeping the local order of each clique, we cannot make the
  // order "worse". And each new call to DecomposeGreedily() can only reduce
  // the number of clique in the cover.
  void ChangeOrder() {
    if (absl::Bernoulli(random_, 0.5)) {
      std::reverse(decomposition_.begin(), decomposition_.end());
    } else {
      std::shuffle(decomposition_.begin(), decomposition_.end(), random_);
    }

    int out_index = 0;
    for (const absl::Span<const int> clique : decomposition_) {
      for (const int i : clique) {
        permutation_[out_index++] = i;
      }
    }
  }

  const std::vector<absl::Span<int>>& decomposition() const {
    return decomposition_;
  }

 private:
  const std::vector<std::vector<int>>& graph_;
  absl::BitGenRef random_;
  std::vector<int>* buffer_;

  std::vector<absl::Span<int>> decomposition_;
  std::vector<int> candidates_;
  std::vector<int> permutation_;
  std::vector<bool> taken_;
  std::vector<bool> temp_;
};

}  // namespace

std::vector<absl::Span<int>> AtMostOneDecomposition(
    const std::vector<std::vector<int>>& graph, absl::BitGenRef random,
    std::vector<int>* buffer) {
  CliqueDecomposition decomposer(graph, random, buffer);
  for (int pass = 0; pass < 10; ++pass) {
    decomposer.DecomposeGreedily();
    if (decomposer.decomposition().size() == 1) break;
    decomposer.ChangeOrder();
  }
  return decomposer.decomposition();
}

absl::Span<const int64_t> SortedSubsetSums::Compute(
    absl::Span<const int64_t> elements, int64_t maximum_sum,
    bool abort_if_maximum_reached) {
  sums_.clear();
  sums_.push_back(0);
  for (const int64_t e : elements) {
    DCHECK_GE(e, 0);
    if (e == 0 || e > maximum_sum) continue;

    // Optimization: If all the sums in [0, maximum_sum] are already reachable
    // we can abort early since no new reachable sum wil be discovered.
    if (sums_.size() == maximum_sum + 1) return sums_;

    // Early abort when asked if we already reached maximum_sum.
    if (abort_if_maximum_reached && sums_.back() == maximum_sum) return sums_;

    // We merge sort sums_ and sums_ + e into new_sums_.
    int i = 0;
    int j = 0;
    new_sums_.clear();
    const int size = sums_.size();
    const int64_t* const data = sums_.data();
    int64_t last_pushed = -1;
    while (i < size) {
      DCHECK_LE(j, i);  // Since we add a positive e.
      const int64_t a = data[i];
      const int64_t b = data[j] + e;
      int64_t to_push;
      if (a <= b) {
        ++i;
        if (a == b) ++j;
        to_push = a;
      } else {
        to_push = b;
        ++j;
      }
      if (to_push == last_pushed) continue;
      if (to_push > maximum_sum) {
        j = size;  // so we don't keep pushing below.
        break;
      }
      last_pushed = to_push;
      new_sums_.push_back(to_push);
    }

    // We are sure of this since we will break only when to_push > maximum_sum
    // and we are guarantee to have pushed all sums below "to_push" before, that
    // includes all the initial sums in sums_.
    DCHECK_EQ(i, size);

    for (; j < size; ++j) {
      const int64_t to_push = data[j] + e;
      if (to_push == last_pushed) continue;
      if (to_push > maximum_sum) break;
      last_pushed = to_push;
      new_sums_.push_back(to_push);
    }
    std::swap(sums_, new_sums_);
  }
  return sums_;
}

double MaxBoundedSubsetSumExact::ComplexityEstimate(int num_elements,
                                                    int64_t bin_size) {
  double estimate = std::numeric_limits<double>::infinity();
  if (num_elements < 100) {
    estimate = 2 * pow(2.0, num_elements / 2);
  }
  return std::min(estimate, static_cast<double>(bin_size) *
                                static_cast<double>(num_elements));
}

int64_t MaxBoundedSubsetSumExact::MaxSubsetSum(
    absl::Span<const int64_t> elements, int64_t bin_size) {
  // Get rid of a few easy to decide corner cases.
  if (elements.empty()) return 0;
  if (elements.size() == 1) {
    if (elements[0] > bin_size) return 0;
    return elements[0];
  }

  // We split the elements in two equal-sized parts.
  const int middle = elements.size() / 2;
  const auto span_a = sums_a_.Compute(elements.subspan(0, middle), bin_size,
                                      /*abort_if_maximum_reached=*/true);
  if (span_a.back() == bin_size) return bin_size;

  const auto span_b = sums_b_.Compute(elements.subspan(middle), bin_size,
                                      /*abort_if_maximum_reached=*/true);
  if (span_b.back() == bin_size) return bin_size;

  // For all possible sum a, we compute the largest sum b that fits.
  // We do that in linear time thanks to the sorted partial sums.
  int64_t result = 0;
  CHECK(!span_a.empty());
  CHECK(!span_b.empty());
  int j = span_b.size() - 1;
  for (int i = 0; i < span_a.size(); ++i) {
    while (j >= 0 && span_a[i] + span_b[j] > bin_size) --j;
    result = std::max(result, span_a[i] + span_b[j]);
    if (result == bin_size) return bin_size;
  }
  return result;
}

std::vector<int> FindMostDiverseSubset(int k, int n,
                                       absl::Span<const int64_t> distances,
                                       std::vector<int64_t>& buffer,
                                       int always_pick_mask) {
  DCHECK_LE(k, n);
  std::vector<int> result;
  result.reserve(k);

  if (k == n) {
    for (int i = 0; i < n; ++i) result.push_back(i);
    return result;
  }

  if (k == n - 1) {
    // We just exclude the one closer to all the other.
    int64_t worse = std::numeric_limits<int64_t>::max();
    int to_exclude = -1;
    for (int i = 0; i < n; ++i) {
      if ((always_pick_mask >> i) & 1) continue;

      int64_t score = 0;
      for (int j = 0; j < n; ++j) {
        if (i != j) score += distances[i * n + j];
      }
      if (score < worse) {
        worse = score;
        to_exclude = i;
      }
    }

    CHECK_NE(to_exclude, -1);
    for (int i = 0; i < n; ++i) {
      if (i != to_exclude) result.push_back(i);
    }
    return result;
  }

  CHECK_LE(n, 25);
  const int limit = 1 << n;
  buffer.assign(limit, 0);
  int best_mask;
  int best_value = -1;
  for (unsigned int mask = 1; mask < limit; ++mask) {
    const int hamming_weight = absl::popcount(mask);

    // TODO(user): Increase mask by more than one ? but counting to 1k is fast
    // anyway.
    if (hamming_weight > k) continue;
    int low_bit = -1;
    int64_t sum = 0;
    for (int i = 0; i < n; ++i) {
      if ((mask >> i) & 1) {
        if (low_bit == -1) {
          low_bit = i;
        } else {
          sum += distances[low_bit * n + i];
        }
      }
    }
    buffer[mask] = buffer[mask ^ (1 << low_bit)] + sum;
    if (hamming_weight == k && buffer[mask] > best_value) {
      if ((mask & always_pick_mask) != always_pick_mask) continue;
      best_value = buffer[mask];
      best_mask = mask;
    }
  }

  for (int i = 0; i < n; ++i) {
    if ((best_mask >> i) & 1) {
      result.push_back(i);
    }
  }
  return result;
}

std::vector<std::pair<int, int>> HeuristicallySplitLongLinear(
    absl::Span<const int64_t> coeffs) {
  std::vector<std::pair<int, int>> result;
  if (coeffs.empty()) return result;

  // Split an interval [0, size) into num_parts mostly equal parts.
  const auto append_splits = [&result](int offset, int size, int num_parts) {
    int previous_start = 0;
    for (int64_t b = 0; b < num_parts; ++b) {
      const int next_start = static_cast<int64_t>(b + 1) * size / num_parts;
      result.push_back({offset + previous_start, next_start - previous_start});
      previous_start = next_start;
    }
  };

  const int num_terms = coeffs.size();
  const int num_buckets = static_cast<int>(std::round(std::sqrt(num_terms)));

  int num_differents = 1;
  for (int i = 1; i < num_terms; ++i) {
    if (coeffs[i - 1] != coeffs[i]) ++num_differents;
  }

  // If we don't have many different coefficients, we always create parts
  // with exactly the same coeffs. We split large part evenly into size /
  // expected_part_size.
  if (num_differents < 20) {
    const int expected_part_size = num_terms / num_buckets;

    for (int i = 0; i < num_terms;) {
      int j = i + 1;
      for (; j < num_terms; ++j) {
        if (coeffs[j] != coeffs[i]) break;
      }

      const int local_size = j - i;
      const int num_local_buckets =
          MathUtil::CeilOfRatio(local_size, expected_part_size);
      append_splits(i, local_size, num_local_buckets);

      i = j;
    }

    return result;
  }

  // Otherwise we just split into num_buckets buckets.
  append_splits(0, num_terms, num_buckets);
  return result;
}

}  // namespace sat
}  // namespace operations_research