Advent of Code 2025 - Day 8 Solution

Core Problem Description

Given $N$ points in three-dimensional space, where each point’s coordinates $(x, y, z)$ are non-negative integers. It is also ensured that each point’s coordinates are unique, and the Euclidean distances between any two points are all distinct. Compute:

1. The first $K$ shortest distances among all point pairs, then output the product of the sizes of the top three connected components after connecting these distances as edges;
2. Following the order of edge addition, output the product of the $x$ coordinates of the two points connected by the last edge added before all points become connected.

Thought Process

Problem Breakdown

The problem’s requirements are complex, so we first need to break it down into several sub-problems. For the first part, we need to:

1. Find the first $K$ shortest Euclidean distances among all point pairs;
2. Use these distances as edges to build a graph and find connected components;
3. Compute the product of the sizes of the top three connected components.

For the second part, we need to:

1. Connect the points step by step in the order of edge addition;
2. Before all points are connected, record the product of the $x$ coordinates of the two points connected by the last added edge.

It’s not difficult to see that for the first part we can use a Disjoint Set Union (DSU) to manage connected components, and for the second part we similarly need to use DSU to check the connection status of points. This is not particularly challenging, so the main difficulty lies in the first step of each problem: how to efficiently find the first $K$ shortest distances among all point pairs, or to connect the points in distance order.

Parallel Optimization

For finding the first $K$ shortest distances and the top $3$ largest connected components, we can use max/min heaps to assist the computation. Of course, the most time-consuming step is calculating the distances between all point pairs. To address this, we can employ some simple parallel optimization techniques:

1. Replace Euclidean Distance: Since we only need to compare the relative sizes of distances and not the exact values, we can use the square of the Euclidean distance instead of the Euclidean distance itself for comparison. This way we only need integer multiplication and addition, avoiding the square root operation for floating-point numbers, thereby improving calculation speed.
2. Avoid Redundant Calculations: Since Euclidean distance is symmetric, we only need to calculate each pair once. We can only calculate pairs $(p_i, p_j)$ where $i < j$.
3. Parallel Processing: Using a parallel computing library, we can split all calculation tasks into $M$ parts, each processed by a thread, and finally merge the results. This can significantly reduce computation time, especially on multi-core processors.

Parallel Reduction

In parallel computing, reducing is a common computational pattern: for $M$ input data of type $A$, reduce them to a single output of type $A$ using a binary operator $f: A\times A\to A$. In parallelized reduction, $f$ must be associative, meaning for any $a, b, c$, $f(f(a, b), c) \equiv f(a, f(b, c))$ holds. For example, in this problem, when merging two heaps, the output is still a heap, and the operation is associative. For such operations, we can partition the input data into several pieces, perform the reduction operation in multiple threads, a common pattern being the tree structure shown below:

If it’s a folding operation, where the function type is $g: A\times B\to A$, then during parallelization, typically only the first layer’s operation is replaced with $g: A\times B\to A$, while the remaining layers still use $f: A\times A\to A$ for reduction. This is also the pattern for implementing parallel folds in rayon.

With this, we can already solve the first part:

1
let mut dsu = DisjointSet::new(self.nodes.nrows());
2
(0..self.nodes.nrows())
3
    .into_par_iter()
4
    .flat_map_iter(|i| {
5

6
        ((i + 1)..self.nodes.nrows())
7
            .map(|j| (i, j))
8
            .collect::<Vec<_>>()
9
    })
10

11
    .fold(BinaryHeap::new, |mut heap, (i, j)| {
12
        let d = self.dist(i, j);
13
        heap.push((d, i, j));
14
        if heap.len() > self.max_steps {
15
            heap.pop();
16
        }
17
        heap
18
    })
19

20
    .reduce(BinaryHeap::new, |mut acc, mut heap| {
21
        acc.extend(heap.drain());
22
        while acc.len() > self.max_steps {
23
            acc.pop();
24
        }
25
        acc
26
    })
