Another factor of 10 in problem size would have required more efficient data structures and the “correct” algorithm. But as it is, we can get by with sets-as-lists and just updating a hash table, if we’re willing to wait a few seconds for the answer.
(ql:quickload :str)
(defun parse-line (line)
(coerce (mapcar #'parse-integer (str:split "," line)) 'vector))
(defun read-inputs (filename)
(let ((input-lines (uiop:read-file-lines filename)))
(make-array (list (length input-lines) 3)
:initial-contents (mapcar #'parse-line input-lines))))
(defun distance (points p q)
(flet ((square (x) (* x x)))
(sqrt (+ (square (- (aref points p 0) (aref points q 0)))
(square (- (aref points p 1) (aref points q 1)))
(square (- (aref points p 2) (aref points q 2)))))))
(defun all-labeled-edges (points)
(loop for j from 0 to (1- (car (array-dimensions points)))
nconcing (loop for i from 0 to (1- j)
collect (list (distance points i j) i j))))
(defun short-labeled-edges (points nedges)
(subseq (sort (all-labeled-edges points) #'< :key #'first) 0 nedges))
(defun adjacency-map (labeled-edges)
(let ((result (make-hash-table)))
(loop for (nil v w) in labeled-edges
do (setf (gethash v result) (cons w (gethash v result)))
(setf (gethash w result) (cons v (gethash w result))))
result))
(defun components (n adjacency-map)
(let ((remaining (loop for i from 0 to (1- n) collect i))
(result nil))
(loop for v = (car remaining)
until (null v)
do (let ((this-component nil)
(next-component (list v)))
(loop until (subsetp next-component this-component)
do (progn
(setf this-component next-component)
(setf next-component
(reduce #'union
(cons this-component
(mapcar #'(lambda (w) (gethash w adjacency-map))
this-component))))))
(setf result (cons this-component result))
(loop for w in this-component
do (setf remaining (delete w remaining)))))
result))
(defun main-1 (filename)
(let* ((points (read-inputs filename))
(adjacency (adjacency-map (short-labeled-edges points 1000)))
(components (components (car (array-dimensions points)) adjacency))
(lengths (sort (mapcar #'length components) #'>)))
(* (car lengths) (cadr lengths) (caddr lengths))))
(defun fusing-edge (n labeled-edges)
(let ((sorted-edges (sort labeled-edges #'< :key #'first))
(components (make-hash-table)))
(loop for i from 0 to (1- n)
do (setf (gethash i components) (list i)))
(loop for (nil v w) in sorted-edges
do (let ((new-component (union (gethash v components) (gethash w components))))
(cond ((= (length new-component) n)
(return (list v w)))
((not (subsetp new-component (gethash v components)))
(loop for x in new-component
do (setf (gethash x components) new-component))))))))
(defun main-2 (filename)
(let* ((points (read-inputs filename))
(n (car (array-dimensions points))))
(destructuring-bind (v w) (fusing-edge n (all-labeled-edges points))
(* (aref points v 0) (aref points w 0)))))


Applied topology! Instead of using the standard even-odd raycasting trick for winding number (which I didn’t remember), I constructed an explicit inward-pointing normal at each edge of the large polygon. This way you can test directly whether edges that lie on the boundary of your candidate rectangle have the right orientation. The code isn’t pretty (
edge-exterior-to-rect?is particularly awful), and I’m sure it could be made faster, but this worked almost the first time.