aoc 2023 day 24 part 2 solver (does not work)

aoc-2023-day24-part2-solver.cs

enum Plane { X = 0, Y = 1, Z = 2 }
struct DimensionTrajectory
{
	public DimensionTrajectory(Hailstone h, Plane p)
	{
		Position = p switch { Plane.X => h.pos.X, Plane.Y => h.pos.Y, Plane.Z => h.pos.Z };
		Velocity = p switch { Plane.X => h.vel.Dx, Plane.Y => h.vel.Dy, Plane.Z => h.vel.Dz };
	}
	public readonly long Position;
	public readonly long Velocity;
}

object TrySolvePart2(List<Hailstone> hailstones)
{

	// okay, this gives us four unknowns.  I can add two more unknowns (z0 and vz0) and have their coefficients be zero.
	/*
Let the rock be located at $ (x_0, y_0, z_0) $ with velocity $ ({v_x}_0, {v_y}_0, {v_z}_0) $

Let hailstone $n, 1 <= n <= 300$ be located at $ (x_n, y_n, z_n) $ with velocity $ ({v_x}_n, {v_y}_n, {v_z}_n ) $

$$
\begin{align*}
x_0 + {v_x}_0 t_1 = x_1 + {v_x}_1 t_1 \\
y_0 + {v_y}_0 t_1 = y_1 + {v_y}_1 t_1 \\
\\
x_0 + {v_x}_0 t_2 = x_2 + {v_x}_2 t_2 \\
y_0 + {v_y}_0 t_2 = y_2 + {v_y}_2 t_2 \\
\end{align*}
$$
Solving for $ t_n $ :

$$
\begin{align*}
t_1 = (x_1 - x_0) / ({v_x}_0 - {v_x}_1) \\
t_1 = (y_1 - y_0) / ({v_y}_0 - {v_y}_1) \\

t_2 = (x_2 - x_0) / ({v_x}_0 - {v_x}_2) \\
t_2 = (y_2 - y_0) / ({v_y}_0 - {v_y}_2) \\
\end{align*}
$$

thus

$$
\begin{align*}
(x_1 - x_0) / ({v_x}_0 - {v_x}_1) = (y_1 - y_0) / ({v_y}_0 - {v_y}_1) \\
(x_2 - x_0) / ({v_x}_0 - {v_x}_2) = (y_2 - y_0) / ({v_y}_0 - {v_y}_2)
\end{align*}
$$

rearranging:

$$
\begin{align*}
x_0 {v_y}_0 - y_0 {v_x}_0 = x_0 {v_y}_1 - y_0 {v_x}_1 + y_1 {v_x}_1 - x_1 {v_y}_1 + x_1 {v_y}_0 - y_1 {v_x}_0 \\
x_0 {v_y}_0 - y_0 {v_x}_0 = x_0 {v_y}_2 - y_0 {v_x}_2 + y_2 {v_x}_2 - x_2 {v_y}_2 + x_2 {v_y}_0 - y_2 {v_x}_0 
\end{align*}
$$
	
*/
	/* 
	0 = vy1 x0 - vx1 y0 + y1 vx1 - x1 vy1 + x1 vy0 - y1 vx0 -
       (vy2 x0 - vx2 y0 + y2 vx2 - x2 vy2 + x2 vy0 - y2 vx0)

   0 = (vy2 - vy1) x0 - (vx2 - vx1) y0 + (y2 vx2 - y1 vx1) - (x2 vy2 - x1 vy1) + (x2 - x1) vy0 - (y2 - y1) vx0

   (vy1 - vy2) x0 + (vx2 - vx1) y0 + (y2 - y1) vx0 + (x1 - x2) vy0 = (y2 vx2 - y1 vx1) - (x2 vy2 - x1 vy1)

	*/

	double[,] matrix = new double[6, 6];
	double[] answers = new double[6];
	void AddEquation(int index, int p1idx, int p2idx, DimensionTrajectory h1x, DimensionTrajectory h1y, DimensionTrajectory h2x, DimensionTrajectory h2y)
	{
		answers[index] = (h2y.Position * h2x.Velocity - h2x.Position - h2y.Velocity) - (h1y.Position * h1x.Velocity - h1x.Position * h1y.Velocity);
		matrix[index, p1idx] = h1y.Velocity - h2y.Velocity;
		matrix[index, p2idx] = h2x.Velocity - h1x.Velocity;
		matrix[index, 3 + p1idx] = h2y.Position - h1y.Position;
		matrix[index, 3 + p2idx] = h1x.Position - h2x.Position;
	}
	void AddEquationByPlane(int index, Plane p1, Plane p2, Hailstone h1, Hailstone h2)
	{
		AddEquation(index, (int)p1, (int)p2, new DimensionTrajectory(h1, p1), new DimensionTrajectory(h1, p2), new DimensionTrajectory(h2, p1), new DimensionTrajectory(h2, p2));
	}
	int AddEquations(int index, Hailstone h1, Hailstone h2)
	{
		AddEquationByPlane(index, Plane.X, Plane.Y, h1, h2);
		AddEquationByPlane(index + 1, Plane.X, Plane.Z, h1, h2);
		AddEquationByPlane(index + 2, Plane.Y, Plane.Z, h1, h2);
		return 3;
	}
	int offset = 0;
	offset += AddEquations(offset, hailstones[0], hailstones[1]);
	offset += AddEquations(offset, hailstones[1], hailstones[2]);
	if (offset != answers.Length) throw new Exception();

	var m = new Matrix<double>(matrix);
	var unknowns = m.GaussSolve(answers);
	unknowns.Dump();
	return unknowns;
}