C#

using QuickGraph;
using QuickGraph.Algorithms.ShortestPath;

namespace aoc24;

[ForDay(16)]
public class Day16 : Solver {
  private string[] data;
  private int width, height;
  private int start_x, start_y;
  private int end_x, end_y;

  private readonly (int, int)[] directions = [(1, 0), (0, 1), (-1, 0), (0, -1)];
  private record class Edge((int, int, int) Source, (int, int, int) Target) : IEdge<(int, int, int)>;

  private DelegateVertexAndEdgeListGraph<(int, int, int), Edge> graph;
  private AStarShortestPathAlgorithm<(int, int, int), Edge> search;

  private long min_distance;
  private List<(int, int, int)> min_distance_targets;

  public void Presolve(string input) {
    data = input.Trim().Split("\n");
    width = data[0].Length;
    height = data.Length;
    for (int i = 0; i < width; i++) {
      for (int j = 0; j < height; j++) {
        if (data[j][i] == 'S') {
          start_x = i;
          start_y = j;
        } else if (data[j][i] == 'E') {
          end_x = i;
          end_y = j;
        }
      }
    }
    graph = MakeGraph();
    var start = (start_x, start_y, 0);
    search = new AStarShortestPathAlgorithm<(int, int, int), Edge>(
      graph,
      edge => edge.Source.Item3 == edge.Target.Item3 ? 1 : 1000,
      vertex => Math.Abs(vertex.Item1 - start_x) + Math.Abs(vertex.Item2 - start_y) + 1000 *
          Math.Min(vertex.Item3, 4 - vertex.Item3)
      );
    Dictionary<(int, int, int), long> distances = [];
    search.SetRootVertex(start);
    search.ExamineVertex += vertex => {
      if (vertex.Item1 == end_x && vertex.Item2 == end_y) {
        distances[vertex] = (long)search.Distances[vertex];
      }
    };
    search.Compute();
    min_distance = distances.Values.Min();
    min_distance_targets = distances.Keys.Where(v => distances[v] == min_distance).ToList();
  }

  private DelegateVertexAndEdgeListGraph<(int, int, int), Edge> MakeGraph() => new(GetAllVertices(), GetOutEdges);

  private bool GetOutEdges((int, int, int) arg, out IEnumerable<Edge> result_enumerable) {
    List<Edge> result = [];
    var (x, y, dir) = arg;
    result.Add(new Edge(arg, (x, y, (dir + 1) % 4)));
    result.Add(new Edge(arg, (x, y, (dir + 3) % 4)));
    var (tx, ty) = (x + directions[dir].Item1, y + directions[dir].Item2);
    if (data[ty][tx] != '#') result.Add(new Edge(arg, (tx, ty, dir)));
    result_enumerable = result;
    return true;
  }

  private IEnumerable<(int, int, int)> GetAllVertices() {
    for (int i = 0; i < width; i++) {
      for (int j = 0; j < height; j++) {
        if (data[j][i] == '#') continue;
        yield return (i, j, 0);
        yield return (i, j, 1);
        yield return (i, j, 2);
        yield return (i, j, 3);
      }
    }
  }

  private HashSet<(int, int, int)> GetMinimumPathNodesTo((int, int, int) vertex) {
    var (x, y, dir) = vertex;
    if (x == start_x && y == start_y && dir == 0) return [vertex];
    if (!search.Distances.TryGetValue(vertex, out var distance_to_me)) return [];
    List<(int, int, int)> candidates = [
          (x, y, (dir + 1) % 4),
          (x, y, (dir + 3) % 4),
          (x - directions[dir].Item1, y - directions[dir].Item2, dir),
      ];
    HashSet<(int, int, int)> result = [vertex];
    foreach (var (cx, cy, cdir) in candidates) {
      if (!search.Distances.TryGetValue((cx, cy, cdir), out var distance_to_candidate)) continue;
      if (distance_to_candidate > distance_to_me - (dir == cdir ? 1 : 1000)) continue;
      result = result.Union(GetMinimumPathNodesTo((cx, cy, cdir))).ToHashSet();
    }
    return result;
  }

