Improved numerical stability by choosing better pivot

1 files changed, 20 insertions, 10 deletions

diff --git a/2023/24/24b.c b/2023/24/24b.c
index 3011108..4881099 100644
--- a/2023/24/24b.c
+++ b/2023/24/24b.c
@@ -4,19 +4,21 @@ 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.
+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).
 */
 
 #include <stdbool.h>
 #include <stdio.h>
 #include <stdlib.h>
 
-#define D    10
-#define N    300
+#define D 10
+#define N 300
 
 #define ABS(x) ((x)>0?(x):-(x))
-#define eps 1e1
-#define EQ(x,y) (ABS((x)-(y))<eps)
+#define eps 1e-10
+#define EQ(x,y) (ABS((x)-(y))<eps*(ABS(x)+ABS(y)))
 
 typedef struct { double x, y, z; } point_t;
 typedef struct { point_t p, v; } line_t;
@@ -57,12 +59,18 @@ 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]);
+                int imax;
+                double maxi = 0.0;
+                for (int j = i; j < d; j++) {
+                        if (A[j][i] > maxi) {
+                                maxi = 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[nz][j]);
+                        swap(&A[i][j], &A[imax][j]);
 
                 /* Reduce rows */
                 for (int ii = i+1; ii < d; ii++) {
@@ -106,9 +114,11 @@ void testsolution(line_t sol) {
                 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");