Problem 212: Combined Volume of Cuboids

Project Euler Problem 212 Solution

EulerSolve provides an optimized solution for Project Euler Problem 212, Combined Volume of Cuboids, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 212 asks for the volume of the union of $N=50000$ axis-aligned cuboids. The cuboids are generated deterministically from a lagged Fibonacci sequence, so the input is large but completely fixed. Each cuboid has the half-open form $$C_i=[x_i,x_i+d_{x,i})\times[y_i,y_i+d_{y,i})\times[z_i,z_i+d_{z,i}),$$ with $x_i,y_i,z_i\in\{0,\dots,9999\}$ and $1\le d_{x,i},d_{y,i},d_{z,i}\le 399$. The challenge is the overlap structure: $50000$ cuboids can intersect in extremely complicated ways, so direct inclusion-exclusion is useless. The successful approach is to exploit the arithmetic structure of the coordinates. Because all coordinates are integers and every $x$-length is at most $399$, the entire 3D problem can be reduced to about $10^4$ independent 2D union-area problems. Mathematical Approach From the sequence to the cuboids The pseudo-random sequence is $$S_k=(100003-200003k+300007k^3)\bmod 10^6\qquad(1\le k\le 55),$$ followed by the recurrence $$S_k=(S_{k-24}+S_{k-55})\bmod 10^6\qquad(k\ge 56).$$ Every cuboid uses six consecutive values: $$x_i=S_{6i-5}\bmod 10000,\qquad y_i=S_{6i-4}\bmod 10000,\qquad z_i=S_{6i-3}\bmod 10000,$$ $$d_{x,i}=1+(S_{6i-2}\bmod 399),\qquad d_{y,i}=1+(S_{6i-1}\bmod 399),\qquad d_{z,i}=1+(S_{6i}\bmod 399).$$ The half-open convention matters....

Detailed mathematical approach

Problem Summary

An axis-aligned cuboid with the parameters $\{(x_0,y_0,z_0),(dx,dy,dz)\}$ consists of all points $(X,Y,Z)$ with $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$; its volume is $dx\cdot dy\cdot dz$. The combined volume of several cuboids is the volume of their union. For $n=1,\dots,50000$, the cuboid $C_n$ has $x_0,y_0,z_0$ equal to $S_{6n-5},S_{6n-4},S_{6n-3}$ modulo $10000$ and $dx,dy,dz$ equal to $1$ plus the remainders of $S_{6n-2},S_{6n-1},S_{6n}$ modulo $399$, where the "Lagged Fibonacci Generator" gives

$$\begin{aligned}S_k&=\bigl(100003-200003k\\&\qquad+300007k^3\bigr)\bmod 10^6\end{aligned}$$

for $1\le k\le55$ and $S_k=(S_{k-24}+S_{k-55})\bmod 10^6$ for $k\ge56$. Thus $C_1$ has the parameters $\{(7,53,183),(94,369,56)\}$ and $C_2$ the parameters $\{(2383,3563,5079),(42,212,344)\}$, and the combined volume of $C_1,\dots,C_{100}$ is $723581599$. We must find the combined volume of all $50000$ cuboids. ^[1]

The whole solution in one sentence: all corners have integer coordinates, so space can be cut into slabs of thickness $1$ in which the union is a flat union of rectangles, whose area a sweeping line finds while a segment tree keeps the covered length on the line up to date; the areas of the slabs add up to the volume.

Mathematical Approach

We reduce the solid to slabs, a slab to a sweeping line and the line to a tree, following one small example throughout.

1. The first terms

The cubic gives $S_1=200007$, $S_2=100053$, $S_3=600183$, $S_4=500439$, $S_5=600863$ and $S_6=701497$; their remainders $7,53,183$ modulo $10000$ and $93,368,55$ modulo $399$ give $C_1$ as stated. From $k=56$ on, each term adds the terms $24$ and $55$ places back modulo $10^6$, as in an additive lagged Fibonacci generator with the classical lags $24$ and $55$. ^[2] The cubic is positive for every $k\ge1$, since $300007k^3\gt200003k$, so no negative remainders arise, but at $k=55$ it reaches $49902764463\gt2^{32}$ and needs 64-bit arithmetic. All coordinates are integers from $0$ to $9999+399=10398$; for our cuboids the largest far ends are $10388$ in $x$, $10387$ in $y$ and $10389$ in $z$.

2. Half-open boxes and unit cells

The unit cubes $\{(0,0,0),(1,1,1)\}$ and $\{(1,0,0),(1,1,1)\}$ touch along the face $X=1$ and have the combined volume $2$: a face has no volume. So we may replace every closed cuboid by the half-open box $[x_0,x_1)\times[y_0,y_1)\times[z_0,z_1)$ with $x_1=x_0+dx$, $y_1=y_0+dy$ and $z_1=z_0+dz$. A closed cuboid and its half-open box differ only by points of three faces, so the two unions differ by a set of volume zero, and touching boxes such as the two cubes become disjoint. With integer corners, every unit cell $[u,u+1)\times[v,v+1)\times[w,w+1)$ with integers $u,v,w$ lies inside a half-open box or outside it, because on each axis $[u,u+1)$ lies inside $[x_0,x_1)$ or misses it. The combined volume is therefore the number of covered unit cells. Counting them one by one, the brute force of the C++ tests, is hopeless here: the box $[0,10388)\times[0,10387)\times[0,10389)$ that holds all cuboids has about $1.12\cdot10^{12}$ cells.

