Dynamic Programming II: Bitmasks, Intervals, Trees and Digits
Four state shapes that cover most hard DP problems (a set as a bitmask, an interval, a subtree and a digit position), each with code, an operation count and a stress test.
When the state has to remember more than a number
In Dynamic Programming I the state was almost always one or two indices: "the best answer for the first i items", "the best answer for prefixes of length i and j". That works when the past can be summarized by how far you got. Many problems refuse to be summarized that way.
A delivery drone has to visit 15 drop points and come back. After visiting some of them, what matters for the rest of the trip is which points are done and where the drone is now, not the order in which it got there. A splicing machine joins neighboring recordings on a tape; the cheapest way to finish a stretch of tape depends on where the stretch starts and ends. A network of pipes forms a tree; what you can do at a junction depends on what was decided in the subtree below it. A question about every integer up to 10^18 cannot loop over them; it has to work one digit at a time.
Those are the four shapes in this lesson: subsets, intervals, subtrees and digit positions. The method from lesson 6 does not change (state, transition, order, base cases, answer); what changes is what a state looks like, and with it the table size and the loop order. Every solution below is checked against a brute force with the seeded Lcg generator and the stress-test habit from Fast I/O, Contest Templates and Stress Testing. The next lesson, Binary Search on the Answer, Two Pointers and Sweeping, turns to techniques that replace a DP when the answer itself can be searched.
A set of up to 20 items fits in one integer
A subset of {0, 1, ..., n-1} can be stored as an integer whose bit i is 1 exactly when item i is in the set. With n up to about 20, there are 2^n subsets (about a million), and every subset is a plain array index. The operations you need are one or two machine instructions each:
- Is item
iinmask?mask >> i & 1 - Add item
i:mask | 1 << i - Remove item
i:mask & ~(1 << i) - How many items?
std::popcount(mask), which accepts only unsigned types - The full set:
(1 << n) - 1
The drone problem in this form, which is the travelling salesman problem on a small graph: stops are numbered 0 to n-1, stop 0 is the depot, t[a][b] is the flying time from a to b (not necessarily equal to t[b][a], because of wind), and the drone must visit every stop exactly once and return to the depot. Constraints: 2 <= n <= 16, times up to 10^6.
Trying every order costs (n-1)! routes, which is about 1.3*10^12 at n = 16. The observation that saves us: two partial routes that have visited the same set of stops and stand at the same stop have the same future. Only the cheaper one can be part of an optimal tour. So the state is
best[mask][v] = the least time to leave the depot, visit exactly the stops in mask, and stand at v.
By hand, with four stops: the base case is best[{0}][0] = 0. One step later, best[{0,2}][2] = t[0][2] and best[{0,1}][1] = t[0][1]. The state best[{0,1,2}][1] has only one predecessor, best[{0,2}][2], because standing at 1 after visiting {0,1,2} means 1 came last. But best[{0,1,2,3}][3] can be reached from best[{0,1,2}][1] and from best[{0,1,2}][2], and each of those already summarizes every order that led to it. We keep the minimum and forget the rest. The transition is "extend by one unvisited stop w": best[mask | 1 << w][w] = min(..., best[mask][v] + t[v][w]). Every transition goes from a mask to a strictly larger integer, so looping over masks in increasing order computes every state before it is used. The answer is min over v of best[full][v] + t[v][0].
This dynamic program is usually credited to Bellman and, independently, Held and Karp, all in 1962. It runs in O(n^2 2^n) time and needs O(n 2^n) memory.
The drone route in code
The program follows the course's convention: solve reads from any std::istream and writes to any std::ostream, and main feeds it an embedded sample. On a judge you call solve(std::cin, std::cout) and nothing else changes. Next to the cost it stores, for every state, the stop we came from, so the route can be walked backwards at the end. The input is n followed by the n by n time matrix:
5
0 4 9 7 3
6 0 2 8 5
9 3 0 4 6
5 7 2 0 8
2 6 7 3 0#include <algorithm>
#include <iostream>
#include <limits>
#include <sstream>
#include <vector>
using ll = long long;
constexpr ll kInf = std::numeric_limits<ll>::max() / 4;
void solve(std::istream& in, std::ostream& out) {
int n;
in >> n;
std::vector<std::vector<ll>> t(n, std::vector<ll>(n));
for (auto& row : t)
for (ll& x : row) in >> x;
const int full = (1 << n) - 1;
// best[mask][v]: least time to leave 0, visit exactly the stops in mask, stand at v
std::vector<std::vector<ll>> best(1 << n, std::vector<ll>(n, kInf));
std::vector<std::vector<int>> from(1 << n, std::vector<int>(n, -1));
best[1][0] = 0;
for (int mask = 1; mask <= full; ++mask) {
if (!(mask & 1)) continue; // every route starts at the depot
for (int v = 0; v < n; ++v) {
if (best[mask][v] == kInf) continue;
for (int w = 0; w < n; ++w) {
if (mask & (1 << w)) continue; // already visited
int next = mask | (1 << w);
ll cand = best[mask][v] + t[v][w];
if (cand < best[next][w]) {
best[next][w] = cand;
from[next][w] = v;
}
}
}
}
ll answer = kInf;
int last = -1;
for (int v = 1; v < n; ++v)
if (best[full][v] + t[v][0] < answer) answer = best[full][v] + t[v][0], last = v;
std::vector<int> route{0};
for (int mask = full, v = last; v != 0;) { // walk the choices backwards
route.push_back(v);
int prev = from[mask][v];
mask ^= 1 << v;
v = prev;
}
route.push_back(0);
std::ranges::reverse(route);
out << answer << '\n';
for (int v : route) out << v << ' ';
out << '\n';
}
int main() {
std::istringstream in(R"(5
0 4 9 7 3
6 0 2 8 5
9 3 0 4 6
5 7 2 0 8
2 6 7 3 0)");
solve(in, std::cout); // on a judge: solve(std::cin, std::cout)
}Press Run to compile and run this example.
17
0 4 3 2 1 0Check it by hand: t[0][4] + t[4][3] + t[3][2] + t[2][1] + t[1][0] = 3 + 3 + 2 + 3 + 6 = 17. The constraint n >= 2 matters: with a single stop there is no v >= 1, last stays -1 and the reconstruction would index out of bounds. Handle tiny cases explicitly whenever the constraints allow them.
Stress test, then count the work
One sample proves little. The next program runs the DP against a brute force that tries every order with std::ranges::next_permutation, on 3,000 random matrices with n from 2 to 7. A third of the cases use weights 0 and 1 only, so there are many ties, which is where "keep the minimum" logic tends to break. Then it counts relaxations (the innermost loop body) for growing n and prints them next to n^2 2^n / 8 and the number of orders the brute force would try.
#include <algorithm>
#include <cstdint>
#include <limits>
#include <numeric>
#include <print>
#include <vector>
using ll = long long;
constexpr ll kInf = std::numeric_limits<ll>::max() / 4;
using Matrix = std::vector<std::vector<ll>>;
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
int range(int lo, int hi) { return lo + int(next() % std::uint32_t(hi - lo + 1)); }
};
long long relaxations = 0;
ll fastTour(const Matrix& t) { // the bitmask DP, without the route
int n = int(t.size());
std::vector<std::vector<ll>> best(1 << n, std::vector<ll>(n, kInf));
best[1][0] = 0;
for (int mask = 1; mask < (1 << n); mask += 2) // odd masks contain the depot
for (int v = 0; v < n; ++v) {
if (best[mask][v] == kInf) continue;
for (int w = 0; w < n; ++w) {
if (mask & (1 << w)) continue;
++relaxations;
ll& slot = best[mask | (1 << w)][w];
slot = std::min(slot, best[mask][v] + t[v][w]);
}
}
ll answer = kInf;
for (int v = 1; v < n; ++v) answer = std::min(answer, best.back()[v] + t[v][0]);
return answer;
}
ll slowTour(const Matrix& t) { // try every order of the stops 1..n-1
int n = int(t.size());
std::vector<int> order(n - 1);
std::iota(order.begin(), order.end(), 1);
ll answer = kInf;
do {
ll cost = t[0][order[0]] + t[order.back()][0];
for (int i = 0; i + 1 < n - 1; ++i) cost += t[order[i]][order[i + 1]];
answer = std::min(answer, cost);
} while (std::ranges::next_permutation(order).found);
return answer;
}
int main() {
Lcg rng{2024};
int checked = 0;
for (int test = 0; test < 3000; ++test) {
int n = rng.range(2, 7);
int maxW = test % 3 == 0 ? 1 : 20; // many ties when weights are 0 or 1
Matrix t(n, std::vector<ll>(n));
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j) t[i][j] = i == j ? 0 : rng.range(0, maxW);
if (fastTour(t) != slowTour(t)) {
std::println("mismatch on test {} (n = {})", test, n);
return 1;
}
++checked;
}
std::println("{} random cases, all agree", checked);
std::println("{:>3} {:>12} {:>12} {:>15}", "n", "relaxations", "n^2 2^n / 8",
"orders");
long long orders = 1; // (n-1)!
for (int n = 2; n <= 16; ++n) {
orders *= n - 1;
if (n % 2 || n < 4) continue;
relaxations = 0;
fastTour(Matrix(n, std::vector<ll>(n, 1)));
std::println("{:>3} {:>12} {:>12} {:>15}", n, relaxations,
(1LL * n * n << n) / 8, orders);
}
}Press Run to compile and run this example.
3000 random cases, all agree
n relaxations n^2 2^n / 8 orders
4 15 32 6
6 165 288 120
8 1351 2048 5040
10 9225 12800 362880
12 56331 73728 39916800
14 319501 401408 6227020800
16 1720335 2097152 1307674368000The relaxation count tracks n^2 2^n / 8 and creeps towards it as n grows (the factor 8 comes from only extending from stops inside the mask to stops outside it, roughly n/2 each, over half of all masks). Each step of 2 in n multiplies the DP's work by about 5 to 11 in this range (the factor shrinks towards 4 as n grows), while the brute force multiplies by about n^2. At n = 8 the brute force is still competitive; at n = 16 it does about 760,000 times more work. That crossover is the constraints-first lesson in one table: n <= 8 says "permutations are fine", n between about 15 and 20 says "bitmask DP".
Assignment: let `popcount` be part of the state
Sometimes the second index is redundant. Suppose n couriers and n parcels, cost[i][j] minutes for courier i to deliver parcel j, each courier takes exactly one parcel, and we want the least total. If couriers are handed parcels in the fixed order 0, 1, 2, ..., then after k hand-outs the couriers 0..k-1 are done, and k is simply the number of bits in the mask of parcels taken. So the state is just best[mask], the next courier is std::popcount(mask), and the whole DP is O(n 2^n).
std::popcount (C++20, header <bit>) accepts only unsigned integer types, which is why the mask below is unsigned. The sample input is three couriers:
3
9 2 7
6 4 3
5 8 1#include <algorithm>
#include <bit>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <print>
#include <sstream>
#include <vector>
using ll = long long;
constexpr ll kInf = std::numeric_limits<ll>::max() / 4;
using Matrix = std::vector<std::vector<ll>>;
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
int range(int lo, int hi) { return lo + int(next() % std::uint32_t(hi - lo + 1)); }
};
ll assignFast(const Matrix& cost) {
int n = int(cost.size());
// best[mask]: couriers 0..popcount(mask)-1 have taken exactly the parcels in mask
std::vector<ll> best(1u << n, kInf);
best[0] = 0;
for (unsigned mask = 0; mask + 1 < (1u << n); ++mask) {
if (best[mask] == kInf) continue;
int courier = std::popcount(mask); // the next courier to give a parcel to
for (int p = 0; p < n; ++p)
if (!(mask >> p & 1))
best[mask | 1u << p] =
std::min(best[mask | 1u << p], best[mask] + cost[courier][p]);
}
return best.back();
}
ll assignSlow(const Matrix& cost) { // every permutation: courier i takes perm[i]
std::vector<int> perm(cost.size());
std::iota(perm.begin(), perm.end(), 0);
ll answer = kInf;
do {
ll total = 0;
for (std::size_t i = 0; i < perm.size(); ++i) total += cost[i][perm[i]];
answer = std::min(answer, total);
} while (std::ranges::next_permutation(perm).found);
return answer;
}
void solve(std::istream& in, std::ostream& out) {
int n;
in >> n;
Matrix cost(n, std::vector<ll>(n));
for (auto& row : cost)
for (ll& x : row) in >> x;
out << assignFast(cost) << '\n';
}
int main() {
std::istringstream in("3\n9 2 7\n6 4 3\n5 8 1\n");
solve(in, std::cout);
Lcg rng{7};
int checked = 0;
for (; checked < 2000; ++checked) {
int n = rng.range(1, 7);
Matrix cost(n, std::vector<ll>(n));
for (auto& row : cost)
for (ll& x : row) x = rng.range(0, checked % 2 ? 3 : 1000);
if (assignFast(cost) != assignSlow(cost)) {
std::println("mismatch at case {}", checked);
return 1;
}
}
std::println("{} random cases, all agree", checked);
}Press Run to compile and run this example.
9
2000 random cases, all agreeThe optimum is courier 0 on parcel 1, courier 1 on parcel 0 and courier 2 on parcel 2: 2 + 6 + 1 = 9. Assignment also has a polynomial algorithm, the Hungarian method (Kuhn, 1955; its usual modern form runs in O(n^3)), so for n = 300 you would need that instead. The bitmask version wins on simplicity and absorbs side conditions that the specialized algorithm cannot.
Iterating over submasks costs `3^n`, not `4^n`
Some transitions remove a whole group at once. A workshop has n jobs with durations dur[i] and a daily capacity cap; jobs cannot be split, and we want the fewest days. A day is any group of jobs whose total fits, so
days[mask] = 1 + min over nonempty submasks s of mask with load[s] <= cap of days[mask ^ s].
Looping s over all 2^n integers for every mask would cost 4^n. The loop for (s = mask; s; s = (s - 1) & mask) visits exactly the nonempty submasks of mask, in decreasing order: subtracting 1 clears the lowest set bit and sets all bits below it, and & mask throws away the bits that are not in mask. Summed over all masks it visits 3^n - 2^n pairs, because each item independently is outside mask, in mask but not s, or in both. The program counts that, and also builds the load of every mask in O(2^n) from the mask without its lowest job, using std::countr_zero to find that job.
#include <algorithm>
#include <bit>
#include <cstdint>
#include <print>
#include <vector>
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
int range(int lo, int hi) { return lo + int(next() % std::uint32_t(hi - lo + 1)); }
};
long long submaskSteps = 0;
// Fewest days to run all jobs if one day holds total duration at most cap.
// Every job must satisfy duration <= cap.
int fewestDays(const std::vector<int>& dur, int cap) {
int n = int(dur.size());
std::vector<int> load(1u << n, 0);
for (unsigned m = 1; m < load.size(); ++m) // load of m = load without its lowest job
load[m] = load[m & (m - 1)] + dur[std::countr_zero(m)];
std::vector<int> days(1u << n, n + 1);
days[0] = 0;
for (unsigned mask = 1; mask < days.size(); ++mask)
for (unsigned s = mask; s; s = (s - 1) & mask) { // nonempty submasks
++submaskSteps;
if (load[s] <= cap) days[mask] = std::min(days[mask], days[mask ^ s] + 1);
}
return days.back();
}
int slowDays(const std::vector<int>& dur, int cap) { // try every day label per job
int n = int(dur.size()), best = n;
std::vector<int> label(n, 0);
while (true) {
std::vector<int> sum(n, 0);
for (int i = 0; i < n; ++i) sum[label[i]] += dur[i];
if (std::ranges::all_of(sum, [&](int x) { return x <= cap; }))
best = std::min(best, n - int(std::ranges::count(sum, 0))); // days used
int i = 0; // next labelling, like counting in base n
while (i < n && ++label[i] == n) label[i++] = 0;
if (i == n) return best;
}
}
int main() {
std::println("jobs 4 7 2 5 3 6, cap 10: {} days", fewestDays({4, 7, 2, 5, 3, 6}, 10));
Lcg rng{99};
int checked = 0;
for (; checked < 2000; ++checked) {
int n = rng.range(1, 5), cap = rng.range(1, 12);
std::vector<int> dur(n);
for (int& d : dur) d = rng.range(1, cap);
if (fewestDays(dur, cap) != slowDays(dur, cap)) {
std::println("mismatch at case {}", checked);
return 1;
}
}
std::println("{} random cases, all agree", checked);
std::println("{:>3} {:>10} {:>10}", "n", "steps", "3^n");
for (int n = 4; n <= 13; n += 3) {
submaskSteps = 0;
fewestDays(std::vector<int>(n, 1), n);
long long pow3 = 1;
for (int i = 0; i < n; ++i) pow3 *= 3;
std::println("{:>3} {:>10} {:>10}", n, submaskSteps, pow3);
}
}Press Run to compile and run this example.
jobs 4 7 2 5 3 6, cap 10: 3 days
2000 random cases, all agree
n steps 3^n
4 65 81
7 2059 2187
10 58025 59049
13 1586131 1594323The total duration is 27, so at least 3 days are needed, and {4, 6}, {7, 3}, {2, 5} achieves it. The counts are exactly 3^n - 2^n. In practice 3^n is fine up to about n = 16 (43 million steps); at n = 20 it is 3.5 billion, so look for a cheaper state instead.
Interval DP: a range that splits into two ranges
The tape problem: n recordings lie side by side with sizes s[0..n-1]. The splicer joins two neighboring recordings into one, and a join costs the size of the result. Join everything into a single recording at the least total cost. (Merging adjacent piles, choosing the parenthesization of a matrix product and building an optimal binary search tree share this structure.)
Look at the last join. It combines some left block l..k with the right block k+1..r, and it costs sum(l..r) no matter how the two blocks were built. The two blocks were built independently, so:
cost[l][r] = sum(l..r) + min over l <= k < r of cost[l][k] + cost[k+1][r], with cost[i][i] = 0.
There are O(n^2) intervals and O(n) split points each, so O(n^3) time, which the constraints table puts at n up to about 500. The trap is the loop order. cost[l][r] reads cost[k+1][r], an interval that starts later than l. If you loop l upwards from 0, that value is still 0 when you read it. Loop by interval length instead: every interval you read is shorter, so it is already final. The program runs both orders against a brute force that tries every sequence of adjacent joins.
#include <algorithm>
#include <cstdint>
#include <limits>
#include <print>
#include <vector>
using ll = long long;
constexpr ll kInf = std::numeric_limits<ll>::max() / 4;
struct Lcg {
std::uint64_t s;
std::uint32_t next() {
s = s * 6364136223846793005ULL + 1442695040888963407ULL;
return std::uint32_t(s >> 33);
}
int range(int lo, int hi) { return lo + int(next() % std::uint32_t(hi - lo + 1)); }
};
// cost[l][r]: cheapest way to splice recordings l..r into one.
ll spliceCost(const std::vector<ll>& s, bool byLength) {
int n = int(s.size());
std::vector<ll> pre(n + 1, 0);
for (int i = 0; i < n; ++i) pre[i + 1] = pre[i] + s[i];
std::vector<std::vector<ll>> cost(n, std::vector<ll>(n, 0));
auto relax = [&](int l, int r) {
cost[l][r] = kInf;
for (int k = l; k < r; ++k) // the last splice joins l..k with k+1..r
cost[l][r] = std::min(cost[l][r], cost[l][k] + cost[k + 1][r]);
cost[l][r] += pre[r + 1] - pre[l];
};
if (byLength) { // shorter intervals first: everything we read is final
for (int len = 2; len <= n; ++len)
for (int l = 0; l + len - 1 < n; ++l) relax(l, l + len - 1);
} else { // the tempting order: reads cost[k+1][r] before it is computed
for (int l = 0; l < n; ++l)
for (int r = l + 1; r < n; ++r) relax(l, r);
}
return cost[0][n - 1];
}
ll slowSplice(std::vector<ll> s) { // try every sequence of adjacent splices
if (s.size() == 1) return 0;
ll best = kInf;
for (std::size_t i = 0; i + 1 < s.size(); ++i) {
std::vector<ll> t = s;
t[i] += t[i + 1];
t.erase(t.begin() + i + 1);
best = std::min(best, t[i] + slowSplice(t));
}
return best;
}
int main() {
std::println("sizes 4 1 3 2: {}", spliceCost({4, 1, 3, 2}, true));
for (bool byLength : {false, true}) {
Lcg rng{31};
int checked = 0;
for (; checked < 2000; ++checked) {
std::vector<ll> s(rng.range(1, 7));
for (ll& x : s) x = rng.range(1, 9);
if (spliceCost(s, byLength) != slowSplice(s)) {
std::print("{}: case {} fails, sizes", byLength ? "by length" : "by l, r",
checked);
for (ll x : s) std::print(" {}", x);
std::println(": got {}, want {}", spliceCost(s, byLength), slowSplice(s));
break;
}
}
if (checked == 2000)
std::println("by length: {} random cases, all agree", checked);
}
}Press Run to compile and run this example.
sizes 4 1 3 2: 20
by l, r: case 0 fails, sizes 7 4 6 1: got 18, want 36
by length: 2000 random cases, all agreeFor 4 1 3 2, joining 4 + 1 (cost 5) and 3 + 2 (cost 5) and then the two halves (cost 10) gives 20; any order that first joins 1 + 3 costs more. The wrong loop order does not crash: it produces a plausible, too-small number, 18 instead of 36, because it adds zeros for intervals it has not solved yet. Only a comparison with a slow solution exposes it.