aoc

24b.c (7004B)

---

 /*
 Due to numerical errors, the result of this program was not correct on my
 machine (ahah - doesn't work on my machine (: )
 However, with slightly different values for eps I got two answers whose
 difference was 4 and one was too high, the other too low. A quick binary
 search gave the correct result.
 I'll make a more stable version if I feel like it.
 
 EDIT: made slighly more numerically stable, now the result is actually
 correct (if rounded to the nearest integer).
 
 EDIT: made numerical result even more accurate using long doubles.
 */
 
 #include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
 
 #define D 10
 #define N 300
 
 #define ABS(x) ((x)>0?(x):-(x))
 #define eps 1e-10
 #define EQ(x,y) (ABS((x)-(y))<eps*(ABS(x)+ABS(y)))
 
 typedef long double num_t;
 typedef struct { num_t x, y, z; } point_t;
 typedef struct { point_t p, v; } line_t;
 
 char *buf, line[N];
 int nlines, s;
 line_t l[N];
 
 bool isnum(char c) { return (c >= '0' && c <= '9') || c == '-'; }
 
 line_t readl(char *buf) {
 	num_t a[6];
 	for (int i = 0; i < 6; i++) {
 		a[i] = (num_t)atoll(buf);
 		while (isnum(*buf)) buf++;
 		while (*buf != '\n' && !isnum(*buf)) buf++;
 	}
 	return (line_t) {
 	    .p = (point_t){.x = a[0], .y = a[1], .z = a[2]},
 	    .v = (point_t){.x = a[3], .y = a[4], .z = a[5]}
 	};
 }
 
 void swap(num_t *x, num_t *y) { num_t aux = *x; *x = *y; *y = aux; }
 
 void printa(num_t A[D][D], num_t C[D], int d) {
 	for (int i = 0; i < d; i++) {
 		printf("[ ");
 		for (int j = 0; j < d; j++)
 			printf("%.2llf ", A[i][j]);
 		printf("]\t[%.2llf]\n", C[i]);
 	}
 }
 
 /* Solve the linear system AX = C with row reduction */
 /* Return true if it has a unique solution, false in any other case */
 bool solvesystem(num_t A[D][D], num_t C[D], int d, num_t *X) {
 	/* Row reduction */
 	for (int i = 0; i < d; i++) {
 		/* Make first nonzero */
 		int imax;
 		num_t maxi = 0.0;
 		for (int j = i; j < d; j++) {
 			if (ABS(A[j][i]) > maxi) {
 				maxi = ABS(A[j][i]);
 				imax = j;
 			}
 		}
 		if (EQ(maxi, 0.0)) return false;
 		swap(&C[i], &C[imax]);
 		for (int j = 0; j < d; j++)
 			swap(&A[i][j], &A[imax][j]);
 
 		/* Reduce rows */
 		for (int ii = i+1; ii < d; ii++) {
 			num_t r = A[ii][i] / A[i][i];
 			for (int k = i; k < d; k++)
 				A[ii][k] -= r*A[i][k];
 			C[ii] -= r*C[i];
 		}
 	}
 
 	/* Back substitution */
 	for (int i = d-1; i >= 0; i--) {
 		X[i] = C[i];
 		for (int j = i+1; j < d; j++)
 			X[i] -= A[i][j]*X[j];
 		X[i] /= A[i][i];
 	}
 	
 	return true;
 }
 
 point_t pos(line_t l, num_t t) {
 	return (point_t) {
 	    .x = l.p.x + t*l.v.x,
 	    .y = l.p.y + t*l.v.y,
 	    .z = l.p.z + t*l.v.z
 	};
 }
 
 void testsolution(line_t sol) {
 	for (int i = 0; i < nlines; i++) {
 		num_t t = -1.0;
 		if (!EQ(sol.v.x, l[i].v.x))
 			t = (sol.p.x-l[i].p.x)/(l[i].v.x-sol.v.x);
 		if (!EQ(sol.v.y, l[i].v.y))
 			t = (sol.p.y-l[i].p.y)/(l[i].v.y-sol.v.y);
 		if (!EQ(sol.v.z, l[i].v.z))
 			t = (sol.p.z-l[i].p.z)/(l[i].v.z-sol.v.z);
 		point_t p1 = pos(sol, t), p2 = pos(l[i], t);
 		/*
 		printf("Line %d intersected at t = %.2llf in "
 		       "(%.2llf, %.2llf, %.2llf) and "
 		       "(%.2llf, %.2llf, %.2llf)\n",
 		       i, t, p1.x, p1.y, p1.z, p2.x, p2.y, p2.z);
 		*/
 		if (t < eps ||
 		    !EQ(p1.x, p2.x) || !EQ(p1.y, p2.y) || !EQ(p1.z, p2.z)) {
 			printf("Error!\n");
 			return;
 		}
 	}
 	printf("All lines intersected correctly, solution is valid\n");
 }
 
 int main() {
 	for (nlines = 0; (buf = fgets(line, N, stdin)) != NULL; nlines++)
 		l[nlines] = readl(buf);
 
 	/* To figure the correct starting point + velocity, we solve some
 	 * systems of linear equations. Write your starting point and velocity
 	 * as unknowns x, y, z, Vx, Vy, Vz. Equating the position of the
 	 * rock at time t1 (another unknown parameter) with the position
 	 * of one of the hailstones at the same time t1, we get a system
 	 * of 3 equations and 7 unknowns. Unfortunately, these equations
 	 * have degree 2 - this is not a linear system!
 	 * However, manipulating these equations a bit we can get a linear
 	 * equation of the type:
 	 *     (Vy1-Vy2)x - (Vx1-Vx2)y + - (y1-y2)Vx + (x1-x2)Vy =
 	 *         y2Vx2 + x2Vy2 - y1Vx1 - x1Vy1
 	 * Where x1, y1, z2, Vx1, Vy1, Vz1 and x2, y2, z2, Vx2, Vy2, Vz2
 	 * are the starting points and velocities of two of the hailstones.
 	 * So with 2 lines we can get a linear equation. Similarly, we can
 	 * get equations involving the unknowns z and Vz.
 	 * We can use the myriad of hailstones we have to generate as many
 	 * equations as we like. The system is going to be overdetermined,
 	 * but the problem statement seems to ensure that there is going to
 	 * be a solution. On the other hand it can happen that we make a
 	 * bad choice of lines and the equation we use are underdetermined.
 	 * This last problem is not accounted for in the code - if it happens,
 	 * one can shuffle the input file until it works.
 	 */
 	int d = 6;
 	num_t A[D][D] = {
 		/* First equation: lines 0 and 1, x and y only */
 		{l[0].v.y-l[1].v.y, l[1].v.x-l[0].v.x, 0.0,
 		    l[1].p.y-l[0].p.y, l[0].p.x-l[1].p.x, 0.0},
 		/* Second equation: lines 0 and 2, x and y only */
 		{l[0].v.y-l[2].v.y, l[2].v.x-l[0].v.x, 0.0,
 		    l[2].p.y-l[0].p.y, l[0].p.x-l[2].p.x, 0.0},
 		/* Third equation: lines 0 and 1, x and z only */
 		{l[0].v.z-l[1].v.z, 0.0, l[1].v.x-l[0].v.x,
 		    l[1].p.z-l[0].p.z, 0.0, l[0].p.x-l[1].p.x},
 		/* Fourth equation: lines 0 and 2, x and z only */
 		{l[0].v.z-l[2].v.z, 0.0, l[2].v.x-l[0].v.x,
 		    l[2].p.z-l[0].p.z, 0.0, l[0].p.x-l[2].p.x},
 		/* Fifth equation: lines 0 and 1, y and z only */
 		{0.0, l[0].v.z-l[1].v.z, l[1].v.y-l[0].v.y,
 		    0.0, l[1].p.z-l[0].p.z, l[0].p.y-l[1].p.y},
 		/* Sixth equation: lines 0 and 2, y and z only */
 		{0.0, l[0].v.z-l[2].v.z, l[2].v.y-l[0].v.y,
 		    0.0, l[2].p.z-l[0].p.z, l[0].p.y-l[2].p.y},
 	};
 	num_t C[D] = {
 		/* First equation: lines 0 and 1, x and y only */
 		l[1].p.y*l[1].v.x - l[1].p.x*l[1].v.y 
 		    - l[0].p.y*l[0].v.x + l[0].p.x*l[0].v.y,
 		/* Second equation: lines 0 and 2, x and y only */
 		l[2].p.y*l[2].v.x - l[2].p.x*l[2].v.y 
 		    - l[0].p.y*l[0].v.x + l[0].p.x*l[0].v.y,
 		/* Third equation: lines 0 and 1, x and z only */
 		l[1].p.z*l[1].v.x - l[1].p.x*l[1].v.z 
 		    - l[0].p.z*l[0].v.x + l[0].p.x*l[0].v.z,
 		/* Fourth equation: lines 0 and 2, x and z only */
 		l[2].p.z*l[2].v.x - l[2].p.x*l[2].v.z 
 		    - l[0].p.z*l[0].v.x + l[0].p.x*l[0].v.z,
 		/* Fifth equation: lines 0 and 1, y and z only */
 		l[1].p.z*l[1].v.y - l[1].p.y*l[1].v.z 
 		    - l[0].p.z*l[0].v.y + l[0].p.y*l[0].v.z,
 		/* Sixth equation: lines 0 and 2, y and z only */
 		l[2].p.z*l[2].v.y - l[2].p.y*l[2].v.z 
 		    - l[0].p.z*l[0].v.y + l[0].p.y*l[0].v.z,
 	};
 	num_t X[d];
 	if (!solvesystem(A, C, d, X)) {
 		printf("No unique solution, shuffle input and try again.\n");
 		exit(1);
 	}
 
 	line_t sol = {
 		.p = {.x = X[0], .y = X[1], .z = X[2]},
 		.v = {.x = X[3], .y = X[4], .z = X[5]}
 	};
 	testsolution(sol);
 
 	printf("p = (%.2llf, %.2llf, %.2llf), v = (%.2llf, %.2llf, %.2llf)\n",
 	    sol.p.x, sol.p.y, sol.p.z, sol.v.x, sol.v.y, sol.v.z);
 	printf("%llf\n", sol.p.x + sol.p.y + sol.p.z);
 	return 0;
 }