3. Inclusion–exclusion and its limits

The running example consists of $E_1=\{(0,0,0),(2,2,3)\}$ and $E_2=\{(1,1,1),(2,2,3)\}$, the boxes $[0,2)\times[0,2)\times[0,3)$ and $[1,3)\times[1,3)\times[1,4)$ of volume $12$ each. Two boxes intersect in a box, bounded on each axis by the larger lower end and the smaller upper end, here $[1,2)\times[1,2)\times[1,3)$ of volume $2$, so inclusion–exclusion gives $12+12-2=22$. ^[3] For $n$ cuboids the formula has $2^n-1$ terms, but only subsets of pairwise overlapping cuboids contribute, and these share a box, since on each axis every lower end lies below every upper end. Among the first $100$ cuboids only $C_{59}$ and $C_{68}$ overlap, in $146\cdot140\cdot129=2636760$ cells; the $100$ volumes add up to $726218359$, and $726218359-2636760$ is the checkpoint $723581599$. Among all $50000$ cuboids, $76031$ pairs overlap and $169415$ subsets, of up to $9$ cuboids, share a box. Restricted to them, inclusion–exclusion would be feasible here, but its cost grows with the overlaps: $k$ cuboids around one point give $2^k-1$ terms. The volume of a union of boxes is Klee's measure problem, whose fastest known general method in three dimensions takes about $n^{3/2}$ steps; ^[4] we use the integer coordinates instead.

4. Cutting space into slabs

For an integer $x$, slab $x$ is the set of points with $x\le X\lt x+1$. A half-open box with integer ends $x_0\lt x_1$ contains the whole thickness of the slab when $x_0\le x\lt x_1$, and we then call it active, with the rectangle $[y_0,y_1)\times[z_0,z_1)$; otherwise it misses the slab. Inside slab $x$ the union is therefore a prism of height $1$ over the union $U_x$ of the active rectangles, and its volume is the area $A_x$ of $U_x$. The slabs $0,\dots,x_{\max}-1$, where $x_{\max}$ is the largest $x_1$, are disjoint and contain all cuboids, so

$$V=\sum_{x=0}^{x_{\max}-1}A_x ,$$

the simplest case of Cavalieri's principle. ^[5] In the example, slab $0$ holds only the rectangle $[0,2)\times[0,3)$ of $E_1$, of area $6$, and slab $2$ only the rectangle $[1,3)\times[1,4)$ of $E_2$, also of area $6$; slab $1$ holds both, which meet in $[1,2)\times[1,3)$, so $A_1=6+6-2=10$, and indeed $6+10+6=22$.

For the real input $x_{\max}=10388$. The cuboid $C_n$ is active in $dx$ slabs, so the slabs hold $M=\sum dx=10030243$ rectangles: about $966$ per slab on average, at most $1091$, and no slab is empty. Cutting only at the $10304$ distinct $x$-coordinates of faces, as one would for real coordinates, would save less than $1\%$ of the slabs.

5. The area of one slab by a sweep

Fix a slab with the rectangles $[a_j,b_j)\times[c_j,d_j)$ in the $(y,z)$-plane, $j=1,\dots,m$. The line of constant $y$ meets rectangle $j$ when $a_j\le y\lt b_j$, in the interval $[c_j,d_j)$ of the $z$-axis; the covered length $L(y)$ is the length of the union of these intervals. Let $y_1\lt\dots\lt y_p$ be the distinct values among all $a_j$ and $b_j$. No rectangle starts or ends strictly between neighbouring values, so $L$ is constant on each strip $y_i\le y\lt y_{i+1}$, and $L=0$ outside $[y_1,y_p)$. Hence

$$A=\sum_{i=1}^{p-1}L(y_i)\,(y_{i+1}-y_i).$$

This is a sweep line: it stops only at the events $+[c_j,d_j)$ at $y=a_j$ and $-[c_j,d_j)$ at $y=b_j$, ^[6] and all events at the same $y$ are applied, in any order, before $L$ is measured again. For unions of rectangles the method goes back to Bentley in 1977. ^[4] In slab $1$ of the example the events are $+[0,3)$ at $y=0$, $+[1,4)$ at $y=1$, $-[0,3)$ at $y=2$ and $-[1,4)$ at $y=3$; $L$ is $3$, $4$ and $3$ on the three strips, so again $A_1=10$.

6. The covered length in a segment tree

