src/base/SkCubics.cpp

/*
 * Copyright 2023 Google LLC
 *
 * Use of this source code is governed by a BSD-style license that can be
 * found in the LICENSE file.
 */

#include "src/base/SkCubics.h"

#include "include/private/base/SkAssert.h"
#include "include/private/base/SkFloatingPoint.h"
#include "include/private/base/SkTPin.h"
#include "src/base/SkQuads.h"

#include <algorithm>
#include <cmath>

static constexpr double PI = 3.141592653589793;

static bool nearly_equal(double x, double y) {
    if (sk_double_nearly_zero(x)) {
        return sk_double_nearly_zero(y);
    }
    return sk_doubles_nearly_equal_ulps(x, y);
}

// When the A coefficient of a cubic is close to 0, there can be floating point error
// that arises from computing a very large root. In those cases, we would rather be
// precise about the smaller 2 roots, so we have this arbitrary cutoff for when A is
// really small or small compared to B.
static bool close_to_a_quadratic(double A, double B) {
    if (sk_double_nearly_zero(B)) {
        return sk_double_nearly_zero(A);
    }
    return std::abs(A / B) < 1.0e-7;
}

int SkCubics::RootsReal(double A, double B, double C, double D, double solution[3]) {
    if (close_to_a_quadratic(A, B)) {
        return SkQuads::RootsReal(B, C, D, solution);
    }
    if (sk_double_nearly_zero(D)) {  // 0 is one root
        int num = SkQuads::RootsReal(A, B, C, solution);
        for (int i = 0; i < num; ++i) {
            if (sk_double_nearly_zero(solution[i])) {
                return num;
            }
        }
        solution[num++] = 0;
        return num;
    }
    if (sk_double_nearly_zero(A + B + C + D)) {  // 1 is one root
        int num = SkQuads::RootsReal(A, A + B, -D, solution);
        for (int i = 0; i < num; ++i) {
            if (sk_doubles_nearly_equal_ulps(solution[i], 1)) {
                return num;
            }
        }
        solution[num++] = 1;
        return num;
    }
    double a, b, c;
    {
        // If A is zero (e.g. B was nan and thus close_to_a_quadratic was false), we will
        // temporarily have infinities rolling about, but will catch that when checking
        // R2MinusQ3.
        double invA = sk_ieee_double_divide(1, A);
        a = B * invA;
        b = C * invA;
        c = D * invA;
    }
    double a2 = a * a;
    double Q = (a2 - b * 3) / 9;
    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
    double R2 = R * R;
    double Q3 = Q * Q * Q;
    double R2MinusQ3 = R2 - Q3;
    // If one of R2 Q3 is infinite or nan, subtracting them will also be infinite/nan.
    // If both are infinite or nan, the subtraction will be nan.
    // In either case, we have no finite roots.
    if (!SkIsFinite(R2MinusQ3)) {
        return 0;
    }
    double adiv3 = a / 3;
    double r;
    double* roots = solution;
    if (R2MinusQ3 < 0) {   // we have 3 real roots
        // the divide/root can, due to finite precisions, be slightly outside of -1...1
        const double theta = acos(SkTPin(R / std::sqrt(Q3), -1., 1.));
        const double neg2RootQ = -2 * std::sqrt(Q);

        r = neg2RootQ * cos(theta / 3) - adiv3;
        *roots++ = r;

        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
        if (!nearly_equal(solution[0], r)) {
            *roots++ = r;
        }
        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
        if (!nearly_equal(solution[0], r) &&
            (roots - solution == 1 || !nearly_equal(solution[1], r))) {
            *roots++ = r;
        }
    } else {  // we have 1 real root
        const double sqrtR2MinusQ3 = std::sqrt(R2MinusQ3);
        A = fabs(R) + sqrtR2MinusQ3;
        A = std::cbrt(A); // cube root
        if (R > 0) {
            A = -A;
        }
        if (!sk_double_nearly_zero(A)) {
            A += Q / A;
        }
        r = A - adiv3;
        *roots++ = r;
        if (!sk_double_nearly_zero(R2) &&
             sk_doubles_nearly_equal_ulps(R2, Q3)) {
            r = -A / 2 - adiv3;
            if (!nearly_equal(solution[0], r)) {
                *roots++ = r;
            }
        }
    }
    return static_cast<int>(roots - solution);
}