27
    .into_iter()
28
    .for_each(|(_, i, j)| dsu.union(i, j));
29
dsu.sizes
30
    .values()
31

32
    .fold(BinaryHeap::new(), |mut heap, &size| {
33
        heap.push(Reverse(size));
34
        if heap.len() > 3 {
35
            heap.pop();
36
        }
37
        heap
38
    })
39
    .iter()
40
    .map(|&Reverse(x)| x)
41
    .product::<u64>()
42
    .to_string()

For the second part, directly calculating all edges and sorting them is clearly impractical, as the number of edges we need to handle is on the order of $O(N^2)$. Therefore, we need some optimization to reduce computation. It’s clear that since we only care about the last edge that connects all points into a single connected component, we don’t actually need to care about edges that do not connect different connected components. Based on this, we can adopt the following strategy:

1. Edge Set Based on Points’ Shortest Edge: For each point, we only calculate the shortest edge connecting it to a point belonging to a different connected component, and add these edges to a candidate edge set. This can greatly reduce the number of edges to compute.
2. Dynamically Update the Edge Set: Each time we take the shortest edge $(u,v)$ from the edge set, we need to recalculate the shortest edge from $u$ to other connected components and add this edge to the edge set. This ensures the edge set is updated only when necessary, avoiding unnecessary computations.
3. Loop Until All Points Are Connected: Repeat the above process until all points belong to the same connected component. During this process, record the last edge $(u,v)$ that connected different connected components and output the product of its $x$ coordinates.

Following this idea, we can implement the following code:

1
let mut closest_neighbor = (0..self.nodes.nrows())
2
    .into_par_iter()
3
    .map(|i| {
4
        (0..self.nodes.nrows())
5
            .filter_map(|j| (j != i).then_some((self.dist(i, j), i, j)))
6
            .min_by_key(|&(dist, _, _)| dist)
7
            .map_or_else(
8
                || unreachable!("There should be at least one other node"),
9
                Reverse,
10
            )
11
    })
12
    .collect::<BinaryHeap<_>>();
13
let mut dsu = DisjointSet::new(self.nodes.nrows());
14
loop {
15
    let Some(Reverse((_, i, j))) = closest_neighbor.pop() else {
16
        panic!("No more edges to process");
17
    };
18

19
    let root_i = dsu.find(i);
20
    let root_j = dsu.find(j);
21
    if root_i != root_j {
22
        dsu.union(i, j);
23
    }
24

25
    if dsu.sizes.len() == 1 {
26
        return (self.nodes[[i, 0]] * self.nodes[[j, 0]]).to_string();
27
    }
28

29
    closest_neighbor.push(
30
        (0..self.nodes.nrows())
31
            .filter_map(|k| (root_i != dsu.find(k)).then_some((self.dist(i, k), i, k)))
32
            .min_by_key(|&(dist, _, _)| dist)
33
            .map_or_else(
34
                || unreachable!("At least one different component should exist"),
35
                Reverse,
36
            ),
37
    );
38
}

Octree Representation

In the above methods, although we use parallel computing and some dynamic update strategies to reduce computation time, the computational complexity is still $O(N^2)$. Facing very large datasets, it may still be impossible to complete within a reasonable time. Even with only $N=1,000$ in this problem, the above algorithm requires $21 ms$ and $281 ms$. Therefore, we need to further optimize the algorithm.

For the problem of distances between point pairs in space, we naturally think of using an octree to partition the space. In an octree, space is recursively divided into eight quadrants (as shown in the figure below), allowing us to quickly locate the spatial region of a point and theoretically query a point’s neighboring points efficiently.

It’s worth noting that in this problem, since the coordinates of each point are integers, we can actually use Morton encoding, also known as the Z-order curve, to partition space. Morton encoding interleaves the binary bits of each coordinate to generate a one-dimensional encoded value, achieving a linear representation of space, as shown below:

Magic Bit Operations

