#include <boost/log/trivial.hpp>
#include <boost/unordered/unordered_flat_set.hpp>
#include <filesystem>
#include <fstream>
#include <stdexcept>

namespace {

struct Point {
  int x, y, z;
};

auto distance(const Point &a, const Point &b) {
  auto &[ax, ay, az] = a;
  auto &[bx, by, bz] = b;
  ssize_t dx = ax - bx;
  ssize_t dy = ay - by;
  ssize_t dz = az - bz;
  return size_t(dx * dx + dy * dy + dz * dz);
}

using PList = std::vector<Point>;
using DMap = std::vector<std::pair<std::pair<size_t, size_t>, size_t>>;
using DSet = boost::unordered::unordered_flat_set<size_t>;

auto read() {
  auto rval = PList{};
  auto ih = std::ifstream{"08.txt"};
  auto line = std::string{};
  while (std::getline(ih, line)) {
    auto c1 = line.find(',');
    auto c2 = line.find(',', c1 + 1);
    rval.emplace_back(
        std::stoi(line.substr(0, c1)),
        std::stoi(line.substr(c1 + 1, c2 - c1 - 1)),
        std::stoi(line.substr(c2 + 1))
    );
  }
  return rval;
}

auto dmap(const PList &t) {
  auto rval = DMap{};
  for (auto a = size_t{}; a < t.size() - 1; ++a)
    for (auto b = a + 1; b < t.size(); ++b)
      rval.emplace_back(std::make_pair(a, b), distance(t.at(a), t.at(b)));
  return rval;
}

auto find(size_t c, std::vector<DSet> &sets) {
  for (auto i = sets.begin(); i != sets.end(); ++i)
    if (i->contains(c))
      return i;
  throw std::runtime_error{"can't find"};
}

auto part1(const PList &t, size_t count, const DMap &d) {
  auto circuits = std::vector<DSet>{};
  for (auto i = size_t{}; i < t.size(); ++i) {
    auto c = DSet{t.size()};
    c.insert(i);
    circuits.push_back(std::move(c));
  }

  for (auto i = size_t{}; i < count; ++i) {
    auto a = d.at(i).first.first;
    auto b = d.at(i).first.second;
    auto ac = find(a, circuits);
    auto bc = find(b, circuits);
    if (ac == bc)
      continue;
    ac->insert(bc->begin(), bc->end());
    circuits.erase(bc);
  }

  std::sort(circuits.begin(), circuits.end(), [](auto &a, auto &b) {
    return a.size() > b.size();
  });
  auto total = size_t{1};
  for (auto i = size_t{}; i < 3; ++i)
    total *= circuits.at(i).size();

  return total;
}

auto part2(const PList &t, const DMap &d) {
  auto circuits = std::vector<DSet>{};
  for (auto i = size_t{}; i < t.size(); ++i) {
    auto c = DSet{t.size()};
    c.insert(i);
    circuits.push_back(std::move(c));
  }

  for (auto i = size_t{}; i < d.size(); ++i) {
    auto a = d.at(i).first.first;
    auto b = d.at(i).first.second;
    auto ac = find(a, circuits);
    auto bc = find(b, circuits);
    if (ac == bc)
      continue;
    ac->insert(bc->begin(), bc->end());
    circuits.erase(bc);
    if (circuits.size() == 1) {
      return ssize_t(t.at(a).x) * ssize_t(t.at(b).x);
    }
  }
  throw std::runtime_error{"never joins?"};
}

} // namespace

auto main() -> int {
  auto t = read();
  BOOST_LOG_TRIVIAL(info) << "Day 8: read " << t.size();
  auto test = std::filesystem::current_path().filename().string() ==
              std::string{"test"};

  auto d = dmap(t);
  std::sort(d.begin(), d.end(), [](auto &a, auto &b) {
    return a.second < b.second;
  });

  BOOST_LOG_TRIVIAL(info) << "1: " << part1(t, test ? 10 : 1000, d);
  BOOST_LOG_TRIVIAL(info) << "2: " << part2(t, d);
}

it produces a textual representation of [][3]u32, try it on your example or input :]

1i [
1,$ {
	s/^/[/
	s/$/]/
}
2,$i,
$i ]
$d

Basic

import "lib/github.com/diku-dk/sorts/radix_sort"

type position = (u32, u32, u32)
def positionFromArray (p: [3]u32): position
  = (p[0], p[1], p[2])
def pair_iota (n: i64): [n](i64, i64)
  = map (\ j -> (n, j)) (iota n)