Recomputing $L$ from scratch at every event would take some $2\cdot10^{10}$ steps for the whole input. A segment tree on the range $[0,z_{\max})$, with $z_{\max}=10389$, needs only $O(\log z_{\max})$ steps per event. ^[7] Its root is $[0,z_{\max})$; a node $[l,r)$ with $r-l\ge2$ has the children $[l,h)$ and $[h,r)$, where $h=l+\lfloor(r-l)/2\rfloor$, and the leaves are unit intervals. To add or remove $[c,d)$, start at the root and at each node $[l,r)$ reached: stop if it misses $[c,d)$; if it lies inside $[c,d)$, add $\pm1$ to its counter $\kappa$ and stop; otherwise go to both children. The nodes of the second kind, the canonical pieces of $[c,d)$, are disjoint and make up $[c,d)$, because every leaf lies inside $[c,d)$ or misses it. A node goes on to its children only if $l\lt c\lt r$ or $l\lt d\lt r$, which holds for at most two nodes of each depth; so at most four nodes per depth are reached, and since each split halves the length, the tree has only $\lceil\log_2 z_{\max}\rceil=14$ levels below the root. ^[7] Every node reached that meets $[c,d)$ then recomputes its length $\lambda$, after its children: $\lambda=r-l$ if $\kappa\gt0$, $\lambda=0$ for a leaf with $\kappa=0$, and otherwise the sum of the lengths of the children. The covered length is the $\lambda$ of the root.

7. Why the root holds the covered length

For a node $t$, let $W_t$ be the union of the canonical pieces of current intervals that lie in the subtree of $t$, $t$ included. Claim: $\lambda_t$ always equals the length $|W_t|$; at the root, $W$ is the union of the current intervals. Removing $[c,d)$ reaches the same nodes as adding it, since the recursion depends only on $[c,d)$, and a rectangle is removed at $b_j\gt a_j$, after it was added; so $\kappa_t$ is always the number of current intervals of which $t$ is a canonical piece, never negative. Now go up from the leaves. If $\kappa_t\gt0$, then $t$ itself lies in $W_t$, so $|W_t|=r-l$. If $\kappa_t=0$, then $W_t$ is the union of the sets of the children, which lie in disjoint halves, so the lengths add, and a leaf has $W_t=\varnothing$. These are the three cases of the rule for $\lambda$. An update recomputes $\lambda$, children first, at every node it reaches that meets $[c,d)$; every other node misses $[c,d)$ or lies below a canonical piece, and nothing in its subtree has changed. No counter is ever pushed down to the children, and at the end of a slab, when every interval has been removed, all counters and lengths are $0$ again.

8. The example in the tree

In the example $z_{\max}=4$: the root $[0,4)$ has the children $[0,2)$ and $[2,4)$, and these have the unit intervals as children. The canonical pieces of $[0,3)$ are $[0,2)$ and $[2,3)$, those of $[1,4)$ are $[1,2)$ and $[2,4)$. In slab $1$:

$y=0$: adding $[0,3)$ sets $\kappa=1$ at $[0,2)$ and $[2,3)$, so $\lambda$ is $2$ at $[0,2)$, $1$ at $[2,4)$ and $3$ at the root;
$y=1$: adding $[1,4)$ sets $\kappa=1$ at $[1,2)$ and $[2,4)$; now $\lambda=2$ at $[2,4)$ and $\lambda=4$ at the root;
$y=2$: removing $[0,3)$ sets $\kappa=0$ at $[0,2)$ and $[2,3)$; the node $[0,2)$ now adds the lengths of its children, $0+1=1$, because the leaf $[1,2)$ still counts $[1,4)$, and the root has $\lambda=1+2=3$;
$y=3$: removing $[1,4)$ brings every counter and length back to $0$.

The root lengths $3$, $4$ and $3$ are those of Section 5. The third step shows why counters stay at the canonical pieces: the leaf $[1,2)$ remembers $[1,4)$ although the counter of its parent has dropped to $0$.

How the Code Works

The C++ program ^[8] runs validation_checks, generates the $50000$ cuboids and prints their combined volume. build_lfg_sequence fills a vector s with $S_1,\dots,S_{300000}$, evaluating the cubic in long long; its correction of negative remainders never fires, by Section 1. generate_cuboids stores each cuboid as a Cuboid with the half-open ends $x_0,x_1,y_0,y_1,z_0,z_1$ of Section 2. combined_volume finds $x_{\max}$ and $z_{\max}$ and counts the rectangles of every slab with a difference array: $+1$ at diff[x0], $-1$ at diff[x1], then running sums. ^[9] It reserves exactly that much space per slab and copies the rectangle of each cuboid, a RectYZ of four 16-bit numbers, into its $dx$ slab lists. The slab areas are added in unsigned 64-bit integers, since the volume is bounded only by the $1.12\cdot10^{12}$ cells of Section 2.

