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 }