Project Euler

Combined Volume of Cuboids

An axis-aligned cuboid with parameters {(x_0, y_0, z_0), (delta x, delta y, delta z)} consists of all points (X, Y, Z) with x_0 <= X <= x_0 + delta x, y_0 <= Y <= y_0 + delta y, z_0 <= Z <= z_0 + d...

Source sync May 21, 2026

Problem #0212

Level Level 12

Solved By 1,635

Languages C++, Python

Answer 328968937309

Length 515 words

modular_arithmeticgeometrylinear_algebra

An axis-aligned cuboid, specified by parameters $\{(x_0, y_0, z_0), (dx, dy, dz)\}$, consists of all points $(X,Y,Z)$ such that $x_0 \le X \le x_0 + dx$, $y_0 \le Y \le y_0 + dy$ and $z_0 \le Z \le z_0 + dz$. The volume of the cuboid is the product, $dx \times dy \times dz$. The combined volume of a collection of cuboids is the volume of their union and will be less than the sum of the individual volumes if any cuboids overlap.

Let $C_1, \dots, C_{50000}$ be a collection of $50000$ axis-aligned cuboids such that $C_n$ has parameters

\begin{align*} x_0 &= S_{6n - 5} \bmod 10000\\ y_0 &= S_{6n - 4} \bmod 10000\\ z_0 &= S_{6n - 3} \bmod 10000\\ dx &= 1 + (S_{6n - 2} \bmod 399)\\ dy &= 1 + (S_{6n - 1} \bmod 399)\\ dz &= 1 + (S_{6n} \bmod 399) \end{align*}

where $S_1,\dots,S_{300000}$ come from the "Lagged Fibonacci Generator":

For $1 \le k \le 55$, $S_k = [100003 - 200003k + 300007k^3] \pmod{1000000}$.
For $56 \le k$, $S_k = [S_{k -24} + S_{k - 55}] \pmod{1000000}$.

Thus, $C_1$ has parameters $\{(7,53,183),(94,369,56)\}$,

$C_2$ has parameters $\{(2383,3563,5079),(42,212,344)\}$, and so on.

The combined volume of the first $100$ cuboids, $C_1, \dots, C_{100}$, is $723581599$.

What is the combined volume of all $50000$ cuboids, $C_1, \dots, C_{50000}$?

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

\operatorname{vol}\!\left(\bigcup_{i=1}^{n} C_i\right) = \sum_{j=1}^{M-1} (z_{j+1} - z_j) \cdot A_j

where $A_j$ is the 2-dimensional measure of the union of the $xy$ -projections of all cuboids active in the slab $[z_j, z_{j+1})$ .

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})$ , no cuboid face lies strictly between $z_j$ and $z_{j+1}$ , so the active set is constant there. Hence the cross-sectional area is constant on that slab, equal to $A_j$ , and integrating gives the stated sum. $\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

A_j = \sum_{\ell=1}^{K-1} (x_{\ell+1} - x_\ell) \cdot L_\ell

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

Proof. Apply Cavalieri’s principle again in dimension 2. Between consecutive $x$ -boundaries the active rectangle set is fixed, so the union length in the $y$ -direction is constant on that strip. $\square$

Lemma 1 (Interval union via greedy merge). Given intervals $[a_1, b_1], \ldots, [a_m, b_m]$ sorted by left endpoint, their union length can be computed in linear time by greedily merging overlaps.

Proof. Because the intervals are sorted by left endpoint, each new interval either starts a new disjoint component or extends the current merged component. Summing the lengths of the resulting maximal merged components gives the union length. $\square$

Lemma 2 (Integer layer reduction). In this problem all face coordinates are integers. Therefore the volume is also the sum over the unit slabs $[z, z+1)$ of the union area of the cuboids active at height $z$ .

Proof. Every cuboid is of the form

[x_1, x_2] \times [y_1, y_2] \times [z_1, z_2]

with integer endpoints, so its contribution is constant on each unit slab $[z, z+1)$ with $z_1 \leq z < z_2$ . Summing those constant cross-sectional areas over all integer $z$ reproduces the volume. $\square$

Editorial

The naive decomposition by distinct $z$ -faces is mathematically correct but still too slow if every slab rescans all 50,000 cuboids and then recomputes the 2D union with nested brute force. The crucial simplification here is that every coordinate is an integer below 10,400, so the volume can be accumulated one unit $z$ -layer at a time.

For each integer layer $[z, z+1)$ we maintain the active cuboids by start and end events. Their projections onto the $xy$ -plane are just rectangles, so the remaining task is a standard rectangle-union problem: sweep in $x$ , and keep a segment tree over the $y$ -axis storing the total covered $y$ -length. That gives the union area of the active rectangles for the current layer, which contributes directly to the total volume because the slab thickness is 1.

