2 files changed, 281 insertions, 0 deletions

diff --git a/2023/24/24a.c b/2023/24/24a.c
new file mode 100644
index 0000000..6c75c08
--- /dev/null
+++ b/2023/24/24a.c
@@ -0,0 +1,76 @@
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define N    300
+#define XMIN 2e14
+#define XMAX 4e14
+#define eps  1e-6
+
+#define ABS(x) ((x)>0?(x):-(x))
+
+typedef struct { double x, y, z; } point_t;
+typedef struct { point_t p, v; } line_t;
+
+char *buf, line[N];
+int n, s;
+line_t l[N];
+
+bool isnum(char c) { return (c >= '0' && c <= '9') || c == '-'; }
+
+line_t readl(char *buf) {
+        double a[6];
+        for (int i = 0; i < 6; i++) {
+                a[i] = (double)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]}
+        };
+}
+
+double det(double a, double b, double c, double d) { return a*d-b*c; }
+
+bool parallel2d(line_t l1, line_t l2) {
+        double d = det(l1.v.x, l2.v.x, l1.v.y, l2.v.y);
+        return ABS(d) < eps;
+}
+
+point_t collide2d(line_t l1, line_t l2, double *t1, double *t2) {
+        double a1 = l1.v.x;
+        double b1 = -l2.v.x;
+        double c1 = l2.p.x - l1.p.x;
+        double a2 = l1.v.y;
+        double b2 = -l2.v.y;
+        double c2 = l2.p.y - l1.p.y;
+        
+        double d = det(a1, b1, a2, b2);
+        *t1 = det(c1, b1, c2, b2) / d;
+        *t2 = det(a1, c1, a2, c2) / d;
+
+        double x = l1.p.x + l1.v.x * (*t1);
+        double y = l1.p.y + l1.v.y * (*t1);
+
+        return (point_t) {.x = x, .y = y};
+}
+
+int main() {
+        for (n = 0; (buf = fgets(line, N, stdin)) != NULL; n++)
+                l[n] = readl(buf);
+
+        for (int i = 0; i < n; i++) {
+                for (int j = i+1; j < n; j++) {
+                        if (parallel2d(l[i], l[j])) continue;
+                        double t1, t2;
+                        point_t p = collide2d(l[i], l[j], &t1, &t2);
+                        if (t1 < -eps || t2 < -eps) continue;
+                        s += p.x > XMIN-eps && p.x < XMAX+eps &&
+                             p.y > XMIN-eps && p.y < XMAX+eps;
+                }
+        }
+
+        printf("%d\n", s);
+        return 0;
+}
diff --git a/2023/24/24b.c b/2023/24/24b.c
new file mode 100644
index 0000000..3011108
--- /dev/null
+++ b/2023/24/24b.c
@@ -0,0 +1,205 @@
+/*
+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.
+*/
+
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define D    10
+#define N    300
+
+#define ABS(x) ((x)>0?(x):-(x))
+#define eps 1e1
+#define EQ(x,y) (ABS((x)-(y))<eps)
+
+typedef struct { double 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) {
+        double a[6];
+        for (int i = 0; i < 6; i++) {
+                a[i] = (double)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(double *x, double *y) { double aux = *x; *x = *y; *y = aux; }
+
+void printa(double A[D][D], double C[D], int d) {
+        for (int i = 0; i < d; i++) {
+                printf("[ ");
+                for (int j = 0; j < d; j++)
+                        printf("%.2lf ", A[i][j]);
+                printf("]\t[%.2lf]\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(double A[D][D], double C[D], int d, double *X) {
+        /* Row reduction */
+        for (int i = 0; i < d; i++) {
+                /* Make first nonzero */
+                int nz;
+                for (nz = i; nz < d && EQ(A[i][nz], 0.0); nz++) ;
+                if (nz == d) return false;
+                swap(&C[i], &C[nz]);
+                for (int j = 0; j < d; j++)
+                        swap(&A[i][j], &A[nz][j]);
+
+                /* Reduce rows */
+                for (int ii = i+1; ii < d; ii++) {
+                        double 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;
+}
+
+bool equal(point_t p, point_t q) {
+        return EQ(p.x, q.x) && EQ(p.y, q.y) && EQ(p.z, q.z);
+}
+
+point_t pos(line_t l, double 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++) {
+                double 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 = %.2lf "
+                       "in (%.2lf, %.2lf, %.2lf) and (%.2lf, %.2lf, %.2lf)\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;
+        double 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},
+        };
+        double 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,
+        };
+        double 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 = (%.2lf, %.2lf, %.2lf), v = (%.2lf, %.2lf, %.2lf)\n",
+            sol.p.x, sol.p.y, sol.p.z, sol.v.x, sol.v.y, sol.v.z);
+        printf("%lf\n", sol.p.x + sol.p.y + sol.p.z);
+        return 0;
+}