Problem 212: Combined Volume of Cuboids

Mathematical Development

Definition 1. Let $\mathcal{C} = \{C_1, \ldots, C_n\}$ be a collection of axis-aligned cuboids in $\mathbb{R}^3$. The combined volume is $\operatorname{vol}\!\bigl(\bigcup_{i=1}^{n} C_i\bigr)$, where $\operatorname{vol}$ denotes 3-dimensional Lebesgue measure.

Theorem 1 (Cavalieri’s principle for cuboid unions). Let $z_1 < z_2 < \cdots < z_M$ be the distinct $z$-coordinates among all faces $\{z_0^{(i)}, z_0^{(i)} + \delta z^{(i)}\}_{i=1}^{n}$. Then

where $A_j$ is the 2-dimensional Lebesgue measure of the union of the $xy$-projections of all cuboids active in the slab $[z_j, z_{j+1})$, i.e., those $C_i$ with $z_0^{(i)} \leq z_j$ and $z_{j+1} \leq z_0^{(i)} + \delta z^{(i)}$.

Proof. By Cavalieri’s principle, $\operatorname{vol}(U) = \int_{-\infty}^{\infty} \operatorname{Area}(U \cap \{z = t\})\,dt$. For $t \in (z_j, z_{j+1})$, the set of cuboids containing the plane $z = t$ is precisely the set of cuboids active in the slab $[z_j, z_{j+1})$, since no cuboid face lies strictly between $z_j$ and $z_{j+1}$. Hence the cross-sectional area $\operatorname{Area}(U \cap \{z = t\})$ is constant on each open interval $(z_j, z_{j+1})$, equal to $A_j$. Integrating over each slab gives $(z_{j+1} - z_j) \cdot A_j$, and summing over all slabs yields the total volume. $\square$

Theorem 2 (2D area via $x$-sweep). For a fixed $z$-slab, let $\mathcal{R}$ be the set of active $xy$-rectangles. Let $x_1 < x_2 < \cdots < x_K$ be the distinct $x$-coordinates among all left and right edges of rectangles in $\mathcal{R}$. Then

where $L_\ell$ is the total length of the union of $y$-intervals from rectangles in $\mathcal{R}$ that cover the strip $[x_\ell, x_{\ell+1})$.

Proof. Apply Cavalieri’s principle in dimension 2: $A_j = \int_{-\infty}^{\infty} \operatorname{Length}(\text{section at } x)\,dx$. Between consecutive $x$-boundaries, the set of active rectangles is constant, so the $y$-union length is constant on each strip. $\square$

Lemma 1 (Interval union via greedy merge). Given intervals $[a_1, b_1], \ldots, [a_m, b_m]$ sorted by left endpoint ($a_1 \leq a_2 \leq \cdots \leq a_m$), the union length is computed by the following algorithm in $\Theta(m)$ time:

1. Initialize $(c, d) \leftarrow (a_1, b_1)$ and $\lambda \leftarrow 0$.
2. For $k = 2, \ldots, m$: if $a_k \geq d$, set $\lambda \leftarrow \lambda + (d - c)$ and $(c, d) \leftarrow (a_k, b_k)$; otherwise set $d \leftarrow \max(d, b_k)$.
3. Set $\lambda \leftarrow \lambda + (d - c)$.

Then $\lambda = \operatorname{Length}\!\bigl(\bigcup_{k=1}^{m} [a_k, b_k]\bigr)$.

Proof. The algorithm maintains the invariant that $[c, d]$ is the current maximal merged interval. Since the intervals are sorted by left endpoint, any new interval $[a_k, b_k]$ either extends the current interval (if $a_k < d$) or starts a disjoint component (if $a_k \geq d$). In the latter case, $[c, d]$ is closed and its length added to $\lambda$. The final addition accounts for the last open interval. Correctness follows by induction on $m$. $\square$

Editorial

We use dynamic programming over the state space implied by the derivation, apply each admissible transition, and read the answer from the final table entry.

Pseudocode

    Z = sorted distinct set of {z0, z0 + dz} for all cuboids
    V = 0

    For j from 1 to |Z| - 1:
        dz = Z[j+1] - Z[j]
        active = {c in cuboids : c.z0 <= Z[j] and Z[j+1] <= c.z0 + c.dz}
        If active is empty then continue

        X = sorted distinct set of {c.x0, c.x0 + c.dx} for c in active
        area = 0

        For l from 1 to |X| - 1:
            dx = X[l+1] - X[l]
            intervals = [(c.y0, c.y0 + c.dy) : c in active, c.x0 <= X[l], X[l+1] <= c.x0 + c.dx]
            If intervals is empty then continue
            sort intervals by left endpoint
            Ly = GreedyMergeLength(intervals)
            area += dx * Ly

        V += dz * area

    Return V

Complexity Analysis

Let $n = 50{,}000$ denote the number of cuboids.

Time. There are $O(n)$ distinct $z$-values, yielding $O(n)$ slabs. For each slab, identifying the active set costs $O(n)$. The $x$-sweep processes $O(n)$ strips, and for each strip the $y$-interval merge costs $O(k \log k)$ where $k$ is the number of active rectangles in that strip. Worst-case total: $O(n^2 \cdot k \log k)$. In practice, the pseudorandom distribution over a $10{,}000^3$ grid produces sparse overlaps, yielding average-case $O(n^2)$ behavior.

Space. $O(n)$ for cuboid storage and working arrays.

Answer

#include <bits/stdc++.h>
using namespace std;