Pseudocode

Generate the lagged sequence and build all cuboids as half-open boxes
[x1, x2) x [y1, y2) x [z1, z2).

For each cuboid:
    append its index to starts[z1]
    append its index to ends[z2]

active = empty set
volume = 0

For z from 0 up to the largest top face minus 1:
    Remove from active every cuboid whose z2 equals z.
    Add to active every cuboid whose z1 equals z.

    If active is empty:
        continue

    Build x-events from the active cuboids:
        (x1, +1, y1, y2) and (x2, -1, y1, y2)
    Sort the x-events by x.

    Clear the segment tree over the y-axis.
    area = 0
    prev_x = first event position

    Sweep through the x-events in order:
        area += covered_y_length * (current_x - prev_x)
        apply every event at current_x to the segment tree
        prev_x = current_x

    volume += area

Return volume

Complexity Analysis

Let $A_z$ be the number of cuboids active in the unit slab $[z, z+1)$ . Then:

Time: building the event buckets is $O(n)$ . For each $z$ -layer, the 2D sweep costs $O(A_z \log A_z + A_z \log 10400)$ . The total is $O\!\left(\sum_z A_z \log A_z\right),$ and here $\sum_z A_z = \sum_i \delta z_i \leq 50{,}000 \cdot 399$ .
Space: $O(n)$ for the cuboids, event buckets, active set, and sweep events, plus $O(10400)$ for the segment tree.

Answer

\boxed{328968937309}

C++ project_euler/problem_212/solution.cpp

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

namespace {

struct Cuboid {
    int x1, y1, z1;
    int x2, y2, z2;
};

struct Event {
    int x;
    int delta;
    int y1;
    int y2;

    bool operator<(const Event& other) const {
        return x < other.x;
    }
};

class SegmentTree {
public:
    explicit SegmentTree(int size) : size_(size), cover_(4 * size, 0), length_(4 * size, 0) {}

    void reset() {
        fill(cover_.begin(), cover_.end(), 0);
        fill(length_.begin(), length_.end(), 0);
    }

    void addInterval(int left, int right, int delta) {
        update(1, 0, size_, left, right, delta);
    }

    int coveredLength() const {
        return length_[1];
    }

private:
    int size_;
    vector<int> cover_;
    vector<int> length_;

    void update(int node, int left, int right, int queryLeft, int queryRight, int delta) {
        if (queryRight <= left || right <= queryLeft) {
            return;
        }

        if (queryLeft <= left && right <= queryRight) {
            cover_[node] += delta;
        } else {
            int mid = (left + right) / 2;
            update(node * 2, left, mid, queryLeft, queryRight, delta);
            update(node * 2 + 1, mid, right, queryLeft, queryRight, delta);
        }

        if (cover_[node] > 0) {
            length_[node] = right - left;
        } else if (right - left == 1) {
            length_[node] = 0;
        } else {
            length_[node] = length_[node * 2] + length_[node * 2 + 1];
        }
    }
};

long long computeUnionArea(const vector<int>& active, const vector<Cuboid>& cuboids, SegmentTree& tree) {
    if (active.empty()) {
        return 0;
    }

    vector<Event> events;
    events.reserve(active.size() * 2);
    for (int idx : active) {
        const Cuboid& cuboid = cuboids[idx];
        events.push_back({cuboid.x1, 1, cuboid.y1, cuboid.y2});
        events.push_back({cuboid.x2, -1, cuboid.y1, cuboid.y2});
    }

    sort(events.begin(), events.end());
    tree.reset();

    long long area = 0;
    int prevX = events.front().x;
    int i = 0;

    while (i < static_cast<int>(events.size())) {
        int x = events[i].x;
        area += 1LL * tree.coveredLength() * (x - prevX);

        while (i < static_cast<int>(events.size()) && events[i].x == x) {
            tree.addInterval(events[i].y1, events[i].y2, events[i].delta);
            i++;
        }

        prevX = x;
    }

    return area;
}

}  // namespace

int main() {
    const int cuboidCount = 50'000;
    const int sequenceLength = 6 * cuboidCount;
    const int maxCoord = 10'400;

    vector<int> sequence(sequenceLength + 1, 0);
    for (int k = 1; k <= 55; k++) {
        long long k3 = 1LL * k * k * k;
        sequence[k] = static_cast<int>((100003LL - 200003LL * k + 300007LL * k3) % 1'000'000LL);
        if (sequence[k] < 0) {
            sequence[k] += 1'000'000;
        }
    }
    for (int k = 56; k <= sequenceLength; k++) {
        sequence[k] = (sequence[k - 24] + sequence[k - 55]) % 1'000'000;
    }

    vector<Cuboid> cuboids(cuboidCount);
    int maxZ = 0;
    for (int n = 1; n <= cuboidCount; n++) {
        int x1 = sequence[6 * n - 5] % 10'000;
        int y1 = sequence[6 * n - 4] % 10'000;
        int z1 = sequence[6 * n - 3] % 10'000;
        int x2 = x1 + 1 + sequence[6 * n - 2] % 399;
        int y2 = y1 + 1 + sequence[6 * n - 1] % 399;
        int z2 = z1 + 1 + sequence[6 * n] % 399;
        cuboids[n - 1] = {x1, y1, z1, x2, y2, z2};
        maxZ = max(maxZ, z2);
    }

    vector<vector<int>> starts(maxZ + 1), ends(maxZ + 1);
    for (int i = 0; i < cuboidCount; i++) {
        starts[cuboids[i].z1].push_back(i);
        ends[cuboids[i].z2].push_back(i);
    }

    vector<int> active;
    active.reserve(cuboidCount);
    vector<int> position(cuboidCount, -1);
    SegmentTree tree(maxCoord);

    long long totalVolume = 0;

    for (int z = 0; z < maxZ; z++) {
        for (int idx : ends[z]) {
            int pos = position[idx];
            if (pos == -1) {
                continue;
            }
            int last = active.back();
            active[pos] = last;
            position[last] = pos;
            active.pop_back();
            position[idx] = -1;
        }

        for (int idx : starts[z]) {
            position[idx] = static_cast<int>(active.size());
            active.push_back(idx);
        }

        if (!active.empty()) {
            totalVolume += computeUnionArea(active, cuboids, tree);
        }
    }

    cout << totalVolume << '\n';
    return 0;
}

Python project_euler/problem_212/solution.py

def solve():
    cuboid_count = 50_000
    sequence_length = 6 * cuboid_count

    sequence = [0] * (sequence_length + 1)
    for k in range(1, 56):
        sequence[k] = (100003 - 200003 * k + 300007 * k * k * k) % 1_000_000
    for k in range(56, sequence_length + 1):
        sequence[k] = (sequence[k - 24] + sequence[k - 55]) % 1_000_000

    cuboids = []
    max_z = 0
    for n in range(1, cuboid_count + 1):
        x1 = sequence[6 * n - 5] % 10_000
        y1 = sequence[6 * n - 4] % 10_000
        z1 = sequence[6 * n - 3] % 10_000
        x2 = x1 + 1 + sequence[6 * n - 2] % 399
        y2 = y1 + 1 + sequence[6 * n - 1] % 399
        z2 = z1 + 1 + sequence[6 * n] % 399
        cuboids.append((x1, y1, z1, x2, y2, z2))
        if z2 > max_z:
            max_z = z2

    starts = [[] for _ in range(max_z + 1)]
    ends = [[] for _ in range(max_z + 1)]
    for idx, (_, _, z1, _, _, z2) in enumerate(cuboids):
        starts[z1].append(idx)
        ends[z2].append(idx)

    max_coord = 10_400
    tree_size = 4 * max_coord

    def compute_union_area(active):
        if not active:
            return 0

        cover = [0] * tree_size
        covered_length = [0] * tree_size
        events = []

        for idx in active:
            x1, y1, _, x2, y2, _ = cuboids[idx]
            events.append((x1, 1, y1, y2))
            events.append((x2, -1, y1, y2))

        events.sort()

        def update(node, left, right, query_left, query_right, delta):
            if query_right <= left or right <= query_left:
                return

            if query_left <= left and right <= query_right:
                cover[node] += delta
            else:
                mid = (left + right) // 2
                update(node * 2, left, mid, query_left, query_right, delta)
                update(node * 2 + 1, mid, right, query_left, query_right, delta)

            if cover[node] > 0:
                covered_length[node] = right - left
            elif right - left == 1:
                covered_length[node] = 0
            else:
                covered_length[node] = covered_length[node * 2] + covered_length[node * 2 + 1]

        area = 0
        prev_x = events[0][0]
        index = 0

        while index < len(events):
            x = events[index][0]
            area += covered_length[1] * (x - prev_x)

            while index < len(events) and events[index][0] == x:
                _, delta, y1, y2 = events[index]
                update(1, 0, max_coord, y1, y2, delta)
                index += 1

            prev_x = x

        return area

    active = []
    position = [-1] * cuboid_count
    total_volume = 0

    for z in range(max_z):
        for idx in ends[z]:
            pos = position[idx]
            if pos != -1:
                last = active[-1]
                active[pos] = last
                position[last] = pos
                active.pop()
                position[idx] = -1

        for idx in starts[z]:
            position[idx] = len(active)
            active.append(idx)

        if active:
            total_volume += compute_union_area(active)

    print(total_volume)


if __name__ == "__main__":
    solve()