Advent of Code 2023 day 24

Day 24 was far too fiddly, and used far too much linear algebra, for my taste. To my shame, I resorted to just looking up standard algorithms and other people's solutions to solve this one.

I used the Linear library a lot here.

Representation and parsing

I did my usual thing of creating a new data type for the Hailstone and a parser for it. The only complication here was the symbolP parser that consumed the varying amounts of whitespace.

type Position2 = V2 Rational
type Position3 = V3 Rational

type Segment = (Position2, Position2)

data Hailstone = Hailstone { _pos :: V3 Rational, _vel :: V3 Rational }
  deriving (Show, Eq, Ord)
makeLenses ''Hailstone

stonesP = stoneP `sepBy` endOfLine
stoneP = Hailstone <$> (vertexP <* symbolP "@") <*> vertexP
vertexP = vecify <$> signed decimal <*> (symbolP "," *> signed decimal) <*> (symbolP "," *> signed decimal)
  where vecify x y z = V3 (x % 1) (y % 1) (z % 1)

symbolP :: Text -> Parser Text
symbolP s = (skipSpace *> string s) <* skipSpace

Part 1

This approach used the test for intersecting line segments, given the two end points of the two line segments.

For each hailstone, I found the times when that stone would cross each of the four (infinitely-extended) boundary lines of the test area. For each of those times, and time 0, I filtered out the times that represented points outside the test area. The segment lies between the positions at the minimum and maximum such times. If the hailstone never entered the test area in positive time, there was no segment.

boundaryPoints :: (Position2, Position2) -> Hailstone -> Maybe Segment
boundaryPoints b@((V2 minX minY), (V2 maxX maxY)) h 
  | null ts = Nothing
  | tMax < 0 = Nothing
  | otherwise = Just (V2 xMin yMin, V2 xMax yMax)
  where V2  x  y = h ^. pos . _xy
        V2 vx vy = h ^. vel . _xy
        tMinX = (minX - x) / vx
        tMaxX = (maxX - x) / vx
        tMinY = (minY - y) / vy
        tMaxY = (maxY - y) / vy
        ts = filter withinT $ filter (>= 0) [tMinX, tMinY, tMaxX, tMaxY, 0]
        tMin = minimum ts
        tMax = maximum ts
        xMin = x + tMin * vx
        xMax = x + tMax * vx
        yMin = y + tMin * vy
        yMax = y + tMax * vy
        withinT t = within b (V2 (x + t * vx) (y + t * vy))

within :: (Position2, Position2) -> Position2 -> Bool
within ((V2 minX minY), (V2 maxX maxY)) (V2 x y) = 
  x >= minX && x <= maxX && y >= minY && y <= maxY

Given a pair of line segments, I can detect if they intersect.

intersects :: Segment -> Segment -> Bool
intersects ((V2 x1 y1), (V2 x2 y2)) ((V2 x3 y3), (V2 x4 y4)) = ad && pd
  where a = (x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)
        p = (x1 - x3) * (y1 - y2) - (y1 - y3) * (x1 - x2)
        d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
        ad = if d < 0 then 0 >= a && a >= d
              else 0 <= a && a <= d
        pd = if d < 0 then 0 >= p && p >= d
              else 0 <= p && p <= d

The whole solution finds the segments of all hailstones (using catMaybes to drop the out-of-bounds ones), finds all possible pairs of segments, and counts how many intersect.

part1, part2 :: [Hailstone] -> Int
part1 stones = length $ filter (uncurry intersects) (pairs bps)
  where boundary = (V2 2e14 2e14 , V2 4e14 4e14) -- ((V2 7 7), (V2 27 27))
        bps = catMaybes $ fmap (boundaryPoints boundary) stones

pairs :: Ord a => [a] -> [(a,a)]
pairs xs = [(x,y) | x <- xs, y <- xs, x < y]

Originally, all this used floating-point Double numbers, but I kept getting rounding errors where the boundary detection was off by one part in 10¹⁴. I then realised that I was only doing basic arithmetic with these numbers, so I quickly converted to using the arbitrary-precision Rational number type. It may be fractionally slower, but it's more accurate.

Part 2

Here, I most definitely skipped any solution and just looked up how others' did it. The best explanation I found was in this comment by ReimannIntegirl, referencing this code by Fleur Kelpin.

Hailstone i is described by the pair \( \langle \mathbf{p}_i, \mathbf{v}_i \rangle \) = \( \langle [ x_i, y_i, z_i ], [ \dot x_i, \dot y_i, \dot z_i] \rangle \). The position of hailstone i at time t is \( [ x_i + \dot x_i t, y_i + \dot y_i t, z_i + \dot z_i t ] \).

The hailstones are numbered 1 to 300. The rock we throw is hailstone 0.

