Day 24 part 2

1 files changed, 150 insertions, 6 deletions

diff --git a/2023/src/bin/day_24.rs b/2023/src/bin/day_24.rs
index 325b4cc..1f22169 100644
--- a/2023/src/bin/day_24.rs
+++ b/2023/src/bin/day_24.rs
@@ -1,4 +1,4 @@
-use nalgebra::{Matrix2, Point2, Point3, Vector2, Vector3};
+use nalgebra::{Matrix2, Matrix6, Point3, RowVector6, Vector2, Vector3, Vector6};
 use nom::{
     bytes::complete::tag,
     character::complete::{i64, line_ending},
@@ -14,6 +14,12 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
     let parsed = Hailstones::parser(&input).unwrap().1;
     dbg!(&parsed.count_intersections_2d(200000000000000., 400000000000000.));
 
+    let magic_rock = parsed.find_magic_throwing_rock();
+    dbg!(&magic_rock);
+    dbg!(
+        magic_rock.position.x as i64 + magic_rock.position.y as i64 + magic_rock.position.z as i64
+    );
+
     Ok(())
 }
 
@@ -22,7 +28,7 @@ struct Hailstones(Vec<Hailstone>);
 
 #[derive(Debug)]
 struct Hailstone {
-    initial_position: Point3<f64>,
+    position: Point3<f64>,
     velocity: Vector3<f64>,
 }
 
@@ -44,6 +50,134 @@ impl Hailstones {
             })
             .sum()
     }