int SkCubics::RootsValidT(double A, double B, double C, double D,
                          double solution[3]) {
    double allRoots[3] = {0, 0, 0};
    int realRoots = SkCubics::RootsReal(A, B, C, D, allRoots);
    int foundRoots = 0;
    for (int index = 0; index < realRoots; ++index) {
        double tValue = allRoots[index];
        if (tValue >= 1.0 && tValue <= 1.00005) {
            // Make sure we do not already have 1 (or something very close) in the list of roots.
            if ((foundRoots < 1 || !sk_doubles_nearly_equal_ulps(solution[0], 1)) &&
                (foundRoots < 2 || !sk_doubles_nearly_equal_ulps(solution[1], 1))) {
                solution[foundRoots++] = 1;
            }
        } else if (tValue >= -0.00005 && (tValue <= 0.0 || sk_double_nearly_zero(tValue))) {
            // Make sure we do not already have 0 (or something very close) in the list of roots.
            if ((foundRoots < 1 || !sk_double_nearly_zero(solution[0])) &&
                (foundRoots < 2 || !sk_double_nearly_zero(solution[1]))) {
                solution[foundRoots++] = 0;
            }
        } else if (tValue > 0.0 && tValue < 1.0) {
            solution[foundRoots++] = tValue;
        }
    }
    return foundRoots;
}

static bool approximately_zero(double x) {
    // This cutoff for our binary search hopefully strikes a good balance between
    // performance and accuracy.
    return std::abs(x) < 0.00000001;
}

static int find_extrema_valid_t(double A, double B, double C,
                                double t[2]) {
    // To find the local min and max of a cubic, we take the derivative and
    // solve when that is equal to 0.
    // d/dt (A*t^3 + B*t^2 + C*t + D) = 3A*t^2 + 2B*t + C
    double roots[2] = {0, 0};
    int numRoots = SkQuads::RootsReal(3*A, 2*B, C, roots);
    int validRoots = 0;
    for (int i = 0; i < numRoots; i++) {
        double tValue = roots[i];
        if (tValue >= 0 && tValue <= 1.0) {
            t[validRoots++] = tValue;
        }
    }
    return validRoots;
}

static double binary_search(double A, double B, double C, double D, double start, double stop) {
    SkASSERT(start <= stop);
    double left = SkCubics::EvalAt(A, B, C, D, start);
    if (approximately_zero(left)) {
        return start;
    }
    double right = SkCubics::EvalAt(A, B, C, D, stop);
    if (!SkIsFinite(left, right)) {
        return -1; // Not going to deal with one or more endpoints being non-finite.
    }
    if ((left > 0 && right > 0) || (left < 0 && right < 0)) {
        return -1; // We can only have a root if one is above 0 and the other is below 0.
    }

    constexpr int maxIterations = 1000; // prevent infinite loop
    for (int i = 0; i < maxIterations; i++) {
        double step = (start + stop) / 2;
        double curr = SkCubics::EvalAt(A, B, C, D, step);
        if (approximately_zero(curr)) {
            return step;
        }
        if ((curr < 0 && left < 0) || (curr > 0 && left > 0)) {
            // go right
            start = step;
        } else {
            // go left
            stop = step;
        }
    }
    return -1;
}

int SkCubics::BinarySearchRootsValidT(double A, double B, double C, double D,
                                      double solution[3]) {
    if (!SkIsFinite(A, B, C, D)) {
        return 0;
    }
    double regions[4] = {0, 0, 0, 1};
    // Find local minima and maxima
    double minMax[2] = {0, 0};
    int extremaCount = find_extrema_valid_t(A, B, C, minMax);
    int startIndex = 2 - extremaCount;
    if (extremaCount == 1) {
        regions[startIndex + 1] = minMax[0];
    }
    if (extremaCount == 2) {
        // While the roots will be in the range 0 to 1 inclusive, they might not be sorted.
        regions[startIndex + 1] = std::min(minMax[0], minMax[1]);
        regions[startIndex + 2] = std::max(minMax[0], minMax[1]);
    }
    // Starting at regions[startIndex] and going up through regions[3], we have
    // an ascending list of numbers in the range 0 to 1.0, between which are the possible
    // locations of a root.
    int foundRoots = 0;
    for (;startIndex < 3; startIndex++) {
        double root = binary_search(A, B, C, D, regions[startIndex], regions[startIndex + 1]);
        if (root >= 0) {
            // Check for duplicates
            if ((foundRoots < 1 || !approximately_zero(solution[0] - root)) &&
                (foundRoots < 2 || !approximately_zero(solution[1] - root))) {
                solution[foundRoots++] = root;
            }
        }
    }
    return foundRoots;
}