knapsack.hpp

#include "noya/knapsack.hpp"

View on GitHub

SourceExpanded

#ifndef NOYA_KNAPSACK_HPP
#define NOYA_KNAPSACK_HPP 1

#include "noya/max_plus_convolution.hpp"

#include <algorithm>
#include <bitset>
#include <limits>
#include <vector>

namespace noya {

/// @brief 0/1 knapsack. items are (weight, value) pairs. O(NlogN + L^2).
template <class P>
std::vector<int64_t> knapsack(int L, const std::vector<P> &items) {
  std::vector<std::vector<int>> buckets(L + 1);
  for (auto &[w, v] : items) {
    assert(w >= 0);
    assert(v >= 0);
    if (w <= L) {
      buckets[w].push_back(v);
    }
  }

  std::vector<int64_t> dp(L + 1, std::numeric_limits<int64_t>::lowest());
  dp[0] = 0;

  for (int w = 1; w <= L; w++) {
    auto &bucket = buckets[w];
    if (bucket.empty()) {
      continue;
    }
    std::sort(bucket.begin(), bucket.end(), std::greater<>());
    const int m = std::min(int(bucket.size()), L / w);
    std::vector<int64_t> sum(m + 1);
    for (int i = 0; i < m; i++) {
      sum[i + 1] = sum[i] + bucket[i];
    }
    // remainder class enumeration
    for (int r = 0; r < w; r++) {
      const int n = int(L - r) / w + 1;
      std::vector<int64_t> v(n);
      for (int i = 0; i < n; i++) {
        v[i] = dp[i * w + r];
      }
      const std::vector<int64_t> &res = concave_maxplus_convolution(v, sum);
      for (int i = 0; i < n; i++) {
        dp[i * w + r] = res[i];
      }
    }
  }
  return dp;
}

// https://qoj.ac/problem/7403
// O(n max{A})
template <class T> T max_subsetsum_leq(const T &C, const std::vector<int> &A) {
  int N = int(A.size());
  int p = -1;
  T cur = 0;
  for (int i = 0; i < N; i++) {
    if (cur + A[i] > C) {
      p = i;
      break;
    } else {
      cur += A[i];
    }
  }
  if (p == -1) {
    return cur;
  }
  const int shift = *std::max_element(A.begin(), A.end());
  std::vector<int> dp0(2 * shift + 1);
  std::vector<int> dp1(2 * shift + 1);

  for (int i = 0; i <= shift; i++)
    dp0[i] = -1;

  dp0[cur - C + shift] = p;

  for (int i = p; i < N; i++) {
    dp1 = dp0;
    for (int j = 0; j <= shift; j++)
      dp1[j + A[i]] = std::max(dp1[j + A[i]], dp0[j]);

    for (int j = shift + A[i]; j >= shift + 1; j--)
      for (int k = dp1[j] - 1; k >= dp0[j]; k--)
        dp1[j - A[k]] = std::max(dp1[j - A[k]], k);
    dp0 = dp1;
  }

  for (int i = shift; i >= 0; i--)
    if (dp0[i] >= 0)
      return i - shift + C;
  return 0;
}

template <int N> std::bitset<N> bool_knapsack(const std::vector<int> &items) {
  std::vector<int> freq(N);
  for (auto &item : items) {
    if (item < N) {
      freq[item] += 1;
    }
  }
  for (int i = 1; i < N; i++) {
    if (freq[i] >= 3) {
      int a = (freq[i] - 1) / 2;
      freq[i * 2] += a;
      freq[i] -= a * 2;
    }
  }
  std::bitset<N> res;
  res[0] = 1;
  for (int i = 1; i < N; i++) {
    for (int j = 0; j < freq[i]; j++) {
      res |= res << i;
    }
  }
  return res;
}

} // namespace noya

#endif // NOYA_KNAPSACK_HPP

#include <algorithm>
#include <bitset>
#include <functional>
#include <limits>
#include <numeric>
#include <vector>

// https://noshi91.github.io/Library/algorithm/smawk.cpp.html