int main() {
    const int NCUBOIDS = 50000;
    const int SLEN = 6 * NCUBOIDS;

    vector<long long> S(SLEN + 1);
    for (int k = 1; k <= 55; k++) {
        long long k3 = (long long)k * k * k;
        S[k] = ((100003LL - 200003LL * k + 300007LL * k3) % 1000000LL + 1000000LL) % 1000000LL;
    }
    for (int k = 56; k <= SLEN; k++)
        S[k] = (S[k - 24] + S[k - 55]) % 1000000LL;

    struct Cuboid { int x1, y1, z1, x2, y2, z2; };
    vector<Cuboid> cuboids(NCUBOIDS);

    for (int n = 1; n <= NCUBOIDS; n++) {
        int i = n - 1;
        cuboids[i] = {
            (int)(S[6*n-5] % 10000), (int)(S[6*n-4] % 10000), (int)(S[6*n-3] % 10000),
            (int)(S[6*n-5] % 10000 + 1 + S[6*n-2] % 399),
            (int)(S[6*n-4] % 10000 + 1 + S[6*n-1] % 399),
            (int)(S[6*n-3] % 10000 + 1 + S[6*n] % 399)
        };
    }

    vector<int> zcoords;
    for (int i = 0; i < NCUBOIDS; i++) {
        zcoords.push_back(cuboids[i].z1);
        zcoords.push_back(cuboids[i].z2);
    }
    sort(zcoords.begin(), zcoords.end());
    zcoords.erase(unique(zcoords.begin(), zcoords.end()), zcoords.end());

    long long totalVolume = 0;

    for (int zi = 0; zi + 1 < (int)zcoords.size(); zi++) {
        int z0 = zcoords[zi], z1 = zcoords[zi + 1];
        int dz = z1 - z0;

        vector<int> active;
        for (int i = 0; i < NCUBOIDS; i++)
            if (cuboids[i].z1 <= z0 && z1 <= cuboids[i].z2)
                active.push_back(i);
        if (active.empty()) continue;

        vector<int> localx;
        for (int i : active) {
            localx.push_back(cuboids[i].x1);
            localx.push_back(cuboids[i].x2);
        }
        sort(localx.begin(), localx.end());
        localx.erase(unique(localx.begin(), localx.end()), localx.end());

        for (int xi = 0; xi + 1 < (int)localx.size(); xi++) {
            int x0 = localx[xi], x1 = localx[xi + 1];
            int dx = x1 - x0;

            vector<pair<int, int>> yintervals;
            for (int i : active)
                if (cuboids[i].x1 <= x0 && x1 <= cuboids[i].x2)
                    yintervals.push_back({cuboids[i].y1, cuboids[i].y2});
            if (yintervals.empty()) continue;

            sort(yintervals.begin(), yintervals.end());
            long long ylen = 0;
            int cy0 = yintervals[0].first, cy1 = yintervals[0].second;
            for (int j = 1; j < (int)yintervals.size(); j++) {
                if (yintervals[j].first >= cy1) {
                    ylen += cy1 - cy0;
                    cy0 = yintervals[j].first;
                    cy1 = yintervals[j].second;
                } else {
                    cy1 = max(cy1, yintervals[j].second);
                }
            }
            ylen += cy1 - cy0;
            totalVolume += (long long)dz * dx * ylen;
        }
    }

    cout << totalVolume << endl;
    return 0;
}

def solve():
    NCUBOIDS = 50000
    SLEN = 6 * NCUBOIDS

    S = [0] * (SLEN + 1)
    for k in range(1, 56):
        S[k] = (100003 - 200003 * k + 300007 * k * k * k) % 1000000
    for k in range(56, SLEN + 1):
        S[k] = (S[k - 24] + S[k - 55]) % 1000000

    cuboids = []
    for n in range(1, NCUBOIDS + 1):
        x0 = S[6 * n - 5] % 10000
        y0 = S[6 * n - 4] % 10000
        z0 = S[6 * n - 3] % 10000
        dx = 1 + S[6 * n - 2] % 399
        dy = 1 + S[6 * n - 1] % 399
        dz = 1 + S[6 * n] % 399
        cuboids.append((x0, y0, z0, x0 + dx, y0 + dy, z0 + dz))

    zset = set()
    for c in cuboids:
        zset.add(c[2])
        zset.add(c[5])
    zcoords = sorted(zset)

    total_volume = 0

    for zi in range(len(zcoords) - 1):
        z0, z1 = zcoords[zi], zcoords[zi + 1]
        dz = z1 - z0
        active = [c for c in cuboids if c[2] <= z0 and z1 <= c[5]]
        if not active:
            continue

        xset = set()
        for c in active:
            xset.add(c[0])
            xset.add(c[3])
        xcoords = sorted(xset)

        for xi in range(len(xcoords) - 1):
            x0, x1 = xcoords[xi], xcoords[xi + 1]
            dx = x1 - x0
            yintervals = []
            for c in active:
                if c[0] <= x0 and x1 <= c[3]:
                    yintervals.append((c[1], c[4]))
            if not yintervals:
                continue

            yintervals.sort()
            ylen = 0
            cy0, cy1 = yintervals[0]
            for j in range(1, len(yintervals)):
                if yintervals[j][0] >= cy1:
                    ylen += cy1 - cy0
                    cy0, cy1 = yintervals[j]
                else:
                    cy1 = max(cy1, yintervals[j][1])
            ylen += cy1 - cy0

            total_volume += dz * dx * ylen

    print(total_volume)

if __name__ == "__main__":
    solve()