diff --git a/src/day21.zig b/src/day21.zig
new file mode 100644
index 0000000..d73662a
--- /dev/null
+++ b/src/day21.zig
@@ -0,0 +1,254 @@
+const aoc = @import("./aoc.zig");
+const std = @import("std");
+const Allocator = std.mem.Allocator;
+const assert = std.debug.assert;
+const PointU32 = @import("./point.zig").PointU32;
+const PointI32 = @import("./point.zig").PointI32;
+
+const Tile = enum {
+    GardenPlot,
+    Rock,
+    Start,
+
+    fn fromU8(char: u8) ?Tile {
+        return switch (char) {
+            '.' => .GardenPlot,
+            '#' => .Rock,
+            'S' => .Start,
+            else => null
+        };
+    }
+};
+
+const Map = struct {
+    allocator: Allocator,
+    tiles: []Tile,
+    width: usize,
+    height: usize,
+
+    fn deinit(self: Map) void {
+        self.allocator.free(self.tiles);
+    }
+};
+
+fn parseInput(allocator: Allocator, lines: []const []const u8) !Map {
+    const width = lines[0].len;
+    const height = lines.len;
+    const parsed = Map {
+        .allocator = allocator,
+        .width = width,
+        .height = height,
+        .tiles = try allocator.alloc(Tile, width * height)
+    };
+    errdefer parsed.deinit();
+
+    for (0.., lines) |y, line| {
+        for (0.., line) |x, char| {
+            const index = y * width + x;
+            parsed.tiles[index] = Tile.fromU8(char) orelse return error.InvalidTile;
+        }
+    }
+
+    return parsed;
+}
+
+fn findStart(map: Map) ?PointU32 {
+    for (0..map.height) |y| {
+        for (0..map.width) |x| {
+            if (map.tiles[y * map.width + x] == .Start) {
+                return PointU32{
+                    .x = @intCast(x),
+                    .y = @intCast(y),
+                };
+            }
+        }
+    }
+
+    return null;
+}
+
+fn simulateSteps(allocator: Allocator, map: Map, from: PointU32, num_steps: usize) !u64 {
+    var answer: usize = 0;
+
+    var seen = std.AutoHashMap(PointU32, void).init(allocator);
+    defer seen.deinit();
+
+    var queue = std.ArrayList(struct { point: PointU32, steps_left: usize }).init(allocator);
+    defer queue.deinit();
+
+    try seen.put(from, {});
+    try queue.append(.{ .point = from, .steps_left = num_steps });
+
+    while (queue.items.len > 0) {
+        const queue_item = queue.orderedRemove(0);
+
+        if (queue_item.steps_left % 2 == 0) {
+            answer += 1;
+        }
+        if (queue_item.steps_left == 0) {
+            continue;
+        }
+
+        const point = queue_item.point;
+
+        const offsets = .{
+            .{  1,  0 },
+            .{  0,  1 },
+            .{ -1,  0 },
+            .{  0, -1 },
+        };
+        inline for (offsets) |offset| {
+            const nx: i32 = @as(i32, @intCast(point.x)) + offset[0];
+            const ny: i32 = @as(i32, @intCast(point.y)) + offset[1];
+            if ((0 <= nx and nx < map.width) and (0 <= ny and ny < map.height)) {
+                const next_point = PointU32{
+                    .x = @intCast(nx),
+                    .y = @intCast(ny),
+                };
+                if (map.tiles[next_point.x * map.width + next_point.y] != .Rock and !seen.contains(next_point)) {
+                    try queue.append(.{ .point = next_point, .steps_left = queue_item.steps_left - 1 });
+                    try seen.put(next_point, {});
+                }
+            }
+        }
+    }
+
+    return answer;
+}
+
+pub fn part1(input: *aoc.Input) !aoc.Result {
+    const allocator = input.allocator;
+    const map = try parseInput(allocator, input.lines);
+    defer map.deinit();
+
+    const start_location = findStart(map) orelse return error.NoStart;
+    const answer = try simulateSteps(allocator, map, start_location, 64);
+
+    return .{ .uint = answer };
+}
+
+fn isColumnEmpty(map: Map, x: usize) bool {
+    for (0..map.height) |y| {
+        const idx = y * map.width + x;
+        if (map.tiles[idx] == .Rock) {
+            return false;
+        }
+    }
+    return true;
+}
+
+fn isRowEmpty(map: Map, y: usize) bool {
+    for (0..map.width) |x| {
+        const idx = y * map.width + x;
+        if (map.tiles[idx] == .Rock) {
+            return false;
+        }
+    }
+    return true;
+}
+
+// Props to this https://github.com/villuna/aoc23/wiki/A-Geometric-solution-to-advent-of-code-2023,-day-21
+// For writting up a nice solution explanation
+// Also this was very useful: https://www.youtube.com/watch?v=9UOMZSL0JTg
+pub fn part2(input: *aoc.Input) !aoc.Result {
+    const allocator = input.allocator;
+    const map = try parseInput(allocator, input.lines);
+    defer map.deinit();
+
+    const start_location = findStart(map) orelse return error.NoStart;
+
+    // "Some" assumptions about the data must be made
+    assert(map.width == map.height); // The grid is a square
+    assert(start_location.x == map.width / 2 and start_location.y == map.height / 2); // Starting location is in the center
+    assert(map.width % 2 == 1); // The size of the square grid is odd
+    assert(isRowEmpty(map, start_location.y) and isColumnEmpty(map, start_location.x));
+    assert(isRowEmpty(map, 0) and isRowEmpty(map, map.height-1));
+    assert(isColumnEmpty(map, 0) and isColumnEmpty(map, map.width-1));
+
+    const map_size: u32 = @intCast(map.width);
+
+    const num_steps = 26501365;
+
+    const n = (num_steps - map_size / 2) / map_size;
+
+    var answer: usize = 0;
+
+    const odd_points = try simulateSteps(allocator, map, start_location, 2*map_size + 1);
+    const odd_squares = std.math.pow(usize, (n - 1) / 2 * 2 + 1, 2);
+    answer += odd_points * odd_squares;
+
+    const even_points = try simulateSteps(allocator, map, start_location, 2*map_size);
+    const even_squares = std.math.pow(usize, n / 2 * 2, 2);
+    answer += even_points * even_squares;
+
+    const cardinal_corners = .{
+        PointU32{ .x = start_location.x, .y = map_size - 1 },
+        PointU32{ .x = start_location.x, .y = 0 },
+        PointU32{ .x = map_size - 1, .y = start_location.y },
+        PointU32{ .x = 0, .y = start_location.y }
+    };
+    inline for(cardinal_corners) |corner| {
+        answer += try simulateSteps(allocator, map, corner, map_size - 1);
+    }
+
+    const small_corners = .{
+        PointU32{ .x = 0            , .y = map_size - 1 } ,
+        PointU32{ .x = map_size - 1 , .y = map_size - 1 } ,
+        PointU32{ .x = 0            , .y = 0 } ,
+        PointU32{ .x = map_size - 1 , .y = 0 } ,
+    };
+    inline for(small_corners) |corner| {
+        const points = try simulateSteps(allocator, map, corner, map_size / 2 - 1);
+        answer += points * n;
+    }
+
+    const large_corners = .{
+        PointU32{ .x = 0            , .y = map_size - 1 } ,
+        PointU32{ .x = map_size - 1 , .y = map_size - 1 } ,
+        PointU32{ .x = 0            , .y = 0 } ,
+        PointU32{ .x = map_size - 1 , .y = 0 } ,
+    };
+    inline for(large_corners) |corner| {
+        const points = try simulateSteps(allocator, map, corner, map_size * 3 / 2 - 1);
+        answer += points * (n - 1);
+    }
+
+    return .{ .uint = answer };
+}
+
+
+test "part 1 example" {
+    const example_input = [_][]const u8{
+        "...........",
+        ".....###.#.",
+        ".###.##..#.",
+        "..#.#...#..",
+        "....#.#....",
+        ".##..S####.",
+        ".##..#...#.",
+        ".......##..",
+        ".##.#.####.",
+        ".##..##.##.",
+        "...........",
+    };
+
+    try aoc.expectAnswerUInt(part1, 42, &example_input);
+}
+
+test "part 2 example" {
+    const modified_example_input = [_][]const u8{
+        "...........",
+        "......##.#.",
+        ".###..#..#.",
+        "..#.#...#..",
+        "....#.#....",
+        ".....S.....",
+        ".##......#.",
+        ".......##..",
+        ".##.#.####.",
+        ".##...#.##.",
+        "...........",
+    };
+
+    try aoc.expectAnswerUInt(part2, 528192639411648, &modified_example_input);
+}