As Jeroen Baert described in his blog, we can use some “magic bit operations” to efficiently interleave binary bits. Specifically, through a series of bit shifts and mask operations, we can scatter the binary bits of each coordinate (and vice versa), achieving interleaving. As shown below:

1
const fn interleave_bits_by_3(x: i64) -> u64 {
2
    let mut x = x.cast_unsigned() & 0x1ff_fff;
3
    x = (x | (x << 32)) & 0x1f_000_000_00f_fff;
4
    x = (x | (x << 16)) & 0x1f_000_0ff_000_0ff;
5
    x = (x | (x << 8)) & 0x1_00f_00f_00f_00f_00f;
6
    x = (x | (x << 4)) & 0x1_0c3_0c3_0c3_0c3_0c3;
7
    x = (x | (x << 2)) & 0x1_249_249_249_249_249;
8
    x
9
}

Using Morton encoding, we can sort the points in space according to their Morton encoded values, resulting in a fractal curve similar to the one below:

The Morton encoded value implicitly contains hierarchical partitioning of a multi-dimensional space. Since we interleave the binary bits of the three-dimensional coordinates, every three bits in the encoding can be seen as a level, and the three-bit binary number at each level indicates one of the eight quadrants at that level. Therefore, by sorting according to Morton encoding, we have implicitly built an octree structure; any octree node can be represented as a prefix of a Morton encoded value and/or a corresponding array interval.

Branch and Bound

With the octree structure, we can leverage spatial partitioning to accelerate nearest neighbor search, which is branch and bound. During the search, after obtaining several nearest neighbor candidate points, we can use the current shortest distance to prune nodes in the octree. Specifically, assuming the current nearest neighbor distance is $d$, when we encounter an octree node during the search, we can calculate the minimum distance $d_{\min}$ between the space covered by that node and the query point. If $d_{\min} \geq d$, it means that all points within that node and its child nodes cannot be closer neighbors, so we can skip the node directly, avoiding unnecessary calculations.

Since in both parts of the problem, we need to consider distances between all points, we can instead use the dual-tree traversal method to batch process calculations between point pairs. For example, for the first part, we can construct all points as the root node of the octree, then traverse:

1. For the current node pair $(L, R)$, calculate the minimum distance $d_{\min}$ between their covered spaces;
2. (Bound) If $d_{\min}$ is greater than the currently known $K$-th largest distance, skip this node pair;
3. (Batch Calculation) Otherwise, if the sizes of $L$ and $R$ are below a certain threshold, directly calculate the distances between all point pairs in $L$ and $R$, and update the nearest neighbor heap;
4. (Branch) Otherwise, generate new node pairs to continue traversal according to the following rules:
   • If $L$ and $R$ are the same node, generate all node pairs combining child nodes without repetition ($i \le j$);
   • Otherwise, split the larger node and generate all combinations of its child nodes with the other node.

When calculating $d_{\min}$, we can use the properties of Morton encoding to quickly estimate the minimum distance between the spaces covered by two nodes. Specifically, each node corresponds to a prefix of a Morton encoded value, and all child nodes’ corresponding Morton encoded values start with that prefix. Therefore, we can extract the corresponding coordinate range from this prefix, obtaining the node’s corresponding Axis-Aligned Bounding Box (AABB). Then, we can use AABB distance calculation methods to get $d_{\min}$ in $O(1)$ time.

AABB Estimated from Prefix vs. Actual AABB

The AABB estimated from the prefix is the theoretically maximum bounding box; the actual point set contained by the node may be distributed more compactly, resulting in a smaller actual AABB, as shown below:

When the point set distribution is relatively uniform, this difference is usually small; but when the point set distribution is extremely uneven, the difference can be significant, affecting pruning efficiency. Therefore, in actual implementation, we can pre-calculate the actual AABB for each node when building the octree, improving the hit rate of pruning paths.