def gaussian_sum (n: i64) = n * (n + 1) / 2

def euclidean_distance (a: position) (b: position): f64
  = f64.sqrt 
    ( (f64.u32 a.0 - f64.u32 b.0) ** 2
    + (f64.u32 a.1 - f64.u32 b.1) ** 2
    + (f64.u32 a.2 - f64.u32 b.2) ** 2
    )

def distance_table [n] (positions: [n]position): [n][n]f64
  = let distance_function = \ i j -> euclidean_distance positions[i] positions[j] in
    tabulate_2d n n distance_function

def existsLength 'a 'b (f: a -> ?[l].[l]b) (x: a): i64
  = length (f x)

def concatMap [n] 'a 'b (f: a -> ?[l].[l]b) (placeholder: b) (xs: [n]a): *[]b
  = let totalLength = reduce (+) 0 <| map (\ x -> length (f x)) xs in
    ( loop (results, offset) = (replicate totalLength placeholder, 0)
      for x in xs
      do
        let bs = f x in
        let scatterIndices = indices bs |> map (+offset) in
        (scatter results scatterIndices bs, offset + length bs)
    ).0

def distance_array [n] (positions: [n]position): []((i64, i64), f64)
  = let table = distance_table positions in
    let triangle_indices = concatMap pair_iota (i64.highest, i64.highest) (iota n |> drop 1) in
    map (\ (i, j) -> ((i, j), table[i, j])) triangle_indices

def sort_distances (distances: []((i64, i64), f64)): []((i64, i64), f64)
  = radix_sort_float_by_key (.1) f64.num_bits f64.get_bit distances

type option 'a
  = #Empty
  | #Present a

def empty 'a : option a = #Empty

def overrideWith (old: u16) (new: u16) (x: option u16): option u16
  = match x
      case #Empty -> #Empty
      case #Present inner -> 
        if inner == old
        then #Present new
        else #Present inner

def orElse 'a (o: option a) (d: a): a
  = match o
      case #Empty -> d
      case #Present x -> x

def is_present 'a (o: option a): bool
  = match o
      case #Empty -> false
      case #Present _ -> true