namespace noya {

/// @brief SMAWK algorithm: compute row minima of a totally monotone matrix.
/// @return Vector where ans[i] is the column index of the minimum in row i.
template <class Select>
std::vector<int> smawk(const int row_size, const int col_size,
                       const Select &select) {
  const std::function<std::vector<int>(const std::vector<int> &,
                                       const std::vector<int> &)>
      solve = [&](const std::vector<int> &row,
                  const std::vector<int> &col) -> std::vector<int> {
    const int n = int(row.size());
    if (n == 0)
      return {};
    std::vector<int> c2;
    for (const int i : col) {
      while (!c2.empty() && select(row[c2.size() - 1], c2.back(), i))
        c2.pop_back();
      if (c2.size() < n)
        c2.push_back(i);
    }
    std::vector<int> r2;
    for (int i = 1; i < n; i += 2)
      r2.push_back(row[i]);
    const std::vector<int> a2 = solve(r2, c2);
    std::vector<int> ans(n);
    for (int i = 0; i != a2.size(); i += 1)
      ans[i * 2 + 1] = a2[i];
    int j = 0;
    for (int i = 0; i < n; i += 2) {
      ans[i] = c2[j];
      const int end = i + 1 == n ? c2.back() : ans[i + 1];
      while (c2[j] != end) {
        j += 1;
        if (select(row[i], ans[i], c2[j]))
          ans[i] = c2[j];
      }
    }
    return ans;
  };
  std::vector<int> row(row_size);
  std::iota(row.begin(), row.end(), 0);
  std::vector<int> col(col_size);
  std::iota(col.begin(), col.end(), 0);
  return solve(row, col);
}

} // namespace noya

namespace noya {

/// @brief Max-plus convolution of two concave sequences.
template <class T>
std::vector<T> two_concave_maxplus_convolution(const std::vector<T> &a,
                                               const std::vector<T> &b) {
  if (a.empty())
    return b;
  if (b.empty())
    return a;
  const int n = int(a.size());
  const int m = int(b.size());
  int p = 0, q = 0;
  std::vector<T> c(n + m - 1);
  c[0] = a[0] + b[0];
  for (int i = 1; i < n + m - 1; i++) {
    if (p + 1 == n) {
      q++;
    } else if (q + 1 == m) {
      p++;
    } else {
      if (a[p + 1] - a[p] > b[q + 1] - b[q]) {
        p++;
      } else {
        q++;
      }
    }
    c[i] = a[p] + b[q];
  }
  return c;
}

/// @brief Min-plus convolution of two concave sequences.
template <class T>
std::vector<T> two_concave_minplus_convolution(const std::vector<T> &a,
                                               const std::vector<T> &b) {
  const int n = int(a.size());
  const int m = int(b.size());
  std::vector<T> negative_a(n);
  std::vector<T> negative_b(m);
  for (int i = 0; i < n; i++) {
    negative_a[i] = -a[i];
  }
  for (int i = 0; i < m; i++) {
    negative_b[i] = -b[i];
  }
  auto negative_c = two_concave_maxplus_convolution(negative_a, negative_b);
  std::vector<T> c(n + m - 1);
  for (int i = 0; i < n + m - 1; i++)
    c[i] = -negative_c[i];
  return c;
}

/// @brief Max-plus convolution where b is concave.
template <class T>
std::vector<T> concave_maxplus_convolution(const std::vector<T> &a,
                                           const std::vector<T> &b) {
  if (a.empty())
    return b;
  if (b.empty())
    return a;
  const int n = int(a.size());
  const int m = int(b.size());
  const auto get = [&](const int &i, const int &j) -> T {
    return a[j] + b[i - j];
  };
  const auto select = [&](const int &i, const int &j, const int &k) -> bool {
    if (i < k)
      return false;
    if (i - j >= m)
      return true;
    return get(i, j) <= get(i, k);
  };
  const auto amax = smawk(n + m - 1, n, select);
  std::vector<T> c(n + m - 1);
  for (int i = 0; i < n + m - 1; i++)
    c[i] = get(i, amax[i]);
  return c;
}

/// @brief Min-plus convolution where b is concave.
template <class T>
std::vector<T> concave_minplus_convolution(const std::vector<T> &a,
                                           const std::vector<T> &b) {
  const int n = int(a.size());
  const int m = int(b.size());
  std::vector<T> negative_a(n);
  std::vector<T> negative_b(m);
  for (int i = 0; i < n; i++) {
    negative_a[i] = -a[i];
  }
  for (int i = 0; i < m; i++) {
    negative_b[i] = -b[i];
  }
  auto negative_c = concave_maxplus_convolution(negative_a, negative_b);
  std::vector<T> c(n + m - 1);
  for (int i = 0; i < n + m - 1; i++)
    c[i] = -negative_c[i];
  return c;
}

} // namespace noya