In this way, we can efficiently calculate the first $K$ shortest distances among all point pairs, thus solving the first part. But what about the second part? Since in the second part we need to dynamically update connected components, we cannot know in advance how many edges need to be calculated like in the first part, so we cannot directly use the dual-tree traversal method. Therefore, further improvement is needed.

Minimum Spanning Tree

Actually, we can re-examine the requirements of the second part. In the previous analysis, we concluded that we only need to focus on edges connecting different connected components. This process is essentially the Kruskal algorithm for constructing a Minimum Spanning Tree (MST). And the last edge is the longest edge in this MST. Due to the uniqueness of the MST, if we can directly construct the MST, then we can find this edge faster.

Indeed, we have faster algorithms for constructing the minimum spanning tree, such as Borůvka’s algorithm. The core idea of this algorithm is:

1. Initially, each point is an independent connected component;
2. For each connected component, find its shortest edge to any other connected component;
3. Add these edges to the minimum spanning tree and merge these connected components;
4. Repeat steps 2 and 3 until all points belong to the same connected component.

By combining the dual-tree traversal method with the octree, we can efficiently find the shortest edge for each connected component, thus quickly constructing the minimum spanning tree. We can also add an extra pruning condition: if the two nodes currently being processed belong to the same connected component, skip that node pair. The complete algorithm flow is as follows:

1. Initialize all points as independent connected components;
2. Initialize the shortest distance for each connected component to infinity;
3. Push the root node pair $(root, root)$ onto the stack;
4. Repeat the following steps until the stack is empty:
   1. For the node pair $(L, R)$ at the top of the stack, calculate the minimum distance $d_{\min}$ between their covered spaces;
   2. If $d_{\min}$ is greater than the maximum of the shortest distances of the connected components contained in $L$ and $R$, skip this node pair;
   3. If $L$ and $R$ belong to the same connected component, skip this node pair;
   4. If the size of $L \times R$ is below a certain threshold, directly calculate the distances between all point pairs in $L$ and $R$, and update the shortest distance for each connected component;
   5. Otherwise, generate new node pairs to push onto the stack according to the following rules:
      • If $L$ and $R$ are the same node, generate all node pairs combining child nodes without repetition ($i \le j$);
      • Otherwise, split the larger node and generate all combinations of its child nodes with the other node.
5. Use the shortest edge of each connected component to merge connected components, and record the longest edge updated during this process.
6. Repeat steps 2 to 5 until all points belong to the same connected component.

Core Algorithm Implementation

In the actual implementation for this problem, due to the small scale of the point set, we can adopt some detailed optimizations:

1. Fully construct the octree, determining the point set interval for each node to avoid frequent point set partitioning during traversal;
2. Pre-calculate the actual bounding box for each node to obtain a tighter bounding box than estimated from the Morton encoding prefix, improving pruning path hit rate.
3. In the second part, pre-calculate the connected components contained in all nodes for each iteration to avoid frequent Disjoint Set Union queries during traversal.

1
impl Puzzle {
2

3
    fn deeper_search_nodes(
4
        &self,
5
        left_node_idx: usize,
6
        right_node_idx: usize,
7
    ) -> Vec<(usize, usize)> {
8

9
        if left_node_idx == right_node_idx {
10
            let children = &self.octree_children[left_node_idx];
11
            return children
12
                .iter()
13
                .enumerate()
14
                .flat_map(|(idx, &i)| children[idx..].iter().map(move |&j| (i, j)))
15
                .collect();
16
        }
17

18
        if self.octree[left_node_idx].size() > self.octree[right_node_idx].size() {
19
            self.octree_children[left_node_idx]
20
                .iter()
21
                .map(|&child| (child, right_node_idx))
22
                .collect()
23
        } else {
24
            self.octree_children[right_node_idx]
25
                .iter()
26
                .map(|&child| (left_node_idx, child))
27
                .collect()
28
        }
29
    }
30
}
31

