From e1e7a4b8e86087da79b86ba2ea394785346396db Mon Sep 17 00:00:00 2001
From: Sebastiano Tronto <sebastiano@tronto.net>
Date: Mon, 25 Dec 2023 00:12:46 +0100
Subject: Improved numerical stability by choosing better pivot

---
 2023/24/24b.c | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

(limited to '2023')

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");
-- 
cgit v1.3