namespace noya {

/// @brief 0/1 knapsack. items are (weight, value) pairs. O(NlogN + L^2).
template <class P>
std::vector<int64_t> knapsack(int L, const std::vector<P> &items) {
  std::vector<std::vector<int>> buckets(L + 1);
  for (auto &[w, v] : items) {
    assert(w >= 0);
    assert(v >= 0);
    if (w <= L) {
      buckets[w].push_back(v);
    }
  }

  std::vector<int64_t> dp(L + 1, std::numeric_limits<int64_t>::lowest());
  dp[0] = 0;

  for (int w = 1; w <= L; w++) {
    auto &bucket = buckets[w];
    if (bucket.empty()) {
      continue;
    }
    std::sort(bucket.begin(), bucket.end(), std::greater<>());
    const int m = std::min(int(bucket.size()), L / w);
    std::vector<int64_t> sum(m + 1);
    for (int i = 0; i < m; i++) {
      sum[i + 1] = sum[i] + bucket[i];
    }
    // remainder class enumeration
    for (int r = 0; r < w; r++) {
      const int n = int(L - r) / w + 1;
      std::vector<int64_t> v(n);
      for (int i = 0; i < n; i++) {
        v[i] = dp[i * w + r];
      }
      const std::vector<int64_t> &res = concave_maxplus_convolution(v, sum);
      for (int i = 0; i < n; i++) {
        dp[i * w + r] = res[i];
      }
    }
  }
  return dp;
}

// https://qoj.ac/problem/7403
// O(n max{A})
template <class T> T max_subsetsum_leq(const T &C, const std::vector<int> &A) {
  int N = int(A.size());
  int p = -1;
  T cur = 0;
  for (int i = 0; i < N; i++) {
    if (cur + A[i] > C) {
      p = i;
      break;
    } else {
      cur += A[i];
    }
  }
  if (p == -1) {
    return cur;
  }
  const int shift = *std::max_element(A.begin(), A.end());
  std::vector<int> dp0(2 * shift + 1);
  std::vector<int> dp1(2 * shift + 1);

  for (int i = 0; i <= shift; i++)
    dp0[i] = -1;

  dp0[cur - C + shift] = p;

  for (int i = p; i < N; i++) {
    dp1 = dp0;
    for (int j = 0; j <= shift; j++)
      dp1[j + A[i]] = std::max(dp1[j + A[i]], dp0[j]);

    for (int j = shift + A[i]; j >= shift + 1; j--)
      for (int k = dp1[j] - 1; k >= dp0[j]; k--)
        dp1[j - A[k]] = std::max(dp1[j - A[k]], k);
    dp0 = dp1;
  }

  for (int i = shift; i >= 0; i--)
    if (dp0[i] >= 0)
      return i - shift + C;
  return 0;
}

template <int N> std::bitset<N> bool_knapsack(const std::vector<int> &items) {
  std::vector<int> freq(N);
  for (auto &item : items) {
    if (item < N) {
      freq[item] += 1;
    }
  }
  for (int i = 1; i < N; i++) {
    if (freq[i] >= 3) {
      int a = (freq[i] - 1) / 2;
      freq[i * 2] += a;
      freq[i] -= a * 2;
    }
  }
  std::bitset<N> res;
  res[0] = 1;
  for (int i = 1; i < N; i++) {
    for (int j = 0; j < freq[i]; j++) {
      res |= res << i;
    }
  }
  return res;
}

} // namespace noya