-
Notifications
You must be signed in to change notification settings - Fork 0
/
Day24.hs
233 lines (192 loc) · 10.3 KB
/
Day24.hs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
{-# LANGUAGE OverloadedRecordDot, DuplicateRecordFields, UndecidableInstances, TypeApplications, NoFieldSelectors #-}
import System.Environment (getArgs)
import System.Exit (exitFailure)
import System.IO (readFile')
import Data.List (dropWhileEnd)
import Data.Char (isSpace)
import Data.Maybe (mapMaybe)
import Control.Monad (forM_)
import Data.Either (fromRight)
-- | Performs an assertion within the Either monad.
assert :: Bool -> e -> Either e ()
assert True _ = Right ()
assert False e = Left e
-- | Unwraps an Either with a string error.
fromRight' :: Show a => Either a b -> b
fromRight' (Left e) = error $ show e
fromRight' (Right x) = x
-- | Trims whitespace from the start/end.
trim :: String -> String
trim = dropWhileEnd isSpace . dropWhile isSpace
-- | Splits on the given value.
split :: Eq a => a -> [a] -> [[a]]
split x (y:ys) | x == y = [] : split x ys
| otherwise = let (y':ys') = split x ys
in (y : y') : ys'
split _ [] = [[]]
-- | Finds all unordered pairs.
pairs :: [a] -> [(a, a)]
pairs (x:xs) = ((x,) <$> xs) ++ pairs xs
pairs [] = []
data Vec2 a = Vec2 { x :: a, y :: a }
deriving (Show, Eq, Functor)
data Vec3 a = Vec3 { x :: a, y :: a, z :: a }
deriving (Show, Eq, Functor)
data Rect2 a = Rect2 { topLeft :: Vec2 a, bottomRight :: Vec2 a }
deriving (Show, Eq, Functor)
-- | A linear system of equations of the form ax = b
data LinearSystem a = LinearSystem { a :: [[a]], b :: [a] }
deriving (Show, Eq, Functor)
data Hailstone a = Hailstone { pos :: a, vel :: a }
deriving (Show, Eq, Functor)
type Hailstone2 = Hailstone (Vec2 Double)
type Hailstone3 = Hailstone (Vec3 Double)
class Functor v => Vec v where
(.+.) :: Num a => v a -> v a -> v a
neg :: Num a => v a -> v a
neg = fmap negate
(.-.) :: Num a => v a -> v a -> v a
(.-.) v w = v .+. (neg w)
(*.) :: Num a => a -> v a -> v a
(*.) l = fmap (l *)
(.*) :: Num a => v a -> a -> v a
(.*) = flip (*.)
(./) :: Fractional a => v a -> a -> v a
(./) v l = fmap (/ l) v
instance Vec Vec2 where
v .+. w = Vec2 (v.x + w.x) (v.y + w.y)
instance Vec Vec3 where
v .+. w = Vec3 (v.x + w.x) (v.y + w.y) (v.z + w.z)
instance Vec [] where
(.+.) = zipWith (+)
(.-.) = zipWith (-)
-- | Creates a vector with both components set to the given value.
both :: a -> Vec2 a
both x = Vec2 x x
-- | Projects the given Vec3 to the xy-plane, i.e. a Vec2.
projectXY :: Vec3 a -> Vec2 a
projectXY v = Vec2 v.x v.y
-- | Checks whether the given rectangle contains the given point.
inRect :: Ord a => Rect2 a -> Vec2 a -> Bool
inRect r v = v.x >= r.topLeft.x && v.x <= r.bottomRight.x
&& v.y >= r.topLeft.y && v.y <= r.bottomRight.y
-- | Computes the 2x2 determinant from the given columns.
det2 :: Vec2 Double -> Vec2 Double -> Double
det2 i j = i.x * j.y - i.y * j.x
-- | Computes the 2D intersection between the given two hailstones using Cramer's rule.
-- See https://math.stackexchange.com/questions/406864/intersection-of-two-lines-in-vector-form
intersect2 :: Hailstone2 -> Hailstone2 -> Maybe (Vec2 Double)
intersect2 a b | abs d > 0.00001 && ta >= 0 && tb >= 0 = Just (a.pos .+. (ta *. a.vel))
| otherwise = Nothing
where
d = det2 a.vel (neg b.vel)
dp = b.pos .-. a.pos
ta = det2 dp (neg b.vel) / d
tb = det2 a.vel dp / d
-- | Transforms a linear system to row-echolon form.
rowEcholonForm :: Fractional a => LinearSystem a -> LinearSystem a
rowEcholonForm sys = rowEcholonForm' 0 sys
where rowEcholonForm' i sys | i < n = let (ra:ras) = sys.a
(rb:rbs) = sys.b
rx = ra !! i
ra' = ra ./ rx
rb' = rb / rx
rxs = (!! i) <$> ras
ras' = zipWith (\r rx -> r .-. (ra' .* rx)) ras rxs
rbs' = zipWith (\r rx -> r - (rb' * rx)) rbs rxs
sys' = rowEcholonForm' (i + 1) $ LinearSystem ras' rbs'
in LinearSystem (ra':sys'.a) (rb':sys'.b)
| otherwise = sys
n = length sys.a
-- | Flips the matrix/vector horizontally and vertically.
flipLinearSystem :: LinearSystem a -> LinearSystem a
flipLinearSystem sys = LinearSystem (reverse (reverse <$> sys.a)) (reverse sys.b)
-- | Trasnforms a linear system to reduced row echolon form.
reducedRowEcholonForm :: Fractional a => LinearSystem a -> LinearSystem a
reducedRowEcholonForm = flipLinearSystem . rowEcholonForm . flipLinearSystem . rowEcholonForm
-- | Solves a linear system using Gaussian elimination.
solveLinearSystem :: Fractional a => LinearSystem a -> Either String [a]
solveLinearSystem sys = do
-- Input validation
let n = length sys.a
assert (length sys.b == n) $ "Matrix a must have same numbers of rows as vector b"
forM_ sys.a $ \row -> do
assert (length row == n) "Matrix a must be quadratic"
let rref = reducedRowEcholonForm sys
Right rref.b
-- | Parsea a Vec3 from the given string.
parseVec3 :: String -> Maybe (Vec3 Double)
parseVec3 raw = do
[x, y, z] <- Just $ read . trim <$> split ',' raw
Just $ Vec3 x y z
-- | Parses a hailstone from the given string.
parseHailstone :: String -> Maybe Hailstone3
parseHailstone raw = do
[pos, vel] <- mapM parseVec3 . map trim $ split '@' raw
Just $ Hailstone pos vel
main :: IO ()
main = do
args <- getArgs
case args of
[path] -> do
raw <- readFile' path
-- Part 1 explanation:
--
-- - Parse lines
-- - Compute 2D line intersections
let input = mapMaybe parseHailstone (lines raw)
xings = mapMaybe (uncurry intersect2) (pairs (fmap projectXY <$> input))
bounds = Rect2 (both 200000000000000) (both 400000000000000)
part1 = length (filter (inRect bounds) xings)
putStrLn $ "Part 1: " ++ show part1
-- Part 2 explanation:
--
-- - Let each hailstone i have position p[i] and velocity v[i], i.e. be a 3D line
-- - We're looking for a new line with position p and velocity v, such that for i = 1, ..., n and some t[i]:
--
-- p + t[i] * v = p[i] + t[i] * v[i] (note that t[i] is the same on both sides due to the time constraint)
-- <=> p - p[i] = t[i] * (v[i] - v) (i.e. (p - p[i]) and (v[i] - v) are collinear...)
-- <=> 0 = (p - p[i]) x (v[i] - v) (...therefore we can rewrite this with the cross product)
-- = (p - p[i]) x (v - v[i])
-- = (p x v) - (p x v[i]) - (p[i] x v) + (p[i] x v[i])
--
-- - Note that (p x v) is the only nonlinear term (because it contains two unknowns)
-- - The clever trick is now that since (p x v) is the same for every i,
-- we can just take two of these vector equations and subtract the
-- nonlinear term (p x v) to get rid of it (i >= 2 henceforth):
--
-- (p x v) - (p x v[1]) - (p[1] x v) + (p[1] x v[1])
-- - ((p x v) - (p x v[i]) - (p[i] x v) + (p[i] x v[i]))
-- = - (p x v[1]) - (p[1] x v) + (p[1] x v[1])
-- + (p x v[i]) + (p[i] x v) - (p[i] x v[i])
--
-- - Remembering that each vector equation expands to three scalar equations and that
-- we have 6 unknowns (px, py, pz, vx, vy, vz), we only have to pick two of these
-- remaining equations (we'll go with i = 2 and i = 3) to obtain 6 equations.
-- - These 6 equations form a linear system of equations, which we can solve via Gaussian elimination.
--
-- Credits for helping me crack this problem go to u/evouga: https://www.reddit.com/r/adventofcode/comments/18pnycy/comment/kepu26z
let p i = (input !! (i - 1)).pos
v i = (input !! (i - 1)).vel
[px,py,pz,vx,vy,vz] = map (round @_ @Integer) . fromRight' $ solveLinearSystem @Double LinearSystem
-- Unknowns: px py pz vx vy vz
{ a = [ [(v 2).y - (v 1).y, (v 1).x - (v 2).x, 0 , (p 1).y - (p 2).y, (p 2).x - (p 1).x, 0 ]
, [(v 3).y - (v 1).y, (v 1).x - (v 3).x, 0 , (p 1).y - (p 3).y, (p 3).x - (p 1).x, 0 ]
, [(v 1).z - (v 2).z, 0 , (v 2).x - (v 1).x, (p 2).z - (p 1).z, 0 , (p 1).x - (p 2).x]
, [(v 1).z - (v 3).z, 0 , (v 3).x - (v 1).x, (p 3).z - (p 1).z, 0 , (p 1).x - (p 3).x]
, [0 , (v 2).z - (v 1).z, (v 1).y - (v 2).y, 0 , (p 1).z - (p 2).z, (p 2).y - (p 1).y]
, [0 , (v 3).z - (v 1).z, (v 1).y - (v 3).y, 0 , (p 1).z - (p 3).z, (p 3).y - (p 1).y]
]
, b = [ (p 1).y * (v 1).x - (p 1).x * (v 1).y + (p 2).x * (v 2).y - (p 2).y * (v 2).x
, (p 1).y * (v 1).x - (p 1).x * (v 1).y + (p 3).x * (v 3).y - (p 3).y * (v 3).x
, (p 1).x * (v 1).z - (p 1).z * (v 1).x + (p 2).z * (v 2).x - (p 2).x * (v 2).z
, (p 1).x * (v 1).z - (p 1).z * (v 1).x + (p 3).z * (v 3).x - (p 3).x * (v 3).z
, (p 1).z * (v 1).y - (p 1).y * (v 1).z + (p 2).y * (v 2).z - (p 2).z * (v 2).y
, (p 1).z * (v 1).y - (p 1).y * (v 1).z + (p 3).y * (v 3).z - (p 3).z * (v 3).y
]
}
part2 = px + py + pz
putStrLn $ "Part 2: " ++ show part2 ++ " (p = " ++ show [px, py, pz] ++ ", v = " ++ show [vx, vy, vz] ++ ")"
_ -> do
putStrLn "Usage: day24 <path to input>"
exitFailure