2023 day 24: fix some rounding issues

yut23 · Jan 31, 2024 · 3d6f8a8 · 3d6f8a8
1 parent 1d05458
commit 3d6f8a8
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 17 deletions.
diff --git a/2023/src/day24.cpp b/2023/src/day24.cpp
@@ -9,8 +9,7 @@
 #include "gauss_elim.hpp" // for gauss_elim, solve_upper_triangular, RowPermuter
 #include "lib.hpp"        // for parse_args, DEBUG
 #include "unit_test/pretty_print.hpp" // for repr
-#include <cmath>                      // for round
-#include <cstdint>                    // for int64_t
+#include <cmath>                      // for llround
 #include <iostream>                   // for cout, cerr
 #include <vector>                     // for vector
 // IWYU pragma: no_include <algorithm>  // for fill_n
@@ -67,13 +66,17 @@ std::int64_t part_2(const std::vector<aoc::day24::Hailstone> &stones) {
                                   << "b: " << pretty_print::repr(b) << "\n";
                     }
                     if (rp) {
+                        if constexpr (aoc::DEBUG) {
+                            std::cerr << "row order: "
+                                      << pretty_print::repr(rp->row_order)
+                                      << "\n";
+                        }
                         auto x = solve_upper_triangular(A, b, *rp);
                         if constexpr (aoc::DEBUG) {
                             std::cerr << "x: " << pretty_print::repr(x) << "\n";
                         }
-                        return std::round<std::int64_t>(x[0]) +
-                               std::round<std::int64_t>(x[1]) +
-                               std::round<std::int64_t>(x[2]);
+                        return std::llround(x[0]) + std::llround(x[1]) +
+                               std::llround(x[2]);
                     }
                 }
             }

diff --git a/2023/src/day24.hpp b/2023/src/day24.hpp
@@ -23,8 +23,8 @@
 namespace aoc::day24 {
 
 struct Hailstone {
-    std::int64_t px, py, pz;
-    std::int64_t vx, vy, vz;
+    long double px, py, pz;
+    long double vx, vy, vz;
 };
 
 std::istream &operator>>(std::istream &is, Hailstone &stone) {
@@ -64,15 +64,11 @@ find_intersection_xy(const Hailstone &a, const Hailstone &b) {
      *  Simplify@First@Solve[y-#<>"py"==(#<>"vy"/#<>"vx")(x-#<>"px")&/@{"a.","b."},{x,y}],
      * "\n"]
      */
-    // cast the large position values to doubles, so we don't get integer
-    // overflow
-    double apx = a.px, apy = a.py;
-    double bpx = b.px, bpy = b.py;
-    double x = (apy * a.vx * b.vx - apx * a.vy * b.vx - a.vx * bpy * b.vx +
-                a.vx * bpx * b.vy) /
+    double x = (a.py * a.vx * b.vx - a.px * a.vy * b.vx - a.vx * b.py * b.vx +
+                a.vx * b.px * b.vy) /
                (-(a.vy * b.vx) + a.vx * b.vy);
-    double y = (a.vy * bpy * b.vx - apy * a.vx * b.vy + apx * a.vy * b.vy -
-                a.vy * bpx * b.vy) /
+    double y = (a.vy * b.py * b.vx - a.py * a.vx * b.vy + a.px * a.vy * b.vy -
+                a.vy * b.px * b.vy) /
                (a.vy * b.vx - a.vx * b.vy);
 
     if (std::copysign(1.0, x - a.px) != std::copysign(1.0, a.vx) ||
@@ -85,7 +81,7 @@ find_intersection_xy(const Hailstone &a, const Hailstone &b) {
 }
 
 template <class T>
-std::pair<aoc::ds::Grid<T>, std::vector<T>>
+std::tuple<aoc::ds::Grid<T>, std::vector<T>>
 make_system(const std::vector<aoc::day24::Hailstone> &stones,
             const std::pair<std::size_t, std::size_t> &pair_1,
             const std::pair<std::size_t, std::size_t> &pair_2) {

diff --git a/2023/src/day24.py b/2023/src/day24.py
@@ -159,7 +159,7 @@ def make_linear_system(
 
 
 def part_2_linalg(stones: list[Hailstone]) -> int:
-    index_pairs = ((0, 1), (2, 3))
+    index_pairs = ((0, 1), (0, 2))
     A, b = make_linear_system(stones, index_pairs)
     # pm1, pm2, pm3 = stones[index_pairs[0][0]].pos
     # vm1, vm2, vm3 = stones[index_pairs[0][0]].vel