From ac31fa35fa1c0b3533fb7242222ef714b6f861cc Mon Sep 17 00:00:00 2001
From: Robin Leroy <egg.robin.leroy@gmail.com>
Date: Wed, 18 Sep 2024 10:47:58 +0200
Subject: [PATCH] That is not how cancellation works.

---
 numerics/double_precision_body.hpp |  1 +
 numerics/root_finders_body.hpp     | 12 ++++++------
 numerics/root_finders_test.cpp     | 14 ++++++++------
 3 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/numerics/double_precision_body.hpp b/numerics/double_precision_body.hpp
index 2ddffcd1d9..6033fd21dc 100644
--- a/numerics/double_precision_body.hpp
+++ b/numerics/double_precision_body.hpp
@@ -223,6 +223,7 @@ constexpr DoublePrecision<Product<T, U>> VeltkampDekkerProduct(T const& a,
   U const hy = y - py + py;
   U const ty = y - hy;
   // Veltkamp’s 1968 algorithm, as given in [Dek71, p. 234].
+  // See also exactmul in [Lin81, p. 278].
   z = x * y;
   zz = (((hx * hy - z) + hx * ty) + tx * hy) + tx * ty;
   // Dekker’s algorithm (5.12) would be
diff --git a/numerics/root_finders_body.hpp b/numerics/root_finders_body.hpp
index 23aa0a6797..71e465c98f 100644
--- a/numerics/root_finders_body.hpp
+++ b/numerics/root_finders_body.hpp
@@ -356,13 +356,13 @@ BoundedArray<Argument, 2> SolveQuadraticEquation(
   // Use double precision for the discriminant because there can be
   // cancellations.  Higham mentions that it is necessary “to use extended
   // precision (or some trick tantamount to the use of extended precision).”
-  DoublePrecision<Discriminant> discriminant = TwoSum(a₁ * a₁, -4.0 * a₀ * a₂);
+  DoublePrecision<Discriminant> discriminant =
+      TwoProduct(a₁, a₁) - TwoProduct(4.0 * a₀, a₂);
 
-  if (discriminant.value == zero && discriminant.error == zero) {
+  if (discriminant.value == zero) {
     // One solution.
     return {origin - 0.5 * a₁ / a₂};
-  } else if (discriminant.value < zero ||
-             (discriminant.value == zero && discriminant.error < zero)) {
+  } else if (discriminant.value < zero) {
     // No solution.
     return {};
   } else {
@@ -371,9 +371,9 @@ BoundedArray<Argument, 2> SolveQuadraticEquation(
     Derivative1 numerator;
     constexpr Derivative1 zero{};
     if (a₁ > zero) {
-      numerator = -a₁ - Sqrt(discriminant.value + discriminant.error);
+      numerator = -a₁ - Sqrt(discriminant.value);
     } else {
-      numerator = -a₁ + Sqrt(discriminant.value + discriminant.error);
+      numerator = -a₁ + Sqrt(discriminant.value);
     }
     auto const x₁ = origin + numerator / (2.0 * a₂);
     auto const x₂ = origin + (2.0 * a₀) / numerator;
diff --git a/numerics/root_finders_test.cpp b/numerics/root_finders_test.cpp
index 0a40a74bc0..4ec4de1114 100644
--- a/numerics/root_finders_test.cpp
+++ b/numerics/root_finders_test.cpp
@@ -469,14 +469,16 @@ TEST_F(RootFindersTest, QuadraticEquations) {
   EXPECT_THAT(s3, ElementsAre(-1.0));
 
   // An ill-conditioned system.  I fart in its general direction.  If done
-  // naively, this yields {-100032., -99968.4} according to Mathematica.
+  // naïvely, the results are -100000.0031'0988855 and -99999.9968'90111535
+  // (separator after the last correct decimal digit), whose error is
+  // approximately 3.5 million ULPs.
   auto const s4 = SolveQuadraticEquation(0.0,
-                                         1.0000001e25,
-                                         2.0000003e20,
-                                         1.0000001e15);
+                                         1.000'000'000'000'001e25,
+                                         2.000'000'000'000'003e20,
+                                         1.000'000'000'000'001e15);
   EXPECT_THAT(s4,
-              ElementsAre(AlmostEquals(-100031.62777541532972762902, 66),
-                          AlmostEquals(-99968.38222458367027247098, 65)));
+              ElementsAre(AlmostEquals(-100000.00316153381194436811, 0),
+                          AlmostEquals(-99999.996838466282967631886, 1)));
 
   // A typed system.
   Instant const t0;