commit | bbd4d3514eec1981e86226b577887c45a4b01d0c | [log] [tgz] |
---|---|---|

author | Herb Derby <herb@google.com> | Wed May 24 16:49:47 2023 -0400 |

committer | SkCQ <skcq-be@skia-corp.google.com.iam.gserviceaccount.com> | Fri May 26 16:02:32 2023 +0000 |

tree | 706fa150d3427241dbdeeb38964b8a6aa0e7aff7 | |

parent | d707d84b126d103b7e6bbef47d53a9c0d3bb0609 [diff] |

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]; + } +}