AreaUnion::compute handles one slab. It makes the two events of Section 5 for every rectangle, an Event with $y$, $z_0$, $z_1$ and $\pm1$, sorts them by $y$ with std::sort, resets the tree and walks through the list. At each new $y$ it first adds length_[1], the root length, times the distance from the previous $y$, and only then applies all events at this $y$, so every strip is measured with its own covered length. update(node, l, r, ql, qr, delta) is the recursion of Section 6 on the arrays cover_ and length_ of $4z_{\max}+8$ entries, in which node $v$ has the children $2v$ and $2v+1$ and the root is node $1$; its last lines are the rule for $\lambda$. By Section 7 the reset before each slab is not even needed.

validation_checks asserts the parameters of $C_1$ and $C_2$ and the checkpoint $723581599$, and compares combined_volume with brute_force_volume, which marks the unit cells of every cuboid in an array and counts them, on $8$ random collections of $30$ cuboids with lower corners from $0$ to $11$ and sides from $1$ to $4$, drawn from a linear congruential generator. Being assert statements, these comparisons are skipped when NDEBUG is defined.

The Python and Java versions compute the same sums without checks. In Python, solve appends one shared tuple $(y_0,y_1,z_0,z_1)$ per cuboid to its slab lists, and area_worker computes the areas of a block of consecutive slabs, with a nested update and fresh lists cover and length for every slab; solve adds the block sums, in unbounded integers. Java does everything in main, with six int arrays for the cuboids, four exactly sized arrays per slab, one reused event array int[2 * N][4] sorted by $y$ with Arrays.sort, and a static update on the range $[0,z_{\max}+1)$, one leaf more than needed; areas and the total are long.

Complexity and Verification

Generation takes $O(n)$ time for $n=50000$ cuboids, and setting up and filling the $x_{\max}$ slab lists takes $O(x_{\max}+M)$ time and memory, with $M=\sum dx\le399n$. A slab with $R$ rectangles costs $O(R\log R)$ for sorting and $O(R\log z_{\max})$ for its $2R$ updates, plus $O(z_{\max})$ for the reset, so the total time is $O(M\log n+M\log z_{\max}+x_{\max}z_{\max})$ and the memory $O(M+x_{\max}+z_{\max})$. For this input the sweeps process $20060486$ events. In an instrumented copy of the C++ sweep, whose tree Python shares, they cause $805108530$ calls of update, about $40$ per event and at most $55$; Java's tree, one leaf larger, needs $805448794$. The C++ resets write $863533664$ array entries in fast block fills. The slab lists dominate the memory, with $8$ bytes per rectangle in C++, about $80$ MB, and $16$ in Java, about $160$ MB. The C++ and Java programs finish within a few seconds on a desktop computer; Python needs more than a minute of computing time.

The C++ checks pass. For this article, outside the published solutions, we also verified:

all three programs print the same volume;
a separate program without sorting or segment tree, which sweeps along $x$ with a counter for every unit square of the $(y,z)$ range, adds the rectangle of a cuboid at $x_0$, removes it at $x_1$ and sums the numbers of nonzero counters, gives $723581599$ for $100$ cuboids and the same volume as the three programs for all $50000$;
a Python version of Sections 4 to 7 agrees with the cell count and with full inclusion–exclusion on $3000$ random collections of up to $9$ small cuboids, and with the separate program for the first $100$, $300$ and $1000$ cuboids;
for the first $100$ cuboids a grid cut at the $197$ distinct face coordinates of each axis gives $723581599$; testing all pairs and extending them to larger sets confirms the overlap counts of Section 3.

Footnotes and References

The references supply the background facts; the slab reduction, the tree invariant and all counts for the generated cuboids are derived or computed in this article.

Project Euler 212 — Combined Volume of Cuboids. The official statement defines closed axis-aligned cuboids by the parameters {(x₀, y₀, z₀), (dx, dy, dz)}, generates C₁, …, C₅₀₀₀₀ with the lagged Fibonacci generator, gives C₁, C₂ and the combined volume 723581599 of the first 100 cuboids, and asks for that of all 50000.
Wikipedia — Lagged Fibonacci generator. The generator forms Sₙ ≡ Sₙ₋ⱼ ⋆ Sₙ₋ₖ (mod m) with 0 < j < k; with addition it is called additive, and the lags (24, 55) are among the commonly used pairs, also in Knuth's list.
Wikipedia — Inclusion–exclusion principle. The size of a union is the alternating sum of the sizes of the intersections over all nonempty index sets; the identity also holds for measurable sets of finite measure.
Wikipedia — Klee's measure problem. It asks for the measure of a union of axis-aligned boxes. In two dimensions Bentley (1977) solved it in optimal O(n log n) time by sweeping a line; for d ≥ 3 the best known bound is Chan's O(n^(d/2)) from 2013.
Wikipedia — Cavalieri's principle. Solids whose cross-sections at every height have equal areas have equal volumes; the principle is a particular case of Fubini's theorem.
Wikipedia — Sweep line algorithm. A conceptual line moves across the plane, stopping at certain points to handle the objects that meet it.
Wikipedia — Segment tree. Invented by Jon Bentley in 1977, it stores an interval at no more than two nodes of each depth, and an insertion visits at most four nodes per level, O(log n) in all.
C++ — Euler212.cpp. The linked immutable C++ revision, with combined_volume, AreaUnion and validation_checks, is the implementation record for this article.
Wikipedia — Prefix sum. Prefix sums and finite differences are inverse operations, so running sums of the differences +1 at x₀ and −1 at x₁ count the active cuboids of each slab.