+
+    fn find_magic_throwing_rock(&self) -> Hailstone {
+        (0..self.0.len())
+            .flat_map(move |h1| {
+                (h1 + 1..self.0.len()).flat_map(move |h2| {
+                    (h2 + 1..self.0.len())
+                        .flat_map(move |h3| (h3 + 1..self.0.len()).map(move |h4| [h1, h2, h3, h4]))
+                })
+            })
+            .take(1000000)
+            .map(|hailstones| {
+                let rock = self.find_magic_throwing_rock_for_hailstones(hailstones);
+                (
+                    // the solution I'm after uses integers. This tries to find the minimum error.
+                    rock.position.x.abs().fract()
+                        + rock.position.y.abs().fract()
+                        + rock.position.z.abs().fract()
+                        + rock.velocity.x.abs().fract()
+                        + rock.velocity.y.abs().fract()
+                        + rock.velocity.z.abs().fract()
+                        + self
+                            .0
+                            .iter()
+                            .map(|h| rock.collition_time(&h).fract())
+                            .sum::<f64>(),
+                    rock,
+                )
+            })
+            .min_by(|(error_a, _), (error_b, _)| error_a.total_cmp(error_b))
+            .unwrap()
+            .1
+    }
+
+    fn find_magic_throwing_rock_for_hailstones(&self, hailstones: [usize; 4]) -> Hailstone {
+        // unknowns are (x, y, z, dx, dy, dz)
+        let h1 = &self.0[hailstones[0]];
+        let h2 = &self.0[hailstones[1]];
+        let h3 = &self.0[hailstones[2]];
+        let h4 = &self.0[hailstones[3]];
+
+        let coefficients: Matrix6<f64> = Matrix6::from_rows(&[
+            RowVector6::new(
+                h2.velocity.y - h1.velocity.y,
+                h1.velocity.x - h2.velocity.x,
+                0.,
+                h1.position.y - h2.position.y,
+                h2.position.x - h1.position.x,
+                0.,
+            ),
+            RowVector6::new(
+                h2.velocity.z - h1.velocity.z,
+                0.,
+                h1.velocity.x - h2.velocity.x,
+                h1.position.z - h2.position.z,
+                0.,
+                h2.position.x - h1.position.x,
+            ),
+            RowVector6::new(
+                0.,
+                h2.velocity.z - h1.velocity.z,
+                h1.velocity.y - h2.velocity.y,
+                0.,
+                h1.position.z - h2.position.z,
+                h2.position.y - h1.position.y,
+            ),
+            RowVector6::new(
+                h4.velocity.y - h3.velocity.y,
+                h3.velocity.x - h4.velocity.x,
+                0.,
+                h3.position.y - h4.position.y,
+                h4.position.x - h3.position.x,
+                0.,
+            ),
+            RowVector6::new(
+                h4.velocity.z - h3.velocity.z,
+                0.,
+                h3.velocity.x - h4.velocity.x,
+                h3.position.z - h4.position.z,
+                0.,
+                h4.position.x - h3.position.x,
+            ),
+            RowVector6::new(
+                0.,
+                h4.velocity.z - h3.velocity.z,
+                h3.velocity.y - h4.velocity.y,
+                0.,
+                h3.position.z - h4.position.z,
+                h4.position.y - h3.position.y,
+            ),
+        ]);
+        let constants: Vector6<f64> = Vector6::new(
+            h1.position.y * h1.velocity.x
+                - h1.position.x * h1.velocity.y
+                - h2.position.y * h2.velocity.x
+                + h2.position.x * h2.velocity.y,
+            h1.position.z * h1.velocity.x
+                - h1.position.x * h1.velocity.z
+                - h2.position.z * h2.velocity.x
+                + h2.position.x * h2.velocity.z,
+            h1.position.z * h1.velocity.y
+                - h1.position.y * h1.velocity.z
+                - h2.position.z * h2.velocity.y
+                + h2.position.y * h2.velocity.z,
+            h3.position.y * h3.velocity.x
+                - h3.position.x * h3.velocity.y
+                - h4.position.y * h4.velocity.x
+                + h4.position.x * h4.velocity.y,
+            h3.position.z * h3.velocity.x
+                - h3.position.x * h3.velocity.z
+                - h4.position.z * h4.velocity.x
+                + h4.position.x * h4.velocity.z,
+            h3.position.z * h3.velocity.y
+                - h3.position.y * h3.velocity.z
+                - h4.position.z * h4.velocity.y
+                + h4.position.y * h4.velocity.z,
+        );
+
+        if let Some(coefficients_inverse) = coefficients.try_inverse() {
+            let unknowns = coefficients_inverse * constants;
+
+            Hailstone {
+                position: Point3::new(unknowns[0], unknowns[1], unknowns[2]),
+                velocity: Vector3::new(unknowns[3], unknowns[4], unknowns[5]),
+            }
+        } else {
+            panic!("No solution found, matrix didn't invert")
+        }
+    }
 }
 
 impl Hailstone {
@@ -61,7 +195,7 @@ impl Hailstone {
                 ),
             ),
             |(initial_position, velocity)| Hailstone {
-                initial_position,
+                position: initial_position,
                 velocity,
             },
         )(input)
@@ -75,8 +209,8 @@ impl Hailstone {
             -other.velocity.y,
         );
         let constants = Vector2::new(
-            other.initial_position.x - self.initial_position.x,
-            other.initial_position.y - self.initial_position.y,
+            other.position.x - self.position.x,
+            other.position.y - self.position.y,
         );
 
         if let Some(variables_inverse) = variables.try_inverse() {
@@ -84,7 +218,7 @@ impl Hailstone {
             let self_t = intersection.x;
             let other_t = intersection.y;
 
-            let intersection = self.initial_position.xy() + self.velocity.xy() * self_t;
+            let intersection = self.position.xy() + self.velocity.xy() * self_t;
             self_t >= 0.
                 && other_t >= 0.
                 && intersection.x >= min
@@ -95,4 +229,14 @@ impl Hailstone {
             false
         }
     }
+
+    /// This is only intended for hail that definitely collides!
+    fn collition_time(&self, other: &Hailstone) -> f64 {
+        let tx = (self.position.x - other.position.x) / (other.velocity.x - self.velocity.x);
+        let ty = (self.position.y - other.position.y) / (other.velocity.y - self.velocity.y);
+        let tz = (self.position.z - other.position.z) / (other.velocity.z - self.velocity.z);
+        let t = tx.max(ty).max(tz); // sometimes one of these is zero!
+        assert!(t > 0.);
+        t
+    }
 }