Add a low level quadratic root solver
This is a solver based on W. Kahan's solver. It will be
used as part of Skia's multiple existing solvers.
Change-Id: Ib95c5516f1260e5394a8cc072987c14559bc84b6
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/705225
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 0997033..8152c77 100644
--- a/src/base/SkQuads.cpp
+++ b/src/base/SkQuads.cpp
@@ -6,6 +6,7 @@
*/
#include "src/base/SkQuads.h"
+#include "include/private/base/SkAssert.h"
#include "include/private/base/SkFloatingPoint.h"
#include <cmath>
@@ -92,6 +93,26 @@
return discriminant;
}
+SkQuads::RootResult SkQuads::Roots(double A, double B, double C) {
+ SkASSERT(A != 0);
+
+ const double discriminant = Discriminant(A, B, C);
+
+ if (discriminant == 0) {
+ const double root = B / A;
+ return {discriminant, root, root};
+ }
+
+ if (discriminant > 0) {
+ const double D = sqrt(discriminant);
+ const double R = B > 0 ? B + D : B - D;
+ return {discriminant, R / A, C / R};
+ }
+
+ // The discriminant is negative or is not finite.
+ return {discriminant, NAN, NAN};
+}
+
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 2cdc4ac..ec9918f 100644
--- a/src/base/SkQuads.h
+++ b/src/base/SkQuads.h
@@ -29,6 +29,32 @@
*/
static double Discriminant(double A, double B, double C);
+ struct RootResult {
+ double discriminant;
+ double root0;
+ double root1;
+ };
+
+ /**
+ * Calculate the roots of a quadratic.
+ * Given
+ * A*t^2 -2*B*t + C = 0,
+ * calculate the roots.
+ *
+ * This does not try to detect a linear configuration of the equation, or detect if the two
+ * roots are the same. It returns the discriminant and the two roots.
+ *
+ * Not this uses a different form the quadratic equation to reduce rounding error. Give
+ * standard A, B, C. You can call this root finder with:
+ * Roots(A, -0.5*B, C)
+ * to find the roots of A*x^2 + B*x + C.
+ *
+ * The method used to calculate the roots is from
+ * "On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic"
+ * by W. Kahan.
+ */
+ static RootResult Roots(double A, double B, double C);
+
/**
* Puts up to 2 real solutions to the equation
* A*t^2 + B*t + C = 0
diff --git a/tests/QuadRootsTest.cpp b/tests/QuadRootsTest.cpp
index 34eac57..3397d59 100644
--- a/tests/QuadRootsTest.cpp
+++ b/tests/QuadRootsTest.cpp
@@ -213,3 +213,55 @@
F[1] = F[0];
}
}
+
+DEF_TEST(QuadRoots_Basic, reporter) {
+ {
+ // (x - 1) (x - 1) normal quadratic form A = 1, B = 2, C =1.
+ auto [discriminant, r0, r1] = SkQuads::Roots(1, -0.5 * -2, 1);
+ REPORTER_ASSERT(reporter, discriminant == 0);
+ REPORTER_ASSERT(reporter, r0 == 1 && r1 == 1);
+ }
+
+ {
+ // (x + 2) (x + 2) normal quadratic form A = 1, B = 4, C = 4.
+ auto [discriminant, r0, r1] = SkQuads::Roots(1, -0.5 * 4, 4);
+ REPORTER_ASSERT(reporter, discriminant == 0);
+ REPORTER_ASSERT(reporter, r0 == -2 && r1 == -2);
+ }
+}
+
+// Test the roots using
+// Use quadratics of the form F_n * x^2 - 2 * F_(n-1) * x + F_(n-2).
+// The roots are (F_(n–1) ± 1)/F_n if n is even otherwise there are no roots.
+DEF_TEST(QuadRoots_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;
+ auto [discriminant, r0, r1] = SkQuads::Roots(F[0], F[1], F[2]);
+ REPORTER_ASSERT(reporter, discriminant == expectedDiscriminant);
+
+ // There are only real roots when i is even.
+ if (i % 2 == 0) {
+ const double expectedLittle = ((double)F[1] - 1) / F[0];
+ const double expectedBig = ((double)F[1] + 1) / F[0];
+ if (r0 <= r1) {
+ REPORTER_ASSERT(reporter, r0 == expectedLittle);
+ REPORTER_ASSERT(reporter, r1 == expectedBig);
+ } else {
+ REPORTER_ASSERT(reporter, r1 == expectedLittle);
+ REPORTER_ASSERT(reporter, r0 == expectedBig);
+ }
+ } else {
+ REPORTER_ASSERT(reporter, std::isnan(r0));
+ REPORTER_ASSERT(reporter, std::isnan(r1));
+ }
+
+ F[2] = F[1];
+ F[1] = F[0];
+ }
+}
+