32
impl Solution for Puzzle {
33
    fn part1(&self) -> String {
34
        let mut min_dist = BinaryHeap::new();
35
        let mut search_stack = vec![(0, 0)];
36
        while let Some((left_node_idx, right_node_idx)) = search_stack.pop() {
37
            let left_node = self.octree[left_node_idx];
38
            let right_node = self.octree[right_node_idx];
39
            if let Some(&(delta, _, _)) = min_dist.peek() {
40

41
                let (left_min, left_max) = self.octree_bounding_boxes[left_node_idx];
42
                let (right_min, right_max) = self.octree_bounding_boxes[right_node_idx];
43
                let mut min_possible_dist = 0;
44
                for dim in 0..3 {
45
                    if left_max[dim] < right_min[dim] {
46
                        min_possible_dist += (right_min[dim] - left_max[dim]).pow(2);
47
                    } else if right_max[dim] < left_min[dim] {
48
                        min_possible_dist += (left_min[dim] - right_max[dim]).pow(2);
49
                    }
50
                }
51
                if min_possible_dist >= delta {
52
                    continue;
53
                }
54
            }
55

56
            if left_node.size() * right_node.size() <= self.max_steps {
57
                let pairs = if left_node_idx == right_node_idx {
58

59
                    (left_node.start..left_node.end)
60
                        .flat_map(|i| (i + 1..left_node.end).map(move |j| (i, j)))
61
                        .collect::<Vec<_>>()
62
                } else {
63

64
                    (left_node.start..left_node.end)
65
                        .flat_map(|i| (right_node.start..right_node.end).map(move |j| (i, j)))
66
                        .collect::<Vec<_>>()
67
                };
68
                for (i, j) in pairs {
69
                    let d = Self::dist(&self.nodes[i], &self.nodes[j]);
70
                    if d > min_dist.peek().map_or(i64::MAX, |&(d, _, _)| d) {
71
                        continue;
72
                    }
73
                    min_dist.push((d, i, j));
74
                    if min_dist.len() > self.max_steps {
75
                        min_dist.pop();
76
                    }
77
                }
78
                continue;
79
            }
80
            search_stack.extend(self.deeper_search_nodes(left_node_idx, right_node_idx));
81
        }
82
        let mut dsu = DisjointSet::new(self.nodes.len());
83
        for (_, i, j) in min_dist {
84
            dsu.union(i, j);
85
        }
86
        dsu.sizes
87
            .values()
88
            .fold(BinaryHeap::new(), |mut heap, &size| {
89
                heap.push(Reverse(size));
90
                if heap.len() > 3 {
91
                    heap.pop();
92
                }
93
                heap
94
            })
95
            .iter()
96
            .map(|&Reverse(x)| x)
97
            .product::<u64>()
98
            .to_string()
99
    }
100

101
    fn part2(&self) -> String {
102
        let mut dsu = DisjointSet::new(self.nodes.len());
103
        let mut last_edge = (i64::MIN, 0, 0);
104
        let mut components = vec![BTreeSet::new(); self.octree.len()];
105
        while dsu.sizes.len() > 1 {
106

107
            let mut min_edges = dsu.sizes.keys().fold(HashMap::new(), |mut map, &comp| {
108
                map.insert(comp, (i64::MAX, 0, 0));
109
                map
110
            });
111

112
            self.compute_components(&mut dsu, &mut components);
113
            let mut search_stack = vec![(0, 0)];
114
            while let Some((left_node_idx, right_node_idx)) = search_stack.pop() {
115
                let left_node = self.octree[left_node_idx];
116
                let right_node = self.octree[right_node_idx];
117

118
                if components[left_node_idx].len() == 1
119
                    && components[left_node_idx] == components[right_node_idx]
120
                {
121
                    continue;
122
                }
123

124
                let max_min_edge = components[left_node_idx]
125
                    .iter()
126
                    .chain(components[right_node_idx].iter())
127
                    .fold(i64::MIN, |max_edge, &comp| {
128
                        let Some((d, _, _)) = min_edges.get(&comp) else {
129
                            unreachable!("Every component should have an entry in min_edges")
130
                        };
131
                        max_edge.max(*d)
132
                    });
133

134
                if max_min_edge < i64::MAX {
135

136
                    let (left_min, left_max) = self.octree_bounding_boxes[left_node_idx];
137
                    let (right_min, right_max) = self.octree_bounding_boxes[right_node_idx];
138
                    let mut min_possible_dist = 0;
139
                    for dim in 0..3 {
140
                        if left_max[dim] < right_min[dim] {
141
                            min_possible_dist += (right_min[dim] - left_max[dim]).pow(2);
142
                        } else if right_max[dim] < left_min[dim] {
143
                            min_possible_dist += (left_min[dim] - right_max[dim]).pow(2);
144
                        }
145
                    }
146
                    if min_possible_dist >= max_min_edge {
147
                        continue;
148
                    }
149
                }
150

151
                if left_node.size() * right_node.size() <= 64 {
152
                    let pairs = if left_node_idx == right_node_idx {
153

154
                        let nodes = (left_node.start..left_node.end)
155
                            .map(|i| (i, dsu.find(self.nodes[i].index)))
156
                            .collect::<Vec<_>>();
157
                        nodes
158
                            .iter()
159
                            .enumerate()
160
                            .flat_map(|(idx, &i)| {
161
                                nodes[idx + 1..]
162
                                    .iter()
163
                                    .filter_map(move |&j| (i.1 != j.1).then_some((i, j)))
164
                            })
165
                            .collect::<Vec<_>>()
166
                    } else {
167

168
                        let left_nodes = (left_node.start..left_node.end)
169
                            .map(|i| (i, dsu.find(self.nodes[i].index)))
170
                            .collect::<Vec<_>>();
171
                        let right_nodes = (right_node.start..right_node.end)
172
                            .map(|i| (i, dsu.find(self.nodes[i].index)))
173
                            .collect::<Vec<_>>();
174

175
                        left_nodes
176
                            .into_iter()
177
                            .flat_map(|i| {
178
                                right_nodes
179
                                    .iter()
180
                                    .filter_map(move |&j| (i.1 != j.1).then_some((i, j)))
181
                            })
182
                            .collect::<Vec<_>>()
183
                    };
184
                    for ((i, ci), (j, cj)) in pairs {
185
                        let d = Self::dist(&self.nodes[i], &self.nodes[j]);
186
                        min_edges.entry(ci).and_modify(|entry| {
187
                            if d < entry.0 {
188
                                *entry = (d, i, j);
189
                            }
190
                        });
191
                        min_edges.entry(cj).and_modify(|entry| {
192
                            if d < entry.0 {
193
                                *entry = (d, i, j);
194
                            }
195
                        });
196
                    }
197
                    continue;
198
                }
199
                search_stack.extend(self.deeper_search_nodes(left_node_idx, right_node_idx));
200
            }
201
            for (_, (d, i, j)) in min_edges {
202
                dsu.union(self.nodes[i].index, self.nodes[j].index);
203
                if d > last_edge.0 {
204
                    last_edge = (d, i, j);
205
                }
206
            }
207
        }
208
        (self.nodes[last_edge.1].coordinate[0] * self.nodes[last_edge.2].coordinate[0]).to_string()
209
    }
210
}

Summary

The complete code can be found here.

Through the above analysis and implementation, we successfully reduced the solution time to $0.75 ms$ and $24 ms$, achieving a 10-20x performance improvement despite not using parallel computing. This demonstrates that leveraging octree spatial partitioning and branch and bound strategies can effectively reduce computational complexity, allowing us to efficiently handle distance calculation problems between point pairs. Of course, the implementation of this method is relatively complex, requiring careful handling of octree construction, node traversal, and pruning conditions, which requires trade-offs and choices based on specific scenarios.