Solution

Answer: 328968937309

Let the cuboids be interpreted as half-open boxes

$$[x_0,x_0+dx)\times [y_0,y_0+dy)\times [z_0,z_0+dz)$$

which has the same volume as the closed version in the statement, but avoids boundary-overlap ambiguity.

The task is to compute the volume of the union of $50{,}000$ axis-aligned cuboids.

Mathematical analysis

A direct inclusion–exclusion approach is impossible: there are far too many cuboids.

The standard way to compute union volume of axis-aligned boxes is a sweep-line algorithm.

1. Sweep in the $x$-direction

Each cuboid contributes:

an entering event at $x=x_0$
a leaving event at $x=x_0+dx$

Between two consecutive $x$-coordinates, the active set of cuboids does not change.

So if the union area of the active rectangles in the $yz$-plane is $A(x)$, then

$$V = \sum (\Delta x)\cdot A(x)$$

In this problem all coordinates are integers and

$$0 \le x < 10400$$

because:

$x_0 < 10000$
$dx \le 399$

Hence we only need to process $10400$ integer slabs.

2. Computing union area in the $yz$-plane

For the currently active cuboids, project them onto the $yz$-plane.

Each becomes a rectangle

$$[y_1,y_2)\times[z_1,z_2)$$

We again use a sweep-line:

sweep in $y$
maintain covered length in $z$

For each rectangle:

add event $(y_1, z_1,z_2,+1)$
add event $(y_2, z_1,z_2,-1)$

A segment tree over the $z$-axis stores:

coverage count
total covered length

If a node’s count is positive, its whole interval is covered.

Thus the union area is accumulated as:

$$\text{area} += (y_i-y_{i-1}) \times (\text{covered z-length})$$

Python implementation

from collections import defaultdict

# ------------------------------------------------------------
# Lagged Fibonacci Generator
# ------------------------------------------------------------

S = [0] * 300001

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

for k in range(56, 300001):
    S[k] = (S[k - 24] + S[k - 55]) % 1000000

# ------------------------------------------------------------
# Build cuboids and x-events
# ------------------------------------------------------------

ADD = defaultdict(list)
REMOVE = defaultdict(list)

for n in range(1, 50001):

    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)

    cuboid = (y0, y0 + dy, z0, z0 + dz)

    ADD[x0].append(cuboid)
    REMOVE[x0 + dx].append(cuboid)

# ------------------------------------------------------------
# Segment tree for union length on z-axis
# ------------------------------------------------------------

ZMAX = 10400

class SegmentTree:

    def __init__(self, n):
        self.n = n
        self.count = [0] * (4 * n)
        self.length = [0] * (4 * n)

    def update(self, node, left, right, ql, qr, val):

        if qr <= left or right <= ql:
            return

        if ql <= left and right <= qr:
            self.count[node] += val
        else:
            mid = (left + right) // 2

            self.update(node * 2, left, mid, ql, qr, val)
            self.update(node * 2 + 1, mid, right, ql, qr, val)

        if self.count[node] > 0:
            self.length[node] = right - left
        else:
            if right - left == 1:
                self.length[node] = 0
            else:
                self.length[node] = (
                    self.length[node * 2]
                    + self.length[node * 2 + 1]
                )

# ------------------------------------------------------------
# Main sweep in x
# ------------------------------------------------------------

active = []
volume = 0

for x in range(10400):

    # add cuboids starting here
    active.extend(ADD[x])

    # remove cuboids ending here
    for r in REMOVE[x]:
        active.remove(r)

    # build y-events
    events = []

    for y1, y2, z1, z2 in active:
        events.append((y1, z1, z2, +1))
        events.append((y2, z1, z2, -1))

    events.sort()

    # sweep in y
    seg = SegmentTree(ZMAX)

    area = 0
    prev_y = 0
    i = 0

    while i < len(events):

        y = events[i][0]

        area += (y - prev_y) * seg.length[1]

        while i < len(events) and events[i][0] == y:
            _, z1, z2, typ = events[i]
            seg.update(1, 0, ZMAX, z1, z2, typ)
            i += 1

        prev_y = y

    volume += area

print(volume)

Code walkthrough

Generating the cuboids

The Lagged Fibonacci Generator produces $300000$ numbers.

Every six numbers define one cuboid:

$$(x_0,y_0,z_0,dx,dy,dz)$$

matching the statement exactly.

Event construction

For each cuboid:

add it at $x_0$
remove it at $x_0+dx$

This lets us maintain the active set while sweeping in $x$.

Segment tree

The segment tree maintains the union length on the $z$-axis.

Each node stores:

count[node] = how many intervals fully cover that node
length[node] = total covered length

If coverage count is positive, the interval is entirely covered.

Otherwise, coverage comes from children.

Inner sweep in $y$

For active rectangles:

entering event at $y_1$
leaving event at $y_2$

Between consecutive $y$-values, covered $z$-length is constant.

So:

$$\Delta \text{area} = (\Delta y)\times(\text{covered z length})$$

Verification

The problem states:

The combined volume of the first 100 cuboids is 723581599.

Running the same algorithm on only the first 100 cuboids reproduces exactly:

$$723581599$$

so the implementation is validated.

The final volume is about $3.29\times 10^{11}$, which is reasonable:

the total bounding box volume is roughly $10^4{}^3 = 10^{12}$
there are many overlaps
each cuboid is relatively small

So the magnitude is consistent.

Answer: 328968937309

Project Euler Problem 212