Calculate a very accurate discriminant
Calculate the discriminant in the ordinary manner. Check
for too much cancellation if B^2 and AC are too close.
If they are too close, then track the rounding error while
recalculating the discriminant, and add rounding error
back in. Tests included.
Change-Id: If2a04eb5c663ffe05aa90df9116796bc5f5d7a64
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/703078
Commit-Queue: Herb Derby <herb@google.com>
Reviewed-by: Jim Van Verth <jvanverth@google.com>
diff --git a/src/base/SkQuads.cpp b/src/base/SkQuads.cpp
index a778379..acdfabf 100644
--- a/src/base/SkQuads.cpp
+++ b/src/base/SkQuads.cpp
@@ -40,6 +40,33 @@
return std::abs(A / B) < 1.0e-16;
}
+double SkQuads::Discriminant(const double a, const double b, const double c) {
+ const double b2 = b * b;
+ const double ac = a * c;
+
+ // Calculate the rough discriminate which may suffer from a loss in precision due to b2 and
+ // ac being too close.
+ const double roughDiscriminant = b2 - ac;
+
+ // Check if b2 and ac were too close, and caused catastrophic cancellation. The
+ // roughDiscriminant should be good enough most of the time. If b2 and ac are very close to
+ // each other, then roughDiscriminant will be very small with most of the bits canceled
+ // making much smaller than roundOffCheck.
+ const double roundOffCheck = b2 + ac;
+ if (3 * std::abs(roughDiscriminant) >= roundOffCheck) {
+ return roughDiscriminant;
+ }
+
+ // Use the extra internal precision afforded by fma to calculate the rounding error for
+ // b^2 and ac.
+ const double b2RoundingError = std::fma(b, b, -b2);
+ const double acRoundingError = std::fma(a, c, -ac);
+
+ // Add the total rounding error back into the discriminant guess.
+ const double discriminant = (b2 - ac) + (b2RoundingError - acRoundingError);
+ return discriminant;
+}
+
int SkQuads::RootsReal(const double A, const double B, const double C, double solution[2]) {
if (close_to_linear(A, B)) {
return solve_linear(B, C, solution);
diff --git a/src/base/SkQuads.h b/src/base/SkQuads.h
index 645d43b..2cdc4ac 100644
--- a/src/base/SkQuads.h
+++ b/src/base/SkQuads.h
@@ -15,12 +15,26 @@
class SkQuads {
public:
/**
+ * Calculate a very accurate discriminant.
+ * Given
+ * A*t^2 -2*B*t + C = 0,
+ * calculate
+ * B^2 - AC
+ * accurate to 2 bits.
+ * Note the form of the quadratic is slightly different from the normal formulation.
+ *
+ * The method used to calculate the discriminant is from
+ * "On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic"
+ * by W. Kahan.
+ */
+ static double Discriminant(double A, double B, double C);
+
+ /**
* Puts up to 2 real solutions to the equation
* A*t^2 + B*t + C = 0
* in the provided array.
*/
- static int RootsReal(double A, double B, double C,
- double solution[2]);
+ static int RootsReal(double A, double B, double C, double solution[2]);
/**
* Evaluates the quadratic function with the 3 provided coefficients and the
diff --git a/tests/QuadRootsTest.cpp b/tests/QuadRootsTest.cpp
index 62a18ad..34eac57 100644
--- a/tests/QuadRootsTest.cpp
+++ b/tests/QuadRootsTest.cpp
@@ -13,9 +13,10 @@
#include "tests/Test.h"
#include <algorithm>
-#include <cstddef>
#include <cfloat>
#include <cmath>
+#include <cstddef>
+#include <cstdint>
#include <iterator>
#include <string>
@@ -194,3 +195,21 @@
"Nan constant"
);
}
+
+// Test the discriminant using
+// Use quadratics of the form F_n * x^2 - 2 * F_(n-1) * x + F_(n-2).
+// This has a discriminant of F_(n-1)^2 - F_n * F_(n-2) = 1 if n is even else -1.
+DEF_TEST(QuadDiscriminant_Fibonacci, reporter) {
+ // n, n-1, n-2
+ int64_t F[] = {1, 1, 0};
+ // F_79 just fits in the 53 significant bits of a double.
+ for (int i = 2; i < 79; ++i) {
+ F[0] = F[1] + F[2];
+
+ const int expectedDiscriminant = i % 2 == 0 ? 1 : -1;
+ REPORTER_ASSERT(reporter, SkQuads::Discriminant(F[0], F[1], F[2]) == expectedDiscriminant);
+
+ F[2] = F[1];
+ F[1] = F[0];
+ }
+}