def connect (circuits: *[](option u16)) (newCircuitId: u16) (connection: (i64, i64)): (u16, *[](option u16))
  = let circuitA = circuits[connection.0] in
    let circuitB = circuits[connection.1] in
    match (circuitA, circuitB)
      case (#Empty, #Empty) -> 
        ( newCircuitId + 1
        , scatter circuits [connection.0, connection.1] (rep (#Present newCircuitId))
        )
      case (#Present a, #Empty) -> 
        ( newCircuitId
        , scatter circuits [connection.1] [#Present a]
        )
      case (#Empty, #Present b) -> 
        ( newCircuitId
        , scatter circuits [connection.0] [#Present b]
        )
      case (#Present a, #Present b) ->
        ( newCircuitId
        , map (b `overrideWith` a) circuits
        )

def countCircuit (counts: *[]u64) (o: option u16): *[]u64 
  = match o
    case #Empty -> counts 
    case #Present i -> scatter counts [i64.u16 i] [counts[i64.u16 i] + 1]

def countCircuits (c: u16) (circuits: [](option u16)): *[i64.u16 c]u64
  = let circuitCounts = replicate (i64.u16 c) 0 in
    loop counts = circuitCounts
    for circuit in circuits
    do countCircuit counts circuit

def exampleConnectionCount = 10i64
def inputConnectionCount = 1000i64

def part1 (positions: i64) (connectionCount: i64) (distances: []((i64, i64), f64))
  = let connections = take connectionCount distances |> map (.0) in
    let circuitMap: *[positions](option u16) = replicate positions empty in
    ( loop (circuitCount, circuits) = (0, circuitMap)
      for connection in connections
      do
        connect circuits circuitCount connection
    ) |> uncurry countCircuits 
      |> radix_sort u64.num_bits u64.get_bit
      |> reverse
      |> take 3
      |> foldl (*) 1

def part2 (positionCount: i64) (distances: []((i64, i64), f64)) (positions: []position)
  = let circuitMap: *[positionCount](option u16) = replicate positionCount empty in
    ( loop (circuitCount, connectionIndex, circuits) = (0, 0, circuitMap)
      while not
        ( and (map is_present circuits)
        && and (map (== circuits[0]) circuits)
        )
      do
        let connection = distances[connectionIndex].0 in
        let (newCircuitId, circuits') = connect circuits circuitCount connection in
        (newCircuitId, connectionIndex+1, circuits')
    ).1
    |> \ i -> distances[i-1].0
    |> \ (a, b) -> positions[a].0 * positions[b].0

def main [n] (position_array: [n][3]u32)
  = let positions = map positionFromArray position_array in
    let unsorted_distances = distance_array positions in
    let sorted_distances = sort_distances unsorted_distances in
    ( part1 n inputConnectionCount sorted_distances
    , part2 n sorted_distances positions
    )

#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <assert.h>

#define LEN(a)		(sizeof(a)/sizeof(*(a)))
#define NPAIRS(x)	((x)*((x)-1)/2)

#define MAXN	1024

struct node { int x,y,z, group; };
struct pair { struct node *a, *b; int64_t dist_sq; };
struct group { int count; };

static struct node nodes[MAXN];
static struct pair pairs[NPAIRS(MAXN)];
static struct group groups[MAXN];

static int nnodes;
static int ngroups;

static int64_t
node_dist_sq(const struct node *a, const struct node *b)
{
	return
	    (int64_t)(a->x - b->x) * (a->x - b->x) +
	    (int64_t)(a->y - b->y) * (a->y - b->y) +
	    (int64_t)(a->z - b->z) * (a->z - b->z);
}

static int
cmp_pairs(const void *va, const void *vb)
{
	const struct pair *a = va;
	const struct pair *b = vb;

	return
	    a->dist_sq < b->dist_sq ? -1 :
	    a->dist_sq > b->dist_sq ?  1 : 0;
}

static int
cmp_groups_asc(const void *va, const void *vb)
{
	const struct group *a = va;
	const struct group *b = vb;

	return b->count - a->count;
}

static void
merge_groups(int group_a, int group_b)
{
	int i;

	if (group_a == group_b)
		return;

	groups[group_a].count += groups[group_b].count;
	groups[group_b].count = 0;

	for (i=0; i<nnodes; i++)
		if (nodes[i].group == group_b)
			nodes[i].group = group_a;
	
	ngroups--;
}

int
main()
{
	int p1=0,p2=0, p1_limit, i,j, n,p;

	for (; ; nnodes++) {
		assert(nnodes < MAXN);
		n = scanf(" %d,%d,%d",
		    &nodes[nnodes].x,
		    &nodes[nnodes].y,
		    &nodes[nnodes].z);
		if (n < 3)
			break;
		nodes[nnodes].group = nnodes;
		groups[nnodes].count = 1;
	}

	ngroups = nnodes;

	for (p=0, i=0; i<nnodes-1; i++)
	for (j=i+1; j<nnodes; j++, p++) {
		pairs[p].a = &nodes[i];
		pairs[p].b = &nodes[j];
		pairs[p].dist_sq = node_dist_sq(&nodes[i], &nodes[j]);
	}

	qsort(pairs, NPAIRS(nnodes), sizeof(*pairs), cmp_pairs);

	p1_limit = nnodes <= 100 ? 10 : 1000;

	for (i=0; ngroups > 1; i++) {
		merge_groups(pairs[i].a->group, pairs[i].b->group);

		if (ngroups == 1)
			p2 = pairs[i].a->x * pairs[i].b->x;

		if (i == p1_limit) {
			qsort(groups, LEN(groups), sizeof(*groups),
			    cmp_groups_asc);

			p1 = groups[0].count *
			     groups[1].count *
			     groups[2].count;
		}
	}

	printf("08: %d %d\n", p1, p2);
}

type Id = usize;
type Connection = (Id, Id);
type Circuit = HashSet<Id>;

#[derive(PartialEq)]
struct Point {
    x: usize,
    y: usize,
    z: usize,
}

impl FromStr for Point {
    type Err = Report;

    fn from_str(s: &str) -> Result<Self> {
        let mut parts = s.split(',');
        Ok(Point {
            x: parts.next().ok_or_eyre("missing x")?.parse()?,
            y: parts.next().ok_or_eyre("missing y")?.parse()?,
            z: parts.next().ok_or_eyre("missing z")?.parse()?,
        })
    }
}

impl Point {
    fn distance(&self, other: &Point) -> R64 {
        let dist = ((
            self.x.abs_diff(other.x).pow(2) +
            self.y.abs_diff(other.y).pow(2) +
            self.z.abs_diff(other.z).pow(2)) as f64)
            .sqrt();
        r64(dist)
    }
}

class Day08 : Puzzle {

    lateinit var boxes: List<Point3D>
    lateinit var edges: List<Pair<List<Point3D>, Double>>

    override fun readFile() {
        val input = readInputFromFile(2025, 8, false)
        boxes = input.lines().filter { it.isNotBlank() }
            .map { it.split(",") }
            .map { Point3D(it[0].toInt(), it[1].toInt(), it[2].toInt()) }
        edges = calculateEdges(boxes)
    }

    private fun calculateEdges(boxes: List<Point3D>): List<Pair<List<Point3D>, Double>> {
        val edges = mutableListOf<Pair<MutableList<Point3D>, Double>>()
        for (i in boxes.indices) {
            for (j in i + 1 until boxes.size) {
                edges.add(Pair(mutableListOf(boxes[i], boxes[j]), boxes[i].dist(boxes[j])))
            }
        }
        return edges.sortedBy { it.second }
    }

    override fun solvePartOne(): String {
        val connections = buildEmptyConnectionMatrix(boxes.size)
        val mutableEdges = edges.toMutableList()
        for (i in 0..<1000) {
            connectNextEdge(mutableEdges, connections)
        }

        val connectionSet = buildConnectionSet(boxes, connections)

        return connectionSet.map { it.size }
            .sortedByDescending { it }
            .take(3)
            .reduce(Int::times)
            .toString()
    }

    override fun solvePartTwo(): String {
        val connections = buildEmptyConnectionMatrix(boxes.size)
        val mutableEdges = edges.toMutableList()

        var result: Long? = null
        while (result == null) {
            result = connectNextEdge(mutableEdges, connections, true)
        }

        return result.toString()
    }

    // size: width & height of (square) matrix
    private fun buildEmptyConnectionMatrix(size: Int): Array<Array<Boolean>> {
        val connections = Array(size) { Array(size) { false } }
        for (i in connections.indices) {
            connections[i][i] = true
        }
        return connections
    }

    private fun connectNextEdge(mutableEdges: MutableList<Pair<List<Point3D>, Double>>, connections: Array<Array<Boolean>>, part2: Boolean = false): Long? {
        if (mutableEdges.isEmpty()) return null
        val next = mutableEdges[0]

        val point = next.first[0]
        val other = next.first[1]
        connectAll(boxes.indexOf(point), boxes.indexOf(other), connections)
        mutableEdges.remove(next)

        // all nodes are connected, assume that this is happening for the first time
        return if (part2 && connections[0].all { it }) {
            next.first[0].x.toLong() * next.first[1].x.toLong()
        } else {
            null
        }
    }

    private fun connectAll(index: Int, other: Int, connections: Array<Array<Boolean>>) {
        fun connectHelper(hIndex: Int) {
            val newConnections = mutableSetOf<Int>()
            for (i in connections[hIndex].indices) {
                if (connections[hIndex][i]) newConnections.add(i)
            }
            for (boxIndex in newConnections.filter { it != hIndex }) {
                for (conn in newConnections.filter { it != boxIndex }) {
                    connections[boxIndex][conn] = true
                }
            }
        }
        connections[index][other] = true

        connectHelper(index) // update matrix with all values from node at [index]
        connectHelper(other) // update matrix with all values from node at [other]
    }

    // returns 2D-list of all indices of currently active connections
    private fun buildConnectionSet(boxes: List<Point3D>, connections: Array<Array<Boolean>>): Set<List<Int>> {
        val connectionSet = mutableSetOf<List<Int>>()
        val indices = (0 until boxes.size).toMutableList()
        while (indices.isNotEmpty()) {
            val list = mutableListOf<Int>()
            val curr = indices.removeFirst()
            val array = connections[curr]
            for (j in array.indices) {
                if (array[j]) list.add(j)
            }
            connectionSet.add(list)
        }
        return connectionSet
    }

    data class Point3D(val x: Int, val y: Int, val z: Int) {
        fun dist(other: Point3D): Double {
            return sqrt(
                (x - other.x).toDouble().pow(2) +
                        (y - other.y).toDouble().pow(2) +
                        (z - other.z).toDouble().pow(2)
            )
        }
    }
}

package main

import (
	"aoc/utils"
	"cmp"
	"math"
	"slices"
	"strconv"
	"strings"
)

type pos3D struct {
	x, y, z int
}

func (p pos3D) Compare(other pos3D) int {
	dx := cmp.Compare(p.x, other.x)
	if dx != 0 {
		return dx
	}

	dy := cmp.Compare(p.y, other.y)
	if dy != 0 {
		return dy
	}

	return cmp.Compare(p.z, other.z)
}

func (p pos3D) distance(other pos3D) float64 {
	dx := float64(other.x - p.x)
	dy := float64(other.y - p.y)
	dz := float64(other.z - p.z)

	d2 := math.Pow(dx, 2.0) + math.Pow(dy, 2.0) + math.Pow(dz, 2.0)
	return math.Sqrt(d2)
}

func getPosChannel(input chan string) chan pos3D {
	ch := make(chan pos3D, cap(input))

	go func() {
		for line := range input {
			parts := strings.Split(line, ",")
			x, _ := strconv.Atoi(parts[0])
			y, _ := strconv.Atoi(parts[1])
			z, _ := strconv.Atoi(parts[2])
			ch <- pos3D{x, y, z}
		}
		close(ch)
	}()

	return ch
}

type circuits struct {
	circuits map[pos3D]int
	nextID   int
}

func (cc *circuits) newCircuit() (id int) {
	id = cc.nextID
	cc.nextID++
	return id
}

func (cc *circuits) mergeCircuits(id1, id2 int) {
	for p, id := range cc.circuits {
		if id == id2 {
			cc.circuits[p] = id1
		}
	}
}

func createSingletonCircuits(points []pos3D) (cc circuits) {
	cc.circuits = make(map[pos3D]int)
	for _, p := range points {
		id := cc.newCircuit()
		cc.circuits[p] = id
	}
	return cc
}

func (cc circuits) reverseMap() (m map[int][]pos3D) {
	m = make(map[int][]pos3D)

	for p, id := range cc.circuits {
		if _, ok := m[id]; !ok {
			m[id] = []pos3D{}
		}

		m[id] = append(m[id], p)
	}

	return m
}

func (cc circuits) sizeMap() map[int]int {
	circuitSizeMap := make(map[int]int)
	for _, id := range cc.circuits {
		circuitSizeMap[id]++
	}
	return circuitSizeMap
}

type targetedDistance struct {
	distance float64
	p1, p2   pos3D
}

func (p pos3D) distanceWithAll(
	points []pos3D,
	alreadyVisited *map[pos3D]any,
) (tds []targetedDistance) {
	tds = []targetedDistance{}
	for _, op := range points {
		if _, ok := (*alreadyVisited)[op]; ok {
			continue
		}

		var td targetedDistance
		td.distance = math.MaxFloat64
		td.p1 = p
		d := p.distance(op)
		if d < td.distance {
			td.distance = d
			td.p2 = op
		}
		tds = append(tds, td)
	}
	return tds
}

type pointPair [2]pos3D

func (pp pointPair) equals(other pointPair) bool {
	pp1 := pp[0]
	pp2 := pp[1]
	o1 := other[0]
	o2 := other[1]
	return (pp1 == o1 && pp2 == o2) || (pp1 == o2 && pp2 == o1)
}

var IterationCount = 1000

func stepOne(input chan string) (int, error) {
	ch := getPosChannel(input)

	points := []pos3D{}
	for p := range ch {
		points = append(points, p)
	}
	slices.SortFunc(points, pos3D.Compare)
	cc := createSingletonCircuits(points)

	alreadyVisited := make(map[pos3D]any)
	tds := []targetedDistance{}
	for _, p := range points {
		alreadyVisited[p] = nil
		dsts := p.distanceWithAll(points, &alreadyVisited)
		tds = append(tds, dsts...)
	}

	slices.SortFunc(tds, func(a, b targetedDistance) int {
		return cmp.Compare(a.distance, b.distance)
	})

	for idx := range IterationCount {
		td := tds[idx]
		cc.mergeCircuits(cc.circuits[td.p1], cc.circuits[td.p2])
	}

	circuitSizeMap := cc.sizeMap()

	circuitSizes := []int{}
	for _, v := range circuitSizeMap {
		circuitSizes = append(circuitSizes, v)
	}
	slices.Sort(circuitSizes)
	largestThree := circuitSizes[len(circuitSizes)-3:]

	product := 1
	for _, v := range largestThree {
		product *= v
	}

	return product, nil
}

func stepTwo(input chan string) (int, error) {
	ch := getPosChannel(input)

	points := []pos3D{}
	for p := range ch {
		points = append(points, p)
	}
	slices.SortFunc(points, pos3D.Compare)
	cc := createSingletonCircuits(points)

	alreadyVisited := make(map[pos3D]any)
	tds := []targetedDistance{}
	for _, p := range points {
		alreadyVisited[p] = nil
		dsts := p.distanceWithAll(points, &alreadyVisited)
		tds = append(tds, dsts...)
	}

	slices.SortFunc(tds, func(a, b targetedDistance) int {
		return cmp.Compare(a.distance, b.distance)
	})

	idx := 0
	var lastConnection pointPair

	for {
		td := tds[idx]
		idx++
		cc.mergeCircuits(cc.circuits[td.p1], cc.circuits[td.p2])

		circuitSizeMap := cc.sizeMap()
		if len(circuitSizeMap) == 1 {
			lastConnection = pointPair{td.p1, td.p2}
			break
		}
	}

	return lastConnection[0].x * lastConnection[1].x, nil
}

func main() {
	inputFile := utils.FilePath("day08.txt")
	utils.RunStep(utils.ONE, inputFile, stepOne)
	utils.RunStep(utils.TWO, inputFile, stepTwo)
}

type
  AOCSolution[T,U] = tuple[part1: T, part2: U]
  Vec3 = tuple[x,y,z: int]
  Node = ref object
    pos: Vec3
    cid: int

proc dist(a,b: Vec3): float =
  sqrt(float((a.x-b.x)^2 + (a.y-b.y)^2 + (a.z-b.z)^2))

proc solve(input: string, p1_limit: int): AOCSolution[int, int] =
  let boxes = input.splitLines().mapIt:
    let parts = it.split(',')
    let pos = Vec3 (parseInt parts[0], parseInt parts[1], parseInt parts[2])
    Node(pos: pos, cid: -1)

  var dists: seq[(float, (Node, Node))]
  for i in 0 .. boxes.high - 1:
    for j in i+1 .. boxes.high:
      dists.add (dist(boxes[i].pos, boxes[j].pos), (boxes[i], boxes[j]))

  var curcuits: Table[int, HashSet[Node]]
  var curcuitID = 0

  dists.sort(cmp = proc(a,b: (float, (Node, Node))): int = cmp(a[0], b[0]))

  for ind, (d, nodes) in dists:
    var (a, b) = nodes
    let (acid, bcid) = (a.cid, b.cid)

    if acid == -1 and bcid == -1: # new curcuit
      a.cid = curcuitID
      b.cid = curcuitID
      curcuits[curcuitId] = [a, b].toHashSet
      inc curcuitID
    elif bcid == -1: # add to a
      b.cid = acid
      curcuits[acid].incl b
    elif acid == -1: # add to b
      a.cid = bcid
      curcuits[bcid].incl a
    elif acid != bcid: # merge two curcuits
      for node in curcuits[bcid]:
        node.cid = acid
      curcuits[acid].incl curcuits[bcid]
      curcuits.del bcid

    if ind+1 == p1_limit:
      result.part1 = curcuits.values.toseq.map(len).sorted()[^3..^1].prod

    if not(acid == bcid and acid != -1): result.part2 = a.pos.x * b.pos.x

use euclid::default::Point3D;
use euclid::point3;

fn parse_input(input: &str) -> Vec<Point3D<i64>> {
    input
        .lines()
        .map(|l| {
            let mut parts = l.split(',').map(|p| p.parse::<i64>().unwrap());
            let (x, y, z) = (
                parts.next().unwrap(),
                parts.next().unwrap(),
                parts.next().unwrap(),
            );
            point3(x, y, z)
        })
        .collect()
}

// Distances between all points. Reflexive and symmetric pairs are skipped,
// so the Vec's have increasing size, starting at 0.
fn dists(points: &[Point3D<i64>]) -> Vec<Vec<i64>> {
    points
        .iter()
        .enumerate()
        .map(|(idx, &p1)| {
            points
                .iter()
                .take(idx)
                .map(|&p2| (p2 - p1).square_length())
                .collect::<Vec<i64>>()
        })
        .collect()
}

fn sorted_distances(dists: &[Vec<i64>]) -> Vec<(usize, usize, i64)> {
    let mut sorted: Vec<(usize, usize, i64)> = dists
        .iter()
        .enumerate()
        .flat_map(|(i, row)| row.iter().enumerate().map(move |(j, d)| (i, j, *d)))
        .collect();
    sorted.sort_by_key(|(_, _, d)| *d);
    sorted
}

fn part1(input: String) {
    let points = parse_input(&input);
    let d = dists(&points);
    let sorted = sorted_distances(&d);

    let mut circuits: Vec<u32> = (0..points.len() as u32).collect();
    for (i, j, _) in sorted.into_iter().take(1000) {
        let new_circuit = circuits[i];
        let old_circuit = circuits[j];
        if new_circuit != old_circuit {
            for c in circuits.iter_mut() {
                if *c == old_circuit {
                    *c = new_circuit;
                }
            }
        }
    }
    let mut sizes: Vec<u32> = vec![0; points.len()];
    for c in circuits {
        sizes[c as usize] += 1
    }
    sizes.sort_unstable();
    let result = sizes.iter().rev().take(3).product::<u32>();
    println!("{result}");
}

fn part2(input: String) {
    let points = parse_input(&input);
    let d = dists(&points);
    let sorted = sorted_distances(&d);

    let mut circuits: Vec<u32> = (0..points.len() as u32).collect();
    for (i, j, _) in sorted.into_iter() {
        let new_circuit = circuits[i];
        let old_circuit = circuits[j];
        if new_circuit != old_circuit {
            let mut all_connected = true;
            for c in circuits.iter_mut() {
                if *c == old_circuit {
                    *c = new_circuit;
                }
                if *c != new_circuit {
                    all_connected = false;
                }
            }
            if all_connected {
                let result = points[i].x * points[j].x;
                println!("{result}");
                return;
            }
        }
    }
}

util::aoc_main!();

I added elements just until the set first hit size 1, rather than continuing until total number of elements equalled number of points.

import heapq as hq
from itertools import combinations
from collections import Counter

# Disjoint Set Union (or Union-Find) data structure
#   with path compression and union by rank
class DSU:
    def __init__(self, size: int) -> None:
        self.parent = list(range(size)) # ith value is the parent node to node i
        self.rank = [1] * size          # max depth of subtree rooted here (used for union by rank)
        
    def find(self, x: int) -> int:
        # if the node is not its own parent, we need to recurse on parent
        if x != self.parent[x]:
            # path compression
            self.parent[x] = self.find(self.parent[x])
        return self.parent[x]
    
    def isConnected(self, x: int, y: int) -> bool:
        return self.find(x) == self.find(y)
    
    # returns a boolean whether or not union was needed
    def union(self, x: int, y: int) -> bool:
        rootX = self.find(x)
        rootY = self.find(y)        
        
        if rootX == rootY:
            # no union needed
            return False
        
        if self.rank[rootX] > self.rank[rootY]:
            # rootX has deeper subtree, therefore set it as parent to rootY (and its subtree)
            self.parent[rootY] = rootX
        elif self.rank[rootX] < self.rank[rootY]:
            # rootY has deeper subtree, therefore set it as parent to rootX (and its subtree)
            self.parent[rootX] = rootY
        else:
            # both subtrees are of equal depth, therefore choose either one of them as the parent of the other
            # here we chose rootX as the parent of rootY, therefore rootX depth increases by 1
            self.parent[rootY] = rootX
            self.rank[rootX] += 1
        
        # union complete
        return True

# parse puzzle input into junction coordinates (x, y, z)
def parse_junctions(data: str):
    return [tuple(map(int, line.split(','))) for line in data.splitlines()]

# get distance between two points in 3D space
# skipping the sqrt because we only care about relative distances
def get_dist(p1, p2):
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    return (x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2

def part1(data: str, max_connections: int = 1000):
    junctions = parse_junctions(data)
    n = len(junctions)

    # <max_connections> shortest connections
    closest_connections = []
    # iterate over all pairs of junctions
    # combinations is just some syntactic sugar to avoid a nested for loop (still O(n^2))
    for i, j in combinations(range(n), 2):
        dist = get_dist(junctions[i], junctions[j])
        # keep <max_connections> closest connections
        #   have to use -dist to simulate max heap
        if len(closest_connections) < max_connections:
            hq.heappush(closest_connections, (-dist, i, j))
        else:
            hq.heappushpop(closest_connections, (-dist, i, j))

    # union all the closest connections
    dsu = DSU(n)
    for _, i, j in closest_connections:
        dsu.union(i, j)

    # count all circuit sizes
    circuit_sizes = Counter()
    for i in range(n):
        circuit_sizes[dsu.find(i)] += 1

    # multiply the sizes of the 3 largest circuits
    ans = 1
    for _, f in circuit_sizes.most_common(3):
        ans *= f
    return ans

def part2(data: str):
    junctions = parse_junctions(data)
    n = len(junctions)

    # all n^2 junction connections
    all_connections = []
    # combinations is just some syntactic sugar to avoid a nested for loop (still O(n^2))
    for i, j in combinations(range(n), 2):
        dist = get_dist(junctions[i], junctions[j])
        all_connections.append((dist, i, j))
    
    # create min heap of all connections
    # Heapify later to benefit from the heapify algorithm (O(n) instead of O(n log n))
    hq.heapify(all_connections)
    
    # connect junctions until all are connected
    dsu = DSU(n)
    unions = 0
    while all_connections:
        _, i, j = hq.heappop(all_connections)
        unions += dsu.union(i, j)
        # if we have n-1 unions, that means all n junctions are connected
        if unions == n-1:
            return junctions[i][0] * junctions[j][0]
    
    assert False, "It's not possible that all junctions never connect"

sample = """162,817,812
57,618,57
906,360,560
592,479,940
352,342,300
466,668,158
542,29,236
431,825,988
739,650,466
52,470,668
216,146,977
819,987,18
117,168,530
805,96,715
346,949,466
970,615,88
941,993,340
862,61,35
984,92,344
425,690,689"""
assert part1(sample, max_connections=10) == 40
assert part2(sample) == 25272

    #[derive(Debug)]
    struct Circuits {
        circuits: Vec<Vec<usize>>,
    }

    impl Circuits {
        fn find_circuit(&mut self, i: usize) -> Option<usize> {
            for (j, c) in &mut self.circuits.iter().enumerate() {
                if c.contains(&i) {
                    return Some(j);
                }
            }
            None
        }

        fn new_circuit(&mut self, i: usize, j: usize) {
            self.circuits.push(vec![i, j]);
        }

        fn join_circuits(&mut self, i: usize, j: usize) {
            let mut other = self.circuits.get(j).unwrap().clone();
            self.circuits.get_mut(i).unwrap().append(&mut other);
            self.circuits.remove(j);
        }

        fn append_circuit(&mut self, c: usize, x: usize) {
            self.circuits.get_mut(c).unwrap().push(x);
        }

        fn top_three(&mut self) -> Vec<usize> {
            let mut sizes = self
                .circuits
                .iter()
                .map(|c| c.len())
                .collect::<Vec<usize>>();
            sizes.sort();
            sizes.reverse();
            sizes[0..3].to_vec()
        }
    }

    type Pole = (isize, isize, isize);

    fn dist_poles(p1: &Pole, p2: &Pole) -> f32 {
        sqrt(
            ((p1.0 - p2.0) * (p1.0 - p2.0)) as f32
                + ((p1.1 - p2.1) * (p1.1 - p2.1)) as f32
                + ((p1.2 - p2.2) * (p1.2 - p2.2)) as f32,
        )
    }

    #[test]
    fn test_y2025_day8_part1() {
        let input = include_str!("../../input/2025/day_8.txt");
        let poles = input
            .lines()
            .map(|l| {
                let [x, y, z] = l
                    .splitn(3, ",")
                    .map(|c| c.parse::<isize>().unwrap())
                    .collect::<Vec<isize>>()[..]
                else {
                    panic!();
                };
                (x, y, z)
            })
            .collect::<Vec<Pole>>();
        let len = poles.len();

        let mut circuits: Circuits = Circuits { circuits: vec![] };

        let mut pairs = vec![];

        for i in 0..len {
            let first = poles.get(i).unwrap();
            for j in i + 1..len {
                if i == j {
                    continue;
                }
                let second = poles.get(j).unwrap();
                let dist = dist_poles(first, second);
                pairs.push((dist, i, j));
            }
        }

        pairs.sort_by(|a, b| {
            if a.0 < b.0 {
                Ordering::Less
            } else {
                Ordering::Greater
            }
        });

        for (dist, a, b) in pairs[0..1000].iter() {
            let first_circuit = circuits.find_circuit(*a);
            let second_circuit = circuits.find_circuit(*b);

            match (first_circuit, second_circuit) {
                (None, None) => {
                    circuits.new_circuit(*a, *b);
                }
                (Some(c), None) => {
                    circuits.append_circuit(c, *b);
                }
                (None, Some(c)) => {
                    circuits.append_circuit(c, *a);
                }
                (Some(c1), Some(c2)) => {
                    if c1 != c2 {
                        circuits.join_circuits(c1, c2);
                    }
                }
            }
        }

        assert_eq!(circuits.top_three().iter().product::<usize>(), 66640)
    }