Problem 212 source code

C++

#include <algorithm>
#include <atomic>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <thread>
#include <vector>

namespace {

constexpr int LFG_MOD = 1000000;
constexpr int POS_MOD = 10000;
constexpr int LEN_MOD = 399;

struct Cuboid {
    int x0;
    int x1;
    int y0;
    int y1;
    int z0;
    int z1;
};

struct RectYZ {
    std::uint16_t y0;
    std::uint16_t y1;
    std::uint16_t z0;
    std::uint16_t z1;
};

struct Event {
    int y;
    int z0;
    int z1;
    int delta;

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

class AreaUnion {
public:
    explicit AreaUnion(int max_z)
        : max_z_(max_z), cover_(static_cast<std::size_t>(4 * max_z + 8), 0),
          length_(static_cast<std::size_t>(4 * max_z + 8), 0) {}

    std::uint64_t compute(const std::vector<RectYZ>& rects) {
        if (rects.empty()) return 0;

        events_.clear();
        events_.reserve(rects.size() * 2);
        for (const RectYZ& r : rects) {
            events_.push_back({static_cast<int>(r.y0), static_cast<int>(r.z0), static_cast<int>(r.z1), +1});
            events_.push_back({static_cast<int>(r.y1), static_cast<int>(r.z0), static_cast<int>(r.z1), -1});
        }
        std::sort(events_.begin(), events_.end());

        std::fill(cover_.begin(), cover_.end(), 0);
        std::fill(length_.begin(), length_.end(), 0);

        std::uint64_t area = 0;
        int prev_y = events_.front().y;
        std::size_t i = 0;

        while (i < events_.size()) {
            const int y = events_[i].y;
            area += static_cast<std::uint64_t>(length_[1]) * static_cast<std::uint64_t>(y - prev_y);

            while (i < events_.size() && events_[i].y == y) {
                update(1, 0, max_z_, events_[i].z0, events_[i].z1, events_[i].delta);
                ++i;
            }
            prev_y = y;
        }

        return area;
    }

private:
    void update(int node, int l, int r, int ql, int qr, int delta) {
        if (qr <= l || r <= ql) return;

        if (ql <= l && r <= qr) {
            cover_[static_cast<std::size_t>(node)] += delta;
        } else {
            const int mid = l + (r - l) / 2;
            update(node * 2, l, mid, ql, qr, delta);
            update(node * 2 + 1, mid, r, ql, qr, delta);
        }

        if (cover_[static_cast<std::size_t>(node)] > 0) {
            length_[static_cast<std::size_t>(node)] = r - l;
        } else if (r - l == 1) {
            length_[static_cast<std::size_t>(node)] = 0;
        } else {
            length_[static_cast<std::size_t>(node)] =
                length_[static_cast<std::size_t>(node * 2)] + length_[static_cast<std::size_t>(node * 2 + 1)];
        }
    }