Hailstone i and hailstone 0 intersect at some time t when

\begin{eqnarray}
x_i + \dot x_i t & = & x_0 + \dot x_0 t \\
y_i + \dot y_i t & = & y_0 + \dot y_0 t \\
z_i + \dot z_i t & = & z_0 + \dot z_0 t
\end{eqnarray}

Rearranging a bit, you get

\begin{eqnarray}
(x_i - x_0) + (\dot x_i - \dot x_0) t & = & 0 \\
(y_i - y_0) + (\dot y_i - \dot y_0) t & = & 0 \\
(z_i - z_0) + (\dot z_i - \dot z_0) t & = & 0
\end{eqnarray}

which leads to

\[ \frac{x_i - x_0}{\dot x_i - \dot x_0} = \frac{y_i - y_0}{\dot y_i - \dot y_0} = \frac{z_i - z_0}{\dot z_i - \dot z_0} \]

If you take only the x and y dimensions and cross-multiply, you get

\[ (x_i - x_0)(\dot y_i - \dot y_0) = (y_i - y_0)(\dot x_i - \dot x_0) \]

Expanding gives

\[ x_i \dot y_i - x_i \dot y_0 - x_0 \dot y_i + x_0 \dot y_0 = y_i \dot x_i - y_i \dot x_0 - y_0 \dot x_i + y_0 \dot x_0 \]

If you take two such hailstones i and j, and subtract the equations and rearrange, you get

\[ (\dot y_j - \dot y_i) x_0 - (\dot x_j - \dot x_i) y_0 + (y_j - y_i) \dot x_0 + (x_j - x_i) \dot y_0 = y_i \dot x_i - y_j \dot x_j +x_j \dot y_i - x_i \dot y_i \]

If you now take three other pairs hailstones, you will get three other equations. You can treat each equation as a row in the vector equation \(\mathbf{A} \mathbf{x} = \mathbf{b} \) and solve for x. That will give \(x_o, y_o, \dot x_0, \dot y_0 \). If you do the same again, but looking at \(z, \dot z \) rather than (y, \dot y \), you get solutions for all the components of the rock.

That's the theory, now back to the code. The Linear library has a module for Linear.Matrix as well as vectors. It represents a matrix as a vector of vectors. It also has a simple linear equation solver.

coeffs takes a pair of hailstones and a lens (either _y or _z ) and returns the row of the A matrix and the element of the b vector. buildMatrix builds the matrix equation.

buildMatrix :: [Hailstone] -> Lens' Position3 Rational -> (M44 Rational, V4 Rational)
buildMatrix hs l = (V4 a1 a2 a3 a4, V4 b1 b2 b3 b4)
  where ps = take 4 $ chunks2 hs
        [(a1, b1), (a2, b2), (a3, b3), (a4, b4)] = [coeffs p l | p <- ps]

coeffs :: (Hailstone, Hailstone) -> Lens' Position3 Rational -> (V4 Rational, Rational)
coeffs (h1, h2) l = 
  (V4 (h2 ^. vel .  l - h1 ^. vel .  l)
      (h1 ^. vel . _x - h2 ^. vel . _x)
      (h2 ^. pos .  l - h1 ^. pos .  l)
      (h2 ^. pos . _x - h1 ^. pos . _x)
  , h1 ^. pos .  l * h1 ^. vel . _x - 
    h2 ^. pos .  l * h2 ^. vel . _x +
    h2 ^. pos . _x * h2 ^. vel .  l -
    h1 ^. pos . _x * h1 ^. vel .  l
  )

chunks2 :: [a] -> [(a, a)]
chunks2 [] = []
chunks2 [_] = []
chunks2 (x:y:xs) = (x,y) : chunks2 xs

The solution takes these two systems of equations, solves them for the position and velocity of the thrown rock, and builds the solution. There's some fiddling around to convert back from the Rational numbers to a humble Int.

part2 stones = fromIntegral $ numerator solution `div` denominator solution
  where rock = Hailstone (pure 0) (pure 0) -- Hailstone (V3 0 0 0) (V3 0 0 0)
        (ay, by) = buildMatrix stones _y
        (az, bz) = buildMatrix stones _z
        solY = luSolveFinite ay by
        solZ = luSolveFinite az bz
        rock' = rock & pos . _x .~ (solY ^. _x)
                     & pos . _y .~ (solY ^. _y)
                     & pos . _z .~ (solZ ^. _y)
                     & vel . _x .~ (solY ^. _z)
                     & vel . _y .~ (solY ^. _w)
                     & vel . _z .~ (solZ ^. _w)
        solution = (rock' ^. pos . _x) + (rock' ^. pos . _y) + (rock' ^. pos . _z)

And that's it. I'm glad that one's over!

Code

You can get the code from my locally-hosted Git repo, or from Gitlab.