--- /dev/null
+-- Writeup at https://work.njae.me.uk/2023/12/29/advent-of-code-2023-day-21/
+
+import AoC
+import Linear (V2(..), (^+^))
+import qualified Data.Set as S
+import Data.Ix
+import Data.List
+import Data.List.Split
+
+type Position = V2 Int -- r, c
+type Grid = S.Set Position
+
+main :: IO ()
+main =
+ do dataFileName <- getDataFileName
+ text <- readFile dataFileName
+ let (rocks, bounds) = mkGrid text
+ let start = findStart text
+ -- print rocks
+ -- print bounds
+ -- print start
+ print $ part1 rocks bounds start
+ -- print $ part2Quadratic rocks bounds start
+ print $ part2 rocks bounds start -- ((26501365 - 65) `div` (131 * 2))
+
+part1, part2, part2Quadratic :: Grid -> (Position, Position) -> Grid -> Int
+part1 rocks bounds start = S.size $ (!! 64) $ iterate (takeSteps rocks bounds) start
+
+part2Quadratic rocks bounds start = a * x ^ 2 + b * x + c -- (p, q, r, a * x ^ 2 + b * x + c)
+ where ps = (!! 65) $ iterate (takeSteps rocks bounds) start
+ qs = (!! 131) $ iterate (takeSteps rocks bounds) ps
+ rs = (!! 131) $ iterate (takeSteps rocks bounds) qs
+ p = S.size ps
+ q = S.size qs
+ r = S.size rs
+ c = p
+ b = (4*q - 3*p - r) `div` 2
+ a = q - p - b
+ s = 26501365
+ x = (s - 65) `div` 131
+
+part2 rocks bounds start =
+ nEvenFilled + nOddFilled + nEdges + upperPoint + lowerPoint + rCap + lCap
+ where (V2 minR _, V2 maxR _) = bounds
+ tileWidth = maxR - minR + 1
+ (doubleTileSteps, extraSteps) = divMod 26501365 (tileWidth * 2)
+ positions = (!! (tileWidth * 2 + extraSteps)) $ iterate (takeSteps rocks bounds) start
+
+ evenFilled = countInRange bounds (V2 0 0) positions
+ oddFilled = countInRange bounds (V2 tileWidth 0) positions
+ upperPoint = countInRange bounds (V2 (tileWidth * -2) 0) positions
+ lowerPoint = countInRange bounds (V2 (tileWidth * 2) 0) positions
+ urEdge = countInRange bounds (V2 (- tileWidth) tileWidth) positions +
+ countInRange bounds (V2 (- tileWidth * 2) tileWidth) positions
+ lrEdge = countInRange bounds (V2 tileWidth tileWidth) positions +
+ countInRange bounds (V2 ( tileWidth * 2) tileWidth) positions
+ ulEdge = countInRange bounds (V2 (- tileWidth) (- tileWidth)) positions +
+ countInRange bounds (V2 (- tileWidth * 2) (- tileWidth)) positions
+ llEdge = countInRange bounds (V2 tileWidth (- tileWidth)) positions +
+ countInRange bounds (V2 ( tileWidth * 2) (- tileWidth)) positions
+ rCap = countInRange bounds (V2 (- tileWidth) ( tileWidth * 2)) positions +
+ countInRange bounds (V2 0 ( tileWidth * 2)) positions +
+ countInRange bounds (V2 tileWidth ( tileWidth * 2)) positions
+ lCap = countInRange bounds (V2 (- tileWidth) (- tileWidth * 2)) positions +
+ countInRange bounds (V2 0 (- tileWidth * 2)) positions +
+ countInRange bounds (V2 tileWidth (- tileWidth * 2)) positions
+ edges = urEdge + lrEdge + ulEdge + llEdge
+
+ nEvenFilled = (doubleTileSteps * 2 - 1) ^ 2 * evenFilled
+ nOddFilled = (doubleTileSteps * 2) ^ 2 * oddFilled
+ nEdges = (doubleTileSteps * 2 - 1) * edges
+
+
+countInRange :: (Position, Position) -> Position -> Grid -> Int
+countInRange bounds delta cells =
+ S.size $ S.filter (inRange (shiftBounds bounds delta)) cells
+
+shiftBounds :: (Position, Position) -> Position -> (Position, Position)
+shiftBounds (a, b) d = (a ^+^ d, b ^+^ d)
+
+takeSteps :: Grid -> (Position, Position) -> Grid -> Grid
+takeSteps rocks bounds places =
+ S.filter (notAtRock rocks bounds) $ S.unions $ S.map adjacents places
+
+notAtRock :: Grid -> (Position, Position) -> Position -> Bool
+notAtRock rocks (_, V2 maxR maxC) (V2 r c) = here' `S.notMember` rocks
+ where here' = V2 (r `mod` (maxR + 1)) (c `mod` (maxC + 1))
+
+adjacents :: Position -> Grid
+adjacents here = S.map (here ^+^) $ S.fromList [ V2 0 1, V2 1 0, V2 0 (-1), V2 (-1) 0 ]
+
+
+showGrid :: Grid -> (Position, Position) -> Grid -> String
+showGrid rocks bounds cells = unlines $ intercalate [" "] $ chunksOf 7 $ fmap (showRow rocks bounds cells) [-28..34]
+ where showRow rocks bounds cells r = intercalate " " $ chunksOf 7 $ fmap ((showCell rocks bounds cells) . V2 r) [-28..34]
+ showCell rocks bounds cells here
+ | not $ notAtRock rocks bounds here = '#'
+ | here `S.member` cells = 'O'
+ | otherwise = '.'
+-- showCell :: Grid -> (Position, Position) -> Grid -> Position -> Char
+
+-- reading the map
+
+mkGrid :: String -> (Grid, (Position, Position))
+mkGrid text = ( S.fromList [ V2 r c | r <- [0..maxR], c <- [0..maxC]
+ , rows !! r !! c == '#'
+ ]
+ , (V2 0 0, V2 maxR maxC)
+ )
+ where rows = lines text
+ maxR = length rows - 1
+ maxC = (length $ head rows) - 1
+
+findStart :: String -> Grid
+findStart text = S.fromList [ V2 r c | r <- [0..maxR], c <- [0..maxC]
+ , rows !! r !! c == 'S'
+ ]
+ where rows = lines text
+ maxR = length rows - 1
+ maxC = (length $ head rows) - 1
+