  public string SolveFirst() => min_distance.ToString();

  public string SolveSecond() => min_distance_targets
    .SelectMany(v => GetMinimumPathNodesTo(v))
    .Select(vertex => (vertex.Item1, vertex.Item2))
    .ToHashSet()
    .Count
    .ToString();
}

Uiua

Uiua’s new builtin path operator makes this a breeze. Given a function that returns valid neighbours for a point and their relative costs, and another function to test whether you have reached a valid goal, it gives the minimal cost, and all relevant paths. We just need to keep track of the current direction as we work through the maze.

Data  ← ≡°□°/$"_\n_" "###############\n#.......#....E#\n#.#.###.#.###.#\n#.....#.#...#.#\n#.###.#####.#.#\n#.#.#.......#.#\n#.#.#####.###.#\n#...........#.#\n###.#.#####.#.#\n#...#.....#.#.#\n#.#.#.###.#.#.#\n#.....#...#.#.#\n#.###.#.#.#.#.#\n#S..#.....#...#\n###############"
D₄    ← [1_0 ¯1_0 0_1 0_¯1]
End   ← ⊢⊚=@EData
Costs ← :∩▽⟜:≡(≠@#⊡:Data⊢).≡⊟⊙⟜(+1×1000¬≡/×=)+⟜:D₄∩¤°⊟
:path(Costs|≍End⊙◌°⊟)⊂:[1_0]⊢⊚=@SData
&p
&p ⧻◴≡⊢/◇⊂

Haskell

Rather busy today so late and somewhat messy! (Probably the same tomorrow…)

import Data.List
import Data.Map (Map)
import Data.Map qualified as Map
import Data.Maybe
import Data.Set (Set)
import Data.Set qualified as Set

readInput :: String -> Map (Int, Int) Char
readInput s = Map.fromList [((i, j), c) | (i, l) <- zip [0 ..] (lines s), (j, c) <- zip [0 ..] l]

bestPath :: Map (Int, Int) Char -> (Int, Set (Int, Int))
bestPath maze = go (Map.singleton start (0, Set.singleton startPos)) (Set.singleton start)
  where
    start = (startPos, (0, 1))
    walls = Map.keysSet $ Map.filter (== '#') maze
    [Just startPos, Just endPos] = map (\c -> fst <$> find ((== c) . snd) (Map.assocs maze)) ['S', 'E']
    go best edge
      | Set.null edge = Map.mapKeysWith mergePaths fst best Map.! endPos
      | otherwise =
          let nodes' =
                filter (\(x, (c, _)) -> maybe True ((c <=) . fst) $ best Map.!? x) $
                  concatMap (step . (\x -> (x, best Map.! x))) (Set.elems edge)
              best' = foldl' (flip $ uncurry $ Map.insertWith mergePaths) best nodes'
           in go best' $ Set.fromList (map fst nodes')
    step ((p@(i, j), d@(di, dj)), (cost, path)) =
      let rots = [((p, d'), (cost + 1000, path)) | d' <- [(-dj, di), (dj, -di)]]
          moves =
            [ ((p', d), (cost + 1, Set.insert p' path))
              | let p' = (i + di, j + dj),
                p `Set.notMember` walls
            ]
       in moves ++ rots
    mergePaths a@(c1, p1) b@(c2, p2) =
      case compare c1 c2 of
        LT -> a
        GT -> b
        EQ -> (c1, Set.union p1 p2)

main = do
  (score, visited) <- bestPath . readInput <$> readFile "input16"
  print score
  print (Set.size visited)

Python

Part 1: Run Dijkstra’s algorithm to find shortest path.

I chose to represent nodes using the location (i, j) as well as the direction dir faced by the reindeer.
Initially I tried creating the complete adjacency graph but that lead to max recursion so I ended up populating graph for only the nodes I was currently exploring.

Part 2: Track paths while performing Dijkstra’s algorithm.

First, I modified the algorithm to look through neighbors with equal cost along with the ones with lesser cost, so that it would go through all shortest paths.
Then, I keep track of the list of previous nodes for every node explored.
Finally, I use those lists to run through the paths backwards, taking note of all unique locations.

Code:

import os

# paths
here = os.path.dirname(os.path.abspath(__file__))
filepath = os.path.join(here, "input.txt")

# read input
with open(filepath, mode="r", encoding="utf8") as f:
    data = f.read()

from collections import defaultdict
from dataclasses import dataclass
import heapq as hq
import math

# up, right, down left
DIRECTIONS = [(-1, 0), (0, 1), (1, 0), (0, -1)]


# Represent a node using its location and the direction
@dataclass(frozen=True)
class Node:
    i: int
    j: int
    dir: int


maze = data.splitlines()
m, n = len(maze), len(maze[0])

# we always start from bottom-left corner (facing east)
start_node = Node(m - 2, 1, 1)
# we always end in top-right corner (direction doesn't matter)
end_node = Node(1, n - 2, -1)

# the graph will be updated lazily because it is too much processing
#   to completely populate it beforehand
graph = defaultdict(list)
# track nodes whose all edges have been explored
visited = set()
# heap to choose next node to explore
# need to add id as middle tuple element so that nodes dont get compared
min_heap = [(0, id(start_node), start_node)]
# min distance from start_node to node so far
# missing values are treated as math.inf
min_dist = {}
min_dist[start_node] = 0
# keep track of all previous nodes for making path
prev_nodes = defaultdict(list)


# utility method for debugging (prints the map)
def print_map(current_node, prev_nodes):
    pns = set((n.i, n.j) for n in prev_nodes)
    for i in range(m):
        for j in range(n):
            if i == current_node.i and j == current_node.j:
                print("X", end="")
            elif (i, j) in pns:
                print("O", end="")
            else:
                print(maze[i][j], end="")
        print()


# Run Dijkstra's algo
while min_heap:
    cost_to_node, _, node = hq.heappop(min_heap)
    if node in visited:
        continue
    visited.add(node)

    # early exit in the case we have explored all paths to the finish
    if node.i == end_node.i and node.j == end_node.j:
        # assign end so that we know which direction end was reached by
        end_node = node
        break

    # update adjacency graph from current node
    di, dj = DIRECTIONS[node.dir]
    if maze[node.i + di][node.j + dj] != "#":
        moved_node = Node(node.i + di, node.j + dj, node.dir)
        graph[node].append((moved_node, 1))
    for x in range(3):
        rotated_node = Node(node.i, node.j, (node.dir + x + 1) % 4)
        graph[node].append((rotated_node, 1000))

    # explore edges
    for neighbor, cost in graph[node]:
        cost_to_neighbor = cost_to_node + cost
        # The following condition was changed from > to >= because we also want to explore
        #   paths with the same cost, not just better cost
        if min_dist.get(neighbor, math.inf) >= cost_to_neighbor:
            min_dist[neighbor] = cost_to_neighbor
            prev_nodes[neighbor].append(node)
            # need to add id as middle tuple element so that nodes dont get compared
            hq.heappush(min_heap, (cost_to_neighbor, id(neighbor), neighbor))

print(f"Part 1: {min_dist[end_node]}")

# PART II: Run through the path backwards, making note of all coords

visited = set([start_node])
path_locs = set([(start_node.i, start_node.j)])  # all unique locations in path
stack = [end_node]

while stack:
    node = stack.pop()
    if node in visited:
        continue
    visited.add(node)

    path_locs.add((node.i, node.j))

    for prev_node in prev_nodes[node]:
        stack.append(prev_node)

print(f"Part 2: {len(path_locs)}")