    int max_z_;
    std::vector<int> cover_;
    std::vector<int> length_;
    std::vector<Event> events_;
};

std::vector<int> build_lfg_sequence(int needed) {
    std::vector<int> s(static_cast<std::size_t>(needed + 1), 0);

    for (int k = 1; k <= std::min(55, needed); ++k) {
        long long value = 100003LL - 200003LL * k + 300007LL * k * k * k;
        value %= LFG_MOD;
        if (value < 0) value += LFG_MOD;
        s[static_cast<std::size_t>(k)] = static_cast<int>(value);
    }

    for (int k = 56; k <= needed; ++k) {
        s[static_cast<std::size_t>(k)] =
            (s[static_cast<std::size_t>(k - 24)] + s[static_cast<std::size_t>(k - 55)]) % LFG_MOD;
    }

    return s;
}

std::vector<Cuboid> generate_cuboids(int n) {
    const int needed = 6 * n;
    std::vector<int> s = build_lfg_sequence(needed);

    std::vector<Cuboid> cuboids;
    cuboids.reserve(static_cast<std::size_t>(n));

    for (int i = 1; i <= n; ++i) {
        const int base = 6 * i;
        const int x0 = s[static_cast<std::size_t>(base - 5)] % POS_MOD;
        const int y0 = s[static_cast<std::size_t>(base - 4)] % POS_MOD;
        const int z0 = s[static_cast<std::size_t>(base - 3)] % POS_MOD;

        const int dx = 1 + (s[static_cast<std::size_t>(base - 2)] % LEN_MOD);
        const int dy = 1 + (s[static_cast<std::size_t>(base - 1)] % LEN_MOD);
        const int dz = 1 + (s[static_cast<std::size_t>(base)] % LEN_MOD);

        cuboids.push_back({x0, x0 + dx, y0, y0 + dy, z0, z0 + dz});
    }

    return cuboids;
}

std::uint64_t combined_volume(const std::vector<Cuboid>& cuboids, unsigned thread_count) {
    if (cuboids.empty()) return 0;

    int max_x = 0;
    int max_z = 0;
    for (const Cuboid& c : cuboids) {
        max_x = std::max(max_x, c.x1);
        max_z = std::max(max_z, c.z1);
    }

    std::vector<int> diff(static_cast<std::size_t>(max_x + 1), 0);
    for (const Cuboid& c : cuboids) {
        ++diff[static_cast<std::size_t>(c.x0)];
        --diff[static_cast<std::size_t>(c.x1)];
    }

    std::vector<int> counts(static_cast<std::size_t>(max_x), 0);
    int active = 0;
    for (int x = 0; x < max_x; ++x) {
        active += diff[static_cast<std::size_t>(x)];
        counts[static_cast<std::size_t>(x)] = active;
    }

    std::vector<std::vector<RectYZ>> slabs(static_cast<std::size_t>(max_x));
    for (int x = 0; x < max_x; ++x) {
        slabs[static_cast<std::size_t>(x)].reserve(static_cast<std::size_t>(counts[static_cast<std::size_t>(x)]));
    }

    for (const Cuboid& c : cuboids) {
        const RectYZ r{static_cast<std::uint16_t>(c.y0), static_cast<std::uint16_t>(c.y1),
                       static_cast<std::uint16_t>(c.z0), static_cast<std::uint16_t>(c.z1)};
        for (int x = c.x0; x < c.x1; ++x) {
            slabs[static_cast<std::size_t>(x)].push_back(r);
        }
    }

    if (thread_count == 0) thread_count = 1;
    thread_count = std::min<unsigned>(thread_count, static_cast<unsigned>(max_x));

    std::atomic<int> next_x{0};
    std::vector<std::uint64_t> partial(static_cast<std::size_t>(thread_count), 0);
    std::vector<std::thread> workers;
    workers.reserve(static_cast<std::size_t>(thread_count));

    for (unsigned tid = 0; tid < thread_count; ++tid) {
        workers.emplace_back([&, tid]() {
            AreaUnion area_union(max_z);
            std::uint64_t local = 0;

            while (true) {
                const int x = next_x.fetch_add(1, std::memory_order_relaxed);
                if (x >= max_x) break;
                local += area_union.compute(slabs[static_cast<std::size_t>(x)]);
            }

            partial[static_cast<std::size_t>(tid)] = local;
        });
    }

    for (auto& th : workers) th.join();

    std::uint64_t total = 0;
    for (std::uint64_t v : partial) total += v;
    return total;
}

std::uint64_t brute_force_volume(const std::vector<Cuboid>& cuboids) {
    int max_x = 0;
    int max_y = 0;
    int max_z = 0;
    for (const Cuboid& c : cuboids) {
        max_x = std::max(max_x, c.x1);
        max_y = std::max(max_y, c.y1);
        max_z = std::max(max_z, c.z1);
    }

    const int xy = max_x * max_y;
    std::vector<unsigned char> occ(static_cast<std::size_t>(xy * max_z), 0);

    for (const Cuboid& c : cuboids) {
        for (int x = c.x0; x < c.x1; ++x) {
            for (int y = c.y0; y < c.y1; ++y) {
                const int base = (x * max_y + y) * max_z;
                for (int z = c.z0; z < c.z1; ++z) {
                    occ[static_cast<std::size_t>(base + z)] = 1;
                }
            }
        }
    }

    std::uint64_t total = 0;
    for (unsigned char b : occ) total += b;
    return total;
}

void validation_checks() {
    {
        const std::vector<Cuboid> c = generate_cuboids(2);
        assert(c[0].x0 == 7 && c[0].y0 == 53 && c[0].z0 == 183);
        assert(c[0].x1 - c[0].x0 == 94 && c[0].y1 - c[0].y0 == 369 && c[0].z1 - c[0].z0 == 56);
        assert(c[1].x0 == 2383 && c[1].y0 == 3563 && c[1].z0 == 5079);
        assert(c[1].x1 - c[1].x0 == 42 && c[1].y1 - c[1].y0 == 212 && c[1].z1 - c[1].z0 == 344);
    }

    {
        const std::vector<Cuboid> c100 = generate_cuboids(100);
        assert(combined_volume(c100, 1) == 723581599ULL);
    }

    {
        std::uint32_t seed = 1;
        auto next_rand = [&seed]() {
            seed = seed * 1664525u + 1013904223u;
            return seed;
        };

        for (int tc = 0; tc < 8; ++tc) {
            std::vector<Cuboid> small;
            for (int i = 0; i < 30; ++i) {
                const int x0 = static_cast<int>(next_rand() % 12);
                const int y0 = static_cast<int>(next_rand() % 12);
                const int z0 = static_cast<int>(next_rand() % 12);
                const int dx = 1 + static_cast<int>(next_rand() % 4);
                const int dy = 1 + static_cast<int>(next_rand() % 4);
                const int dz = 1 + static_cast<int>(next_rand() % 4);
                small.push_back({x0, x0 + dx, y0, y0 + dy, z0, z0 + dz});
            }

            const std::uint64_t fast = combined_volume(small, 1);
            const std::uint64_t brute = brute_force_volume(small);
            assert(fast == brute);
        }
    }
}

}  // namespace

int main() {
    validation_checks();

    const std::vector<Cuboid> cuboids = generate_cuboids(50000);

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) threads = 1;

