aoc

My solutions for the Advent of Code
git clone https://git.tronto.net/aoc
Download | Log | Files | Refs | README

24b.c (7004B)


      1 /*
      2 Due to numerical errors, the result of this program was not correct on my
      3 machine (ahah - doesn't work on my machine (: )
      4 However, with slightly different values for eps I got two answers whose
      5 difference was 4 and one was too high, the other too low. A quick binary
      6 search gave the correct result.
      7 I'll make a more stable version if I feel like it.
      8 
      9 EDIT: made slighly more numerically stable, now the result is actually
     10 correct (if rounded to the nearest integer).
     11 
     12 EDIT: made numerical result even more accurate using long doubles.
     13 */
     14 
     15 #include <stdbool.h>
     16 #include <stdio.h>
     17 #include <stdlib.h>
     18 
     19 #define D 10
     20 #define N 300
     21 
     22 #define ABS(x) ((x)>0?(x):-(x))
     23 #define eps 1e-10
     24 #define EQ(x,y) (ABS((x)-(y))<eps*(ABS(x)+ABS(y)))
     25 
     26 typedef long double num_t;
     27 typedef struct { num_t x, y, z; } point_t;
     28 typedef struct { point_t p, v; } line_t;
     29 
     30 char *buf, line[N];
     31 int nlines, s;
     32 line_t l[N];
     33 
     34 bool isnum(char c) { return (c >= '0' && c <= '9') || c == '-'; }
     35 
     36 line_t readl(char *buf) {
     37 	num_t a[6];
     38 	for (int i = 0; i < 6; i++) {
     39 		a[i] = (num_t)atoll(buf);
     40 		while (isnum(*buf)) buf++;
     41 		while (*buf != '\n' && !isnum(*buf)) buf++;
     42 	}
     43 	return (line_t) {
     44 	    .p = (point_t){.x = a[0], .y = a[1], .z = a[2]},
     45 	    .v = (point_t){.x = a[3], .y = a[4], .z = a[5]}
     46 	};
     47 }
     48 
     49 void swap(num_t *x, num_t *y) { num_t aux = *x; *x = *y; *y = aux; }
     50 
     51 void printa(num_t A[D][D], num_t C[D], int d) {
     52 	for (int i = 0; i < d; i++) {
     53 		printf("[ ");
     54 		for (int j = 0; j < d; j++)
     55 			printf("%.2llf ", A[i][j]);
     56 		printf("]\t[%.2llf]\n", C[i]);
     57 	}
     58 }
     59 
     60 /* Solve the linear system AX = C with row reduction */
     61 /* Return true if it has a unique solution, false in any other case */
     62 bool solvesystem(num_t A[D][D], num_t C[D], int d, num_t *X) {
     63 	/* Row reduction */
     64 	for (int i = 0; i < d; i++) {
     65 		/* Make first nonzero */
     66 		int imax;
     67 		num_t maxi = 0.0;
     68 		for (int j = i; j < d; j++) {
     69 			if (ABS(A[j][i]) > maxi) {
     70 				maxi = ABS(A[j][i]);
     71 				imax = j;
     72 			}
     73 		}
     74 		if (EQ(maxi, 0.0)) return false;
     75 		swap(&C[i], &C[imax]);
     76 		for (int j = 0; j < d; j++)
     77 			swap(&A[i][j], &A[imax][j]);
     78 
     79 		/* Reduce rows */
     80 		for (int ii = i+1; ii < d; ii++) {
     81 			num_t r = A[ii][i] / A[i][i];
     82 			for (int k = i; k < d; k++)
     83 				A[ii][k] -= r*A[i][k];
     84 			C[ii] -= r*C[i];
     85 		}
     86 	}
     87 
     88 	/* Back substitution */
     89 	for (int i = d-1; i >= 0; i--) {
     90 		X[i] = C[i];
     91 		for (int j = i+1; j < d; j++)
     92 			X[i] -= A[i][j]*X[j];
     93 		X[i] /= A[i][i];
     94 	}
     95 	
     96 	return true;
     97 }
     98 
     99 point_t pos(line_t l, num_t t) {
    100 	return (point_t) {
    101 	    .x = l.p.x + t*l.v.x,
    102 	    .y = l.p.y + t*l.v.y,
    103 	    .z = l.p.z + t*l.v.z
    104 	};
    105 }
    106 
    107 void testsolution(line_t sol) {
    108 	for (int i = 0; i < nlines; i++) {
    109 		num_t t = -1.0;
    110 		if (!EQ(sol.v.x, l[i].v.x))
    111 			t = (sol.p.x-l[i].p.x)/(l[i].v.x-sol.v.x);
    112 		if (!EQ(sol.v.y, l[i].v.y))
    113 			t = (sol.p.y-l[i].p.y)/(l[i].v.y-sol.v.y);
    114 		if (!EQ(sol.v.z, l[i].v.z))
    115 			t = (sol.p.z-l[i].p.z)/(l[i].v.z-sol.v.z);
    116 		point_t p1 = pos(sol, t), p2 = pos(l[i], t);
    117 		/*
    118 		printf("Line %d intersected at t = %.2llf in "
    119 		       "(%.2llf, %.2llf, %.2llf) and "
    120 		       "(%.2llf, %.2llf, %.2llf)\n",
    121 		       i, t, p1.x, p1.y, p1.z, p2.x, p2.y, p2.z);
    122 		*/
    123 		if (t < eps ||
    124 		    !EQ(p1.x, p2.x) || !EQ(p1.y, p2.y) || !EQ(p1.z, p2.z)) {
    125 			printf("Error!\n");
    126 			return;
    127 		}
    128 	}
    129 	printf("All lines intersected correctly, solution is valid\n");
    130 }
    131 
    132 int main() {
    133 	for (nlines = 0; (buf = fgets(line, N, stdin)) != NULL; nlines++)
    134 		l[nlines] = readl(buf);
    135 
    136 	/* To figure the correct starting point + velocity, we solve some
    137 	 * systems of linear equations. Write your starting point and velocity
    138 	 * as unknowns x, y, z, Vx, Vy, Vz. Equating the position of the
    139 	 * rock at time t1 (another unknown parameter) with the position
    140 	 * of one of the hailstones at the same time t1, we get a system
    141 	 * of 3 equations and 7 unknowns. Unfortunately, these equations
    142 	 * have degree 2 - this is not a linear system!
    143 	 * However, manipulating these equations a bit we can get a linear
    144 	 * equation of the type:
    145 	 *     (Vy1-Vy2)x - (Vx1-Vx2)y + - (y1-y2)Vx + (x1-x2)Vy =
    146 	 *         y2Vx2 + x2Vy2 - y1Vx1 - x1Vy1
    147 	 * Where x1, y1, z2, Vx1, Vy1, Vz1 and x2, y2, z2, Vx2, Vy2, Vz2
    148 	 * are the starting points and velocities of two of the hailstones.
    149 	 * So with 2 lines we can get a linear equation. Similarly, we can
    150 	 * get equations involving the unknowns z and Vz.
    151 	 * We can use the myriad of hailstones we have to generate as many
    152 	 * equations as we like. The system is going to be overdetermined,
    153 	 * but the problem statement seems to ensure that there is going to
    154 	 * be a solution. On the other hand it can happen that we make a
    155 	 * bad choice of lines and the equation we use are underdetermined.
    156 	 * This last problem is not accounted for in the code - if it happens,
    157 	 * one can shuffle the input file until it works.
    158 	 */
    159 	int d = 6;
    160 	num_t A[D][D] = {
    161 		/* First equation: lines 0 and 1, x and y only */
    162 		{l[0].v.y-l[1].v.y, l[1].v.x-l[0].v.x, 0.0,
    163 		    l[1].p.y-l[0].p.y, l[0].p.x-l[1].p.x, 0.0},
    164 		/* Second equation: lines 0 and 2, x and y only */
    165 		{l[0].v.y-l[2].v.y, l[2].v.x-l[0].v.x, 0.0,
    166 		    l[2].p.y-l[0].p.y, l[0].p.x-l[2].p.x, 0.0},
    167 		/* Third equation: lines 0 and 1, x and z only */
    168 		{l[0].v.z-l[1].v.z, 0.0, l[1].v.x-l[0].v.x,
    169 		    l[1].p.z-l[0].p.z, 0.0, l[0].p.x-l[1].p.x},
    170 		/* Fourth equation: lines 0 and 2, x and z only */
    171 		{l[0].v.z-l[2].v.z, 0.0, l[2].v.x-l[0].v.x,
    172 		    l[2].p.z-l[0].p.z, 0.0, l[0].p.x-l[2].p.x},
    173 		/* Fifth equation: lines 0 and 1, y and z only */
    174 		{0.0, l[0].v.z-l[1].v.z, l[1].v.y-l[0].v.y,
    175 		    0.0, l[1].p.z-l[0].p.z, l[0].p.y-l[1].p.y},
    176 		/* Sixth equation: lines 0 and 2, y and z only */
    177 		{0.0, l[0].v.z-l[2].v.z, l[2].v.y-l[0].v.y,
    178 		    0.0, l[2].p.z-l[0].p.z, l[0].p.y-l[2].p.y},
    179 	};
    180 	num_t C[D] = {
    181 		/* First equation: lines 0 and 1, x and y only */
    182 		l[1].p.y*l[1].v.x - l[1].p.x*l[1].v.y 
    183 		    - l[0].p.y*l[0].v.x + l[0].p.x*l[0].v.y,
    184 		/* Second equation: lines 0 and 2, x and y only */
    185 		l[2].p.y*l[2].v.x - l[2].p.x*l[2].v.y 
    186 		    - l[0].p.y*l[0].v.x + l[0].p.x*l[0].v.y,
    187 		/* Third equation: lines 0 and 1, x and z only */
    188 		l[1].p.z*l[1].v.x - l[1].p.x*l[1].v.z 
    189 		    - l[0].p.z*l[0].v.x + l[0].p.x*l[0].v.z,
    190 		/* Fourth equation: lines 0 and 2, x and z only */
    191 		l[2].p.z*l[2].v.x - l[2].p.x*l[2].v.z 
    192 		    - l[0].p.z*l[0].v.x + l[0].p.x*l[0].v.z,
    193 		/* Fifth equation: lines 0 and 1, y and z only */
    194 		l[1].p.z*l[1].v.y - l[1].p.y*l[1].v.z 
    195 		    - l[0].p.z*l[0].v.y + l[0].p.y*l[0].v.z,
    196 		/* Sixth equation: lines 0 and 2, y and z only */
    197 		l[2].p.z*l[2].v.y - l[2].p.y*l[2].v.z 
    198 		    - l[0].p.z*l[0].v.y + l[0].p.y*l[0].v.z,
    199 	};
    200 	num_t X[d];
    201 	if (!solvesystem(A, C, d, X)) {
    202 		printf("No unique solution, shuffle input and try again.\n");
    203 		exit(1);
    204 	}
    205 
    206 	line_t sol = {
    207 		.p = {.x = X[0], .y = X[1], .z = X[2]},
    208 		.v = {.x = X[3], .y = X[4], .z = X[5]}
    209 	};
    210 	testsolution(sol);
    211 
    212 	printf("p = (%.2llf, %.2llf, %.2llf), v = (%.2llf, %.2llf, %.2llf)\n",
    213 	    sol.p.x, sol.p.y, sol.p.z, sol.v.x, sol.v.y, sol.v.z);
    214 	printf("%llf\n", sol.p.x + sol.p.y + sol.p.z);
    215 	return 0;
    216 }