Fix findCubicRoots for cases where coefficients are close to zero

The equation is a*x^3 + b*x^2 + c*x + d = 0.
Previously, we would divide by zero if a = ~0.
This change also makes sure that we return zero no roots in the case
where a = ~0 and c*c - 4*b*d < 0, and the case where a = b = c = ~0.
Finally, we return 0 or 1 if we're close enough to assume that it could
be a numerical error.

This change also adds tests for the above cases.

Change-Id: I426d2fc6175b3aff6fe099845bf63d433c158536
Co-authored-by: Christian Strømme <christian.stromme@qt.io>
Reviewed-by: Christian Stromme <christian.stromme@qt.io>

2 files changed, 136 insertions, 5 deletions

diff --git a/src/animation/backend/bezierevaluator.cpp b/src/animation/backend/bezierevaluator.cpp
index 9d81a2f2f..2752904fc 100644
--- a/src/animation/backend/bezierevaluator.cpp
+++ b/src/animation/backend/bezierevaluator.cpp
@@ -161,6 +161,13 @@ float BezierEvaluator::parameterForTime(float time) const
     return 0.0f;
 }
 
+bool almostZero(float value, float threshold=1e-3f)
+{
+    // 1e-3 might seem excessively fuzzy, but any smaller value will make the
+    // factors a, b, and c large enough to knock out the cubic solver.
+    return value > -threshold && value < threshold;
+}
+
 /*!
     \internal
 
@@ -171,17 +178,45 @@ float BezierEvaluator::parameterForTime(float time) const
  */
 int BezierEvaluator::findCubicRoots(const float coeffs[4], float roots[3])
 {
+    const float a = coeffs[3];
+    const float b = coeffs[2];
+    const float c = coeffs[1];
+    const float d = coeffs[0];
+
+    // Simple cases with linear, quadratic or invalid equations
+    if (almostZero(a)) {
+        if (almostZero(b)) {
+            if (almostZero(c))
+                return 0;
+
+            roots[0] = -d / c;
+            return 1;
+        }
+        const float discriminant = c * c - 4.f * b * d;
+        if (discriminant < 0.f)
+            return 0;
+
+        if (discriminant == 0.f) {
+            roots[0] = -c / (2.f * b);
+            return 1;
+        }
+
+        roots[0] = (-c + std::sqrt(discriminant)) / (2.f * b);
+        roots[1] = (-c - std::sqrt(discriminant)) / (2.f * b);
+        return 2;
+    }
+
     // See https://en.wikipedia.org/wiki/Cubic_function#General_solution_to_the_cubic_equation_with_real_coefficients
     // for a description. We depress the general cubic to a form that can more easily be solved. Solve it and then
     // substitue the results back to get the roots of the original cubic.
-    int numberOfRoots;
+    int numberOfRoots = 0;
     const double oneThird = 1.0 / 3.0;
     const double piByThree = M_PI / 3.0;
 
     // Put cubic into normal format: x^3 + Ax^2 + Bx + C = 0
-    const double A = coeffs[ 2 ] / coeffs[ 3 ];
-    const double B = coeffs[ 1 ] / coeffs[ 3 ];
-    const double C = coeffs[ 0 ] / coeffs[ 3 ];
+    const double A = double(b / a);
+    const double B = double(c / a);
+    const double C = double(d / a);
 
     // Substitute x = y - A/3 to eliminate quadratic term (depressed form):
     // x^3 + px + q = 0
@@ -226,8 +261,14 @@ int BezierEvaluator::findCubicRoots(const float coeffs[4], float roots[3])
 
     // Substitute back in
     const double sub = oneThird * A;
-    for (int i = 0; i < numberOfRoots; ++i)
+    for (int i = 0; i < numberOfRoots; ++i) {
         roots[i] -= sub;
+        // Take care of cases where we are close to zero or one
+        if (almostZero(roots[i], 1e-6f))
+            roots[i] = 0.f;
+        if (almostZero(roots[i] - 1.f, 1e-6f))
+            roots[i] = 1.f;
+    }
 
     return numberOfRoots;
 }
diff --git a/tests/auto/animation/bezierevaluator/tst_bezierevaluator.cpp b/tests/auto/animation/bezierevaluator/tst_bezierevaluator.cpp
index 5dc971ab7..b2449ea1d 100644
--- a/tests/auto/animation/bezierevaluator/tst_bezierevaluator.cpp
+++ b/tests/auto/animation/bezierevaluator/tst_bezierevaluator.cpp
@@ -97,6 +97,96 @@ private Q_SLOTS:
         roots[2] = 2.0f - std::sqrt(7.0f);
         QTest::newRow("a=1, b=-5, c=1, d=3") << a << b << c << d << rootCount << roots;
         roots.clear();
+
+        // quadratic equation
+        a = 0.0f;
+        b = 9.0f;
+        c = 11.0f;
+        d = 3.0f;
+        roots.resize(2);
+        roots[0] = -11.0f/18.0f + std::sqrt(13.0f) / 18.0f;
+        roots[1] = -11.0f/18.0f - std::sqrt(13.0f) / 18.0f;
+        QTest::newRow("a=0, b=9, c=11, d=3") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // quadratic equation with discriminant = 0
+        a = 0.0f;
+        b = 1.0f;
+        c = 2.0f;
+        d = 1.0f;
+        roots.resize(1);
+        roots[0] = -1.f;
+        QTest::newRow("a=0, b=1, c=2, d=1") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // quadratic equation with discriminant < 0
+        a = 0.0f;
+        b = 1.0f;
+        c = 4.0f;
+        d = 8.0f;
+        roots.resize(0);
+        QTest::newRow("a=0, b=1, c=4, d=8") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // linear equation
+        a = 0.0f;
+        b = 0.0f;
+        c = 2.0f;
+        d = 1.0f;
+        roots.resize(1);
+        roots[0] = -0.5f;
+        QTest::newRow("a=0, b=0, c=2, d=1") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // linear equation
+        a = 0.0f;
+        b = 0.0f;
+        c = 8.0f;
+        d = -5.0f;
+        roots.resize(1);
+        roots[0] = -d/c;
+        QTest::newRow("a=0, b=0, c=8, d=-5") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // invalid equation
+        a = 0.0f;
+        b = 0.0f;
+        c = 0.0f;
+        d = -5.0f;
+        roots.resize(0);
+        QTest::newRow("a=0, b=0, c=0, d=-5") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // Invalid equation
+        a = 0.0f;
+        b = 0.0f;
+        c = 0.0f;
+        d = 42.0f;
+        roots.resize(0);
+        QTest::newRow("a=0, b=0, c=0, d=42") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // almost linear equation
+        a = 1.90735e-06f;
+        b = -2.86102e-06f;
+        c = 5.0;
+        d = 0.0;
+        roots.resize(1);
+        roots[0] = -d/c;
+        QTest::newRow("a=~0, b=~0, c=5, d=0") << a << b << c << d << roots.size() << roots;
+        roots.clear();
+
+        // case that produces a result just below zero, that should be evaluated as zero
+        a = -0.75f;
+        b = 0.75f;
+        c = 2.5;
+        d = 0.0;
+        roots.resize(3);
+        roots[0] = 2.39297f;
+        roots[1] = 0.f;
+        roots[2] = -1.39297f;
+        QTest::newRow("a=-0.75, b=0.75, c=2.5, d=0") << a << b << c << d << roots.size() << roots;
+        roots.clear();
     }
 
     void checkFindCubicRoots()