    const std::uint64_t answer = combined_volume(cuboids, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

from __future__ import annotations

import multiprocessing as mp
import os
from concurrent.futures import ProcessPoolExecutor


LFG_MOD = 1_000_000
POS_MOD = 10_000
LEN_MOD = 399

G_SLABS: list[list[tuple[int, int, int, int]]] = []
G_MAX_Z = 0


def build_lfg_sequence(needed: int) -> list[int]:
    s = [0] * (needed + 1)
    for k in range(1, min(56, needed + 1)):
        s[k] = (100003 - 200003 * k + 300007 * k * k * k) % LFG_MOD
    for k in range(56, needed + 1):
        s[k] = (s[k - 24] + s[k - 55]) % LFG_MOD
    return s


def generate_cuboids(n: int) -> list[tuple[int, int, int, int, int, int]]:
    needed = 6 * n
    s = build_lfg_sequence(needed)
    cuboids: list[tuple[int, int, int, int, int, int]] = []
    cuboids_append = cuboids.append
    for i in range(1, n + 1):
        base = 6 * i
        x0 = s[base - 5] % POS_MOD
        y0 = s[base - 4] % POS_MOD
        z0 = s[base - 3] % POS_MOD
        dx = 1 + (s[base - 2] % LEN_MOD)
        dy = 1 + (s[base - 1] % LEN_MOD)
        dz = 1 + (s[base] % LEN_MOD)
        cuboids_append((x0, x0 + dx, y0, y0 + dy, z0, z0 + dz))
    return cuboids


def init_area_worker(slabs: list[list[tuple[int, int, int, int]]], max_z: int) -> None:
    global G_SLABS, G_MAX_Z
    G_SLABS = slabs
    G_MAX_Z = max_z


def area_worker(bounds: tuple[int, int]) -> int:
    start_x, end_x = bounds
    slabs = G_SLABS
    max_z = G_MAX_Z
    tree_size = 4 * max_z + 8

    def update(node: int, l: int, r: int, ql: int, qr: int, delta: int, cover: list[int], length: list[int]) -> None:
        if qr <= l or r <= ql:
            return
        if ql <= l and r <= qr:
            cover[node] += delta
        else:
            mid = l + ((r - l) >> 1)
            update(node * 2, l, mid, ql, qr, delta, cover, length)
            update(node * 2 + 1, mid, r, ql, qr, delta, cover, length)

        if cover[node] > 0:
            length[node] = r - l
        elif r - l == 1:
            length[node] = 0
        else:
            length[node] = length[node * 2] + length[node * 2 + 1]

    total = 0
    for x in range(start_x, end_x):
        rects = slabs[x]
        if not rects:
            continue

        events: list[tuple[int, int, int, int]] = []
        events_extend = events.extend
        for y0, y1, z0, z1 in rects:
            events_extend(((y0, z0, z1, 1), (y1, z0, z1, -1)))
        events.sort(key=lambda e: e[0])

        cover = [0] * tree_size
        length = [0] * tree_size
        area = 0
        prev_y = events[0][0]
        i = 0
        m = len(events)
        while i < m:
            y = events[i][0]
            area += length[1] * (y - prev_y)
            while i < m and events[i][0] == y:
                _, z0, z1, d = events[i]
                update(1, 0, max_z, z0, z1, d, cover, length)
                i += 1
            prev_y = y
        total += area

    return total


def solve() -> int:
    n = 50_000
    cuboids = generate_cuboids(n)

    max_x = max(c[1] for c in cuboids)
    max_z = max(c[5] for c in cuboids)

    slabs: list[list[tuple[int, int, int, int]]] = [[] for _ in range(max_x)]
    for x0, x1, y0, y1, z0, z1 in cuboids:
        r = (y0, y1, z0, z1)
        for x in range(x0, x1):
            slabs[x].append(r)

    threads = max(1, min(8, os.cpu_count() or 1))
    if threads == 1:
        init_area_worker(slabs, max_z)
        return area_worker((0, max_x))

    block = (max_x + threads - 1) // threads
    jobs = []
    for t in range(threads):
        start_x = t * block
        end_x = min(max_x, (t + 1) * block)
        if start_x >= end_x:
            break
        jobs.append((start_x, end_x))

    ctx = mp.get_context("fork") if "fork" in mp.get_all_start_methods() else mp.get_context()
    total = 0
    with ProcessPoolExecutor(
        max_workers=len(jobs),
        mp_context=ctx,
        initializer=init_area_worker,
        initargs=(slabs, max_z),
    ) as pool:
        for part in pool.map(area_worker, jobs):
            total += part
    return total


if __name__ == "__main__":
    print(solve())

Java

import java.util.*;

public class Euler212 {
    static int[] cover;
    static int[] len;

    public static void main(String[] args) {
        int LFG_MOD = 1000000, POS_MOD = 10000, LEN_MOD = 399;
        int N = 50000, needed = 6 * N;
        int[] s = new int[needed + 1];
        for (int k = 1; k <= Math.min(55, needed); k++) {
            long v = 100003L - 200003L * k + 300007L * k * k * k;
            v %= LFG_MOD;
            if (v < 0)
                v += LFG_MOD;
            s[k] = (int) v;
        }
        for (int k = 56; k <= needed; k++)
            s[k] = (s[k - 24] + s[k - 55]) % LFG_MOD;

        int[] x0a = new int[N], x1a = new int[N];
        int[] y0a = new int[N], y1a = new int[N];
        int[] z0a = new int[N], z1a = new int[N];
        int maxX = 0, maxZ = 0;
        for (int i = 0; i < N; i++) {
            int b = 6 * (i + 1);
            x0a[i] = s[b - 5] % POS_MOD;
            y0a[i] = s[b - 4] % POS_MOD;
            z0a[i] = s[b - 3] % POS_MOD;
            int dx = 1 + s[b - 2] % LEN_MOD, dy = 1 + s[b - 1] % LEN_MOD, dz = 1 + s[b] % LEN_MOD;
            x1a[i] = x0a[i] + dx;
            y1a[i] = y0a[i] + dy;
            z1a[i] = z0a[i] + dz;
            maxX = Math.max(maxX, x1a[i]);
            maxZ = Math.max(maxZ, z1a[i]);
        }
        int sz = maxZ + 1;
        cover = new int[4 * sz];
        len = new int[4 * sz];

        // Diff array for cuboid count at each x
        int[] diff = new int[maxX + 1];
        for (int i = 0; i < N; i++) {
            diff[x0a[i]]++;
            diff[x1a[i]]--;
        }
        int[] counts = new int[maxX];
        int act = 0;
        for (int x = 0; x < maxX; x++) {
            act += diff[x];
            counts[x] = act;
        }

        // Build slab lists per x
        int[][] slabY0 = new int[maxX][], slabY1 = new int[maxX][];
        int[][] slabZ0 = new int[maxX][], slabZ1 = new int[maxX][];
        int[] fill = new int[maxX];
        for (int x = 0; x < maxX; x++) {
            slabY0[x] = new int[counts[x]];
            slabY1[x] = new int[counts[x]];
            slabZ0[x] = new int[counts[x]];
            slabZ1[x] = new int[counts[x]];
        }
        for (int i = 0; i < N; i++)
            for (int x = x0a[i]; x < x1a[i]; x++) {
                int f = fill[x];
                slabY0[x][f] = y0a[i];
                slabY1[x][f] = y1a[i];
                slabZ0[x][f] = z0a[i];
                slabZ1[x][f] = z1a[i];
                fill[x]++;
            }

        long total = 0;
        int[][] events = new int[2 * N][4];
        for (int x = 0; x < maxX; x++) {
            int nR = fill[x];
            if (nR == 0)
                continue;
            int ne = 0;
            for (int i = 0; i < nR; i++) {
                events[ne][0] = slabY0[x][i];
                events[ne][1] = slabZ0[x][i];
                events[ne][2] = slabZ1[x][i];
                events[ne][3] = 1;
                ne++;
                events[ne][0] = slabY1[x][i];
                events[ne][1] = slabZ0[x][i];
                events[ne][2] = slabZ1[x][i];
                events[ne][3] = -1;
                ne++;
            }
            Arrays.sort(events, 0, ne, (a, b) -> a[0] - b[0]);
            Arrays.fill(cover, 0);
            Arrays.fill(len, 0);
            long area = 0;
            int prevY = events[0][0], ei = 0;
            while (ei < ne) {
                int y = events[ei][0];
                area += (long) len[1] * (y - prevY);
                while (ei < ne && events[ei][0] == y) {
                    update(1, 0, sz, events[ei][1], events[ei][2], events[ei][3]);
                    ei++;
                }
                prevY = y;
            }
            total += area;
        }
        System.out.println(total);
    }

    static void update(int node, int l, int r, int ql, int qr, int delta) {
        if (qr <= l || r <= ql)
            return;
        if (ql <= l && r <= qr)
            cover[node] += delta;
        else {
            int mid = l + (r - l) / 2;
            update(2 * node, l, mid, ql, qr, delta);
            update(2 * node + 1, mid, r, ql, qr, delta);
        }
        if (cover[node] > 0)
            len[node] = r - l;
        else if (r - l == 1)
            len[node] = 0;
        else
            len[node] = len[2 * node] + len[2 * node + 1];
    }
}