knapsack.hpp¶
计算 0/1 背包在每个容量下的最优价值，并对同重量物品用凹 max-plus 卷积加速。
Complexity: Time: O(N log N + L^2), where L is the weight limit. Space: O(L).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(N log N + L^2), where L is the weight limit.
/// Space: O(L).
#include "noya/max_plus_convolution.hpp"
#include <algorithm>
#include <bitset>
#include <limits>
#include <vector>
namespace noya {
/// @brief Return best 0/1-knapsack values for every capacity through L.
/// Items are grouped by equal weight and their values are sorted. Prefix sums
/// inside a group form a concave sequence, so each residue class of capacities
/// is updated by one concave max-plus convolution instead of per-item DP.
template <class P>
std::vector<int64_t> knapsack(int L, const std::vector<P> &xs) {
std::vector<std::vector<int>> bkt(L + 1);
for (auto &[w, v] : xs) {
assert(w >= 0);
assert(v >= 0);
if (w <= L) {
bkt[w].push_back(v);
}
}
std::vector<int64_t> dp(L + 1, std::numeric_limits<int64_t>::lowest());
dp[0] = 0;
for (int val : bkt[0]) {
dp[0] += val;
}
for (int w = 1; w <= L; w++) {
auto &bk1 = bkt[w];
if (bk1.empty()) {
continue;
}
std::sort(bk1.begin(), bk1.end(), std::greater<>());
const int m = std::min(int(bk1.size()), L / w);
std::vector<int64_t> sum(m + 1);
for (int i = 0; i < m; i++) {
sum[i + 1] = sum[i] + bk1[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;
}
/// @brief Largest attainable subset sum not exceeding C in O(n max(A)).
/// Start from the first prefix that exceeds C and keep only deviations within
/// max(A) of C. The DP exchanges later chosen items with earlier unchosen ones;
/// any optimal repair crosses this narrow boundary, so larger deficits need
/// not be represented.
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 sh = *std::max_element(A.begin(), A.end());
std::vector<int> dp0(2 * sh + 1);
std::vector<int> dp1(2 * sh + 1);
for (int i = 0; i <= sh; i++)
dp0[i] = -1;
dp0[cur - C + sh] = p;
for (int i = p; i < N; i++) {
dp1 = dp0;
for (int j = 0; j <= sh; j++)
dp1[j + A[i]] = std::max(dp1[j + A[i]], dp0[j]);
for (int j = sh + A[i]; j >= sh + 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 = sh; i >= 0; i--)
if (dp0[i] >= 0)
return i - sh + C;
return 0;
}
template <int N> std::bitset<N> bool_knapsack(const std::vector<int> &xs) {
std::vector<int> cnt(N);
for (auto &x : xs) {
if (x < N) {
cnt[x] += 1;
}
}
for (int i = 1; i < N; i++) {
if (cnt[i] >= 3) {
int a = (cnt[i] - 1) / 2;
cnt[i * 2] += a;
cnt[i] -= a * 2;
}
}
std::bitset<N> res;
res[0] = 1;
for (int i = 1; i < N; i++) {
for (int j = 0; j < cnt[i]; j++) {
res |= res << i;
}
}
return res;
}
} // namespace noya
#ifndef NOYA_KNAPSACK_HPP
#define NOYA_KNAPSACK_HPP 1
/// @complexity Time: O(N log N + L^2), where L is the weight limit.
/// Space: O(L).
#include "noya/max_plus_convolution.hpp"
#include <algorithm>
#include <bitset>
#include <limits>
#include <vector>
namespace noya {
/// @brief Return best 0/1-knapsack values for every capacity through L.
/// Items are grouped by equal weight and their values are sorted. Prefix sums
/// inside a group form a concave sequence, so each residue class of capacities
/// is updated by one concave max-plus convolution instead of per-item DP.
template <class P>
std::vector<int64_t> knapsack(int L, const std::vector<P> &xs) {
std::vector<std::vector<int>> bkt(L + 1);
for (auto &[w, v] : xs) {
assert(w >= 0);
assert(v >= 0);
if (w <= L) {
bkt[w].push_back(v);
}
}
std::vector<int64_t> dp(L + 1, std::numeric_limits<int64_t>::lowest());
dp[0] = 0;
for (int val : bkt[0]) {
dp[0] += val;
}
for (int w = 1; w <= L; w++) {
auto &bk1 = bkt[w];
if (bk1.empty()) {
continue;
}
std::sort(bk1.begin(), bk1.end(), std::greater<>());
const int m = std::min(int(bk1.size()), L / w);
std::vector<int64_t> sum(m + 1);
for (int i = 0; i < m; i++) {
sum[i + 1] = sum[i] + bk1[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;
}
/// @brief Largest attainable subset sum not exceeding C in O(n max(A)).
/// Start from the first prefix that exceeds C and keep only deviations within
/// max(A) of C. The DP exchanges later chosen items with earlier unchosen ones;
/// any optimal repair crosses this narrow boundary, so larger deficits need
/// not be represented.
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 sh = *std::max_element(A.begin(), A.end());
std::vector<int> dp0(2 * sh + 1);
std::vector<int> dp1(2 * sh + 1);
for (int i = 0; i <= sh; i++)
dp0[i] = -1;
dp0[cur - C + sh] = p;
for (int i = p; i < N; i++) {
dp1 = dp0;
for (int j = 0; j <= sh; j++)
dp1[j + A[i]] = std::max(dp1[j + A[i]], dp0[j]);
for (int j = sh + A[i]; j >= sh + 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 = sh; i >= 0; i--)
if (dp0[i] >= 0)
return i - sh + C;
return 0;
}
template <int N> std::bitset<N> bool_knapsack(const std::vector<int> &xs) {
std::vector<int> cnt(N);
for (auto &x : xs) {
if (x < N) {
cnt[x] += 1;
}
}
for (int i = 1; i < N; i++) {
if (cnt[i] >= 3) {
int a = (cnt[i] - 1) / 2;
cnt[i * 2] += a;
cnt[i] -= a * 2;
}
}
std::bitset<N> res;
res[0] = 1;
for (int i = 1; i < N; i++) {
for (int j = 0; j < cnt[i]; j++) {
res |= res << i;
}
}
return res;
}
} // namespace noya
#endif // NOYA_KNAPSACK_HPP
#include <algorithm>
#include <bitset>
#include <cassert>
#include <cstdint>
#include <functional>
#include <limits>
#include <numeric>
#include <vector>
/// @complexity Time: O(N log N + L^2), where L is the weight limit.
/// Space: O(L).
/// @complexity Time: O(n + m) for the supported convex/concave cases.
/// Space: O(n + m) output and monotone-optimum workspace.
/// @complexity Time: O(rows + columns) matrix probes.
/// Space: O(rows + columns).
namespace noya {
/// @brief SMAWK algorithm: compute row minima of a totally monotone matrix.
/// A stack reduction leaves at most one candidate column per row, recursion
/// solves the odd rows, and monotone argmins bound the scan that interpolates
/// each even row. Every row and column is discarded or scanned only O(1) times.
/// @return Vector where ans[i] is the column index of the minimum in row i.
template <class Select>
std::vector<int> smawk(const int nr, const int nc, const Select &sel) {
const std::function<std::vector<int>(const std::vector<int> &,
const std::vector<int> &)>
sol = [&](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() && sel(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 = sol(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 (sel(row[i], ans[i], c2[j]))
ans[i] = c2[j];
}
}
return ans;
};
std::vector<int> row(nr);
std::iota(row.begin(), row.end(), 0);
std::vector<int> col(nc);
std::iota(col.begin(), col.end(), 0);
return sol(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 convex sequences.
/// The first differences of a convex sequence are nondecreasing. Merging the
/// two difference sequences therefore describes, in order, every step of the
/// lower boundary of their min-plus convolution.
template <class T>
std::vector<T> convex_convex_minplus_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());
std::vector<T> na(n);
std::vector<T> nb(m);
for (int i = 0; i < n; i++) {
na[i] = -a[i];
}
for (int i = 0; i < m; i++) {
nb[i] = -b[i];
}
auto nc = two_concave_maxplus_convolution(na, nb);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = -nc[i];
return c;
}
/// @brief Backward-compatible name for convex-convex min-plus convolution.
template <class T>
std::vector<T> two_concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return convex_convex_minplus_convolution(a, b);
}
/// @brief Max-plus convolution where b is concave.
/// Concavity makes the implicit matrix a[i] + b[row-i] totally monotone, so
/// SMAWK finds all row maxima after only O(n + m) value comparisons.
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 sel = [&](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 amx = smawk(n + m - 1, n, sel);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = get(i, amx[i]);
return c;
}
/// @brief Min-plus convolution of an arbitrary sequence and a convex sequence.
/// Convexity makes the implicit matrix a[i] + b[row-i] totally monotone.
/// Negating both inputs turns row minima into row maxima, which SMAWK finds in
/// linear time without materializing the matrix.
template <class T>
std::vector<T> arbitrary_convex_minplus_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());
std::vector<T> na(n);
std::vector<T> nb(m);
for (int i = 0; i < n; i++) {
na[i] = -a[i];
}
for (int i = 0; i < m; i++) {
nb[i] = -b[i];
}
auto nc = concave_maxplus_convolution(na, nb);
std::vector<T> c(n + m - 1);
for (int i = 0; i < n + m - 1; i++)
c[i] = -nc[i];
return c;
}
/// @brief Min-plus convolution of a convex sequence and an arbitrary sequence.
template <class T>
std::vector<T> convex_arbitrary_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return arbitrary_convex_minplus_convolution(b, a);
}
namespace max_plus_convolution_internal {
template <class Value>
std::vector<int> monotone_row_minima(int rs, int cs, Value val) {
std::vector<int> res(rs);
int str = 1;
while (str < rs) {
str <<= 1;
}
for (; str > 0; str >>= 1) {
for (int row = str - 1; row < rs; row += 2 * str) {
int lhs = row >= str ? res[row - str] : 0;
int las = row + str < rs ? res[row + str] : cs - 1;
res[row] = lhs;
for (int col = lhs + 1; col <= las; col++) {
if (val(row, col) < val(row, res[row])) {
res[row] = col;
}
}
}
}
return res;
}
template <class T>
std::vector<T>
arbitrary_concave_minplus_convolution(const std::vector<T> &arb,
const std::vector<T> &con) {
if (arb.empty()) {
return con;
}
if (con.empty()) {
return arb;
}
for (int i = 0; i + 2 < int(con.size()); i++) {
assert(con[i + 1] - con[i] >= con[i + 2] - con[i + 1]);
}
const int wid = int(arb.size());
const int sz = int(con.size());
const int hei = wid + sz - 1;
std::vector<int> lc(hei), mc(hei, wid - 1);
for (int row = sz; row < hei; row++) {
lc[row] = row - sz + 1;
}
for (int row = 0; row <= hei - sz; row++) {
mc[row] = row;
}
std::vector<int> lr(wid), mr(wid);
for (int col = 0; col < wid; col++) {
lr[col] = col;
mr[col] = sz - 1 + col;
}
std::vector<T> res(hei, std::numeric_limits<T>::max());
std::function<void(int, int, int, int)> sol = [&](int rl, int lro, int cl,
int cr) {
if (mc[rl] >= cr && cl >= lc[lro]) {
auto val = [&](int ri, int rc) {
int col = cr - rc;
int row = rl + ri;
return arb[col] + con[row - col];
};
auto mnm = monotone_row_minima(lro - rl + 1, cr - cl + 1, val);
for (int row = rl; row <= lro; row++) {
res[row] = std::min(res[row], val(row - rl, mnm[row - rl]));
}
return;
}
if (std::int64_t(lro - rl) * (cr - cl) < 1024) {
for (int row = rl; row <= lro; row++) {
int fro = std::max(lc[row], cl);
int to = std::min(mc[row], cr);
for (int col = fro; col <= to; col++) {
res[row] = std::min(res[row], arb[col] + con[row - col]);
}
}
return;
}
if (lro - rl > cr - cl) {
int mid = (rl + lro) / 2;
int nl = std::min(mc[mid], cr);
if (cl <= nl) {
sol(rl, mid, cl, nl);
}
int fst = std::max(lc[mid], cl);
if (fst <= cr) {
sol(mid + 1, lro, fst, cr);
}
} else {
int mid = (cl + cr) / 2;
int nl = std::min(mr[mid], lro);
if (rl <= nl) {
sol(rl, nl, cl, mid);
}
int fst = std::max(lr[mid], rl);
if (fst <= lro) {
sol(fst, lro, mid + 1, cr);
}
}
};
sol(0, hei - 1, 0, wid - 1);
return res;
}
} // namespace max_plus_convolution_internal
/// @brief Min-plus convolution of a concave sequence and an arbitrary
/// sequence. Valid pairs form a diagonal staircase rather than one rectangular
/// Monge matrix. Recursively splitting that staircase produces fully valid
/// rectangles; after reversing their columns, concavity makes row minima
/// monotone and they are found together. Small boundary rectangles are scanned
/// directly.
template <class T>
std::vector<T> concave_arbitrary_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return max_plus_convolution_internal::arbitrary_concave_minplus_convolution(
b, a);
}
/// @brief Backward-compatible name for arbitrary-convex min-plus convolution.
template <class T>
std::vector<T> concave_minplus_convolution(const std::vector<T> &a,
const std::vector<T> &b) {
return arbitrary_convex_minplus_convolution(a, b);
}
} // namespace noya
namespace noya {
/// @brief Return best 0/1-knapsack values for every capacity through L.
/// Items are grouped by equal weight and their values are sorted. Prefix sums
/// inside a group form a concave sequence, so each residue class of capacities
/// is updated by one concave max-plus convolution instead of per-item DP.
template <class P>
std::vector<int64_t> knapsack(int L, const std::vector<P> &xs) {
std::vector<std::vector<int>> bkt(L + 1);
for (auto &[w, v] : xs) {
assert(w >= 0);
assert(v >= 0);
if (w <= L) {
bkt[w].push_back(v);
}
}
std::vector<int64_t> dp(L + 1, std::numeric_limits<int64_t>::lowest());
dp[0] = 0;
for (int val : bkt[0]) {
dp[0] += val;
}
for (int w = 1; w <= L; w++) {
auto &bk1 = bkt[w];
if (bk1.empty()) {
continue;
}
std::sort(bk1.begin(), bk1.end(), std::greater<>());
const int m = std::min(int(bk1.size()), L / w);
std::vector<int64_t> sum(m + 1);
for (int i = 0; i < m; i++) {
sum[i + 1] = sum[i] + bk1[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;
}
/// @brief Largest attainable subset sum not exceeding C in O(n max(A)).
/// Start from the first prefix that exceeds C and keep only deviations within
/// max(A) of C. The DP exchanges later chosen items with earlier unchosen ones;
/// any optimal repair crosses this narrow boundary, so larger deficits need
/// not be represented.
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 sh = *std::max_element(A.begin(), A.end());
std::vector<int> dp0(2 * sh + 1);
std::vector<int> dp1(2 * sh + 1);
for (int i = 0; i <= sh; i++)
dp0[i] = -1;
dp0[cur - C + sh] = p;
for (int i = p; i < N; i++) {
dp1 = dp0;
for (int j = 0; j <= sh; j++)
dp1[j + A[i]] = std::max(dp1[j + A[i]], dp0[j]);
for (int j = sh + A[i]; j >= sh + 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 = sh; i >= 0; i--)
if (dp0[i] >= 0)
return i - sh + C;
return 0;
}
template <int N> std::bitset<N> bool_knapsack(const std::vector<int> &xs) {
std::vector<int> cnt(N);
for (auto &x : xs) {
if (x < N) {
cnt[x] += 1;
}
}
for (int i = 1; i < N; i++) {
if (cnt[i] >= 3) {
int a = (cnt[i] - 1) / 2;
cnt[i * 2] += a;
cnt[i] -= a * 2;
}
}
std::bitset<N> res;
res[0] = 1;
for (int i = 1; i < N; i++) {
for (int j = 0; j < cnt[i]; j++) {
res |= res << i;
}
}
return res;
}
} // namespace noya