blob: 80bd772b4b92e7bca3972cdae3ce5d5b0e9b5f7b [file] [log] [blame]
/*
* 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 "include/core/SkSpan.h"
#include "include/core/SkTypes.h"
#include "include/private/base/SkFloatingPoint.h"
#include "src/base/SkCubics.h"
#include "src/base/SkUtils.h"
#include "src/pathops/SkPathOpsCubic.h"
#include "tests/Test.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iterator>
#include <string>
static void testCubicRootsReal(skiatest::Reporter* reporter, std::string name,
double A, double B, double C, double D,
SkSpan<const double> expectedRoots,
bool skipPathops = false,
bool skipRootValidation = false) {
skiatest::ReporterContext subtest(reporter, name);
// Validate test case
REPORTER_ASSERT(reporter, expectedRoots.size() <= 3,
"Invalid test case, up to 3 roots allowed");
for (size_t i = 0; i < expectedRoots.size(); i++) {
double x = expectedRoots[i];
// A*x^3 + B*x^2 + C*x + D should equal 0 (unless floating point error causes issues)
double y = A * x * x * x + B * x * x + C * x + D;
if (!skipRootValidation) {
REPORTER_ASSERT(reporter, sk_double_nearly_zero(y),
"Invalid test case root %zu. %.16f != 0", i, y);
}
if (i > 0) {
REPORTER_ASSERT(reporter, expectedRoots[i-1] <= expectedRoots[i],
"Invalid test case root %zu. Roots should be sorted in ascending order", i);
}
}
// The old pathops implementation sometimes gives incorrect solutions. We can opt
// our tests out of checking that older implementation if that causes issues.
if (!skipPathops) {
skiatest::ReporterContext subsubtest(reporter, "Pathops Implementation");
double roots[3] = {0, 0, 0};
int rootCount = SkDCubic::RootsReal(A, B, C, D, roots);
REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
"Wrong number of roots returned %zu != %d", expectedRoots.size(),
rootCount);
// We don't care which order the roots are returned from the algorithm.
// For determinism, we will sort them (and ensure the provided solutions are also sorted).
std::sort(std::begin(roots), std::begin(roots) + rootCount);
for (int i = 0; i < rootCount; i++) {
if (sk_double_nearly_zero(expectedRoots[i])) {
REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
"0 != %.16f at index %d", roots[i], i);
} else {
REPORTER_ASSERT(reporter,
sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
"%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
}
}
}
{
skiatest::ReporterContext subsubtest(reporter, "SkCubics Analytic Implementation");
double roots[3] = {0, 0, 0};
int rootCount = SkCubics::RootsReal(A, B, C, D, roots);
REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
"Wrong number of roots returned %zu != %d", expectedRoots.size(),
rootCount);
// We don't care which order the roots are returned from the algorithm.
// For determinism, we will sort them (and ensure the provided solutions are also sorted).
std::sort(std::begin(roots), std::begin(roots) + rootCount);
for (int i = 0; i < rootCount; i++) {
if (sk_double_nearly_zero(expectedRoots[i])) {
REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
"0 != %.16f at index %d", roots[i], i);
} else {
REPORTER_ASSERT(reporter,
sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
"%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
}
}
}
}
DEF_TEST(CubicRootsReal_ActualCubics, reporter) {
// All answers are given with 16 significant digits (max for a double) or as an integer
// when the answer is exact.
testCubicRootsReal(reporter, "one root 1x^3 + 2x^2 + 3x + 4",
1, 2, 3, 4,
{-1.650629191439388});
//-1.650629191439388218880801 from Wolfram Alpha
// (3x-5)(6x-10)(x+4) = 18x^3 + 12x^2 - 190x + 200
testCubicRootsReal(reporter, "touches y axis 18x^3 + 12x^2 - 190x + 200",
18, 12, -190, 200,
{-4.,
1.666666666666667, // 5/3
});
testCubicRootsReal(reporter, "three roots 10x^3 - 20x^2 - 30x + 40",
10, -20, -30, 40,
{-1.561552812808830,
//-1.561552812808830274910705 from Wolfram Alpha
1.,
2.561552812808830,
// 2.561552812808830274910705 from Wolfram Alpha
});
testCubicRootsReal(reporter, "three roots -10x^3 + 200x^2 + 300x - 400",
-10, 200, 300, -400,
{-2.179884793243323,
//-2.179884793243323422232630 from Wolfram Alpha
0.8607083693981839,
// 0.8607083693981838897123320 from Wolfram Alpha
21.31917642384514,
//21.31917642384513953252030 from Wolfram Alpha
});
testCubicRootsReal(reporter, "one root -x^3 + 0x^2 + 5x - 7",
-1, 0, 5, -7,
{-2.747346540307211,
//-2.747346540307210849971436 from Wolfram Alpha
});
testCubicRootsReal(reporter, "one root 2x^3 - 3x^2 + 0x + 3",
2, -3, 0, 3,
{-0.806443932358772,
//-0.8064439323587723772036250 from Wolfram Alpha
});
testCubicRootsReal(reporter, "one root x^3 + 0x^2 + 0x - 9",
1, 0, 0, -9,
{2.080083823051904,
//2.0800838230519041145300568 from Wolfram Alpha
});
testCubicRootsReal(reporter, "three roots 2x^3 - 3x^2 - 4x + 0",
2, -3, -4, 0,
{-0.8507810593582122,
//-0.8507810593582121716220544 from Wolfram Alpha
0.,
2.350781059358212
//2.350781059358212171622054 from Wolfram Alpha
});
testCubicRootsReal(reporter, "R^2 and Q^3 are near zero",
-0.33790159225463867, -0.81997990608215332,
-0.66327774524688721, -0.17884063720703125,
{-0.7995944894729731});
// The following three cases fallback to treating the cubic as a quadratic.
// Otherwise, floating point error mangles the solutions near +- 1
// This means we don't find all the roots, but usually we only care about roots
// in the range [0, 1], so that is ok.
testCubicRootsReal(reporter, "oss-fuzz:55625 Two roots near zero, one big root",
sk_bit_cast<double>(0xbf1a8de580000000), // -0.00010129655
sk_bit_cast<double>(0x4106c0c680000000), // 186392.8125
0.0,
sk_bit_cast<double>(0xc104c0ce80000000), // -170009.8125
{ -0.9550418733785169, // Wolfram Alpha puts the root at X = 0.955042
0.9550418733785169, // (~2e7 error)
// 1.84007e9 is the other root, which we do not find.
},
true /* == skipPathops */, true /* == skipRootValidation */);
testCubicRootsReal(reporter, "oss-fuzz:55625 Two roots near zero, one big root, near linear",
sk_bit_cast<double>(0x3c04040400000000), // -1.3563156-19
sk_bit_cast<double>(0x4106c0c680000000), // 186392.8125
0.0,
sk_bit_cast<double>(0xc104c0ce80000000), // -170009.8125
{ -0.9550418733785169,
0.9550418733785169,
// 1.84007e9 is the other root, which we do not find.
},
true /* == skipPathops */);
testCubicRootsReal(reporter, "oss-fuzz:55680 A nearly zero, C is zero",
sk_bit_cast<double>(0x3eb0000000000000), // 9.5367431640625000e-07
sk_bit_cast<double>(0x409278a560000000), // 1182.1614990234375
0.0,
sk_bit_cast<double>(0xc092706160000000), // -1180.0950927734375
{ -0.9991256228290017,
// -0.9991256232316570469050229 according to Wolfram Alpha (~1e-09 error)
0.9991256228290017,
// 0.9991256224263463476403026 according to Wolfram Alpha (~1e-09 error)
// 1.239586176×10^9 is the other root, which we do not find.
},
true, true /* == skipRootValidation */);
}
DEF_TEST(CubicRootsReal_Quadratics, reporter) {
testCubicRootsReal(reporter, "two roots -2x^2 + 3x + 4",
0, -2, 3, 4,
{-0.8507810593582122,
//-0.8507810593582121716220544 from Wolfram Alpha
2.350781059358212,
// 2.350781059358212171622054 from Wolfram Alpha
});
testCubicRootsReal(reporter, "touches y axis -x^2 + 3x + 4",
0, -2, 3, 4,
{-0.8507810593582122,
//-0.8507810593582121716220544 from Wolfram Alpha
2.350781059358212,
// 2.350781059358212171622054 from Wolfram Alpha
});
testCubicRootsReal(reporter, "no roots x^2 + 2x + 7",
0, 1, 2, 7,
{});
// similar to oss-fuzz:55680
testCubicRootsReal(reporter, "two roots one small one big (and ignored)",
0, -0.01, 200000000000000, -120000000000000,
{ 0.6 },
true /* == skipPathops */);
}
DEF_TEST(CubicRootsReal_Linear, reporter) {
testCubicRootsReal(reporter, "positive slope 3x + 4",
0, 0, 3, 4,
{-1.333333333333333});
testCubicRootsReal(reporter, "negative slope -2x - 8",
0, 0, -2, -8,
{-4.});
}
DEF_TEST(CubicRootsReal_Constant, reporter) {
testCubicRootsReal(reporter, "No intersections y = 4",
0, 0, 0, 4,
{});
testCubicRootsReal(reporter, "Infinite solutions y = 0",
0, 0, 0, 0,
{0.});
}
DEF_TEST(CubicRootsReal_NonFiniteNumbers, reporter) {
// The Pathops implementation does not check for infinities nor nans in all cases.
double roots[3] = {0, 0, 0};
REPORTER_ASSERT(reporter,
SkCubics::RootsReal(NAN, 1, 2, 3, roots) == 0,
"Nan A"
);
REPORTER_ASSERT(reporter,
SkCubics::RootsReal(1, NAN, 2, 3, roots) == 0,
"Nan B"
);
REPORTER_ASSERT(reporter,
SkCubics::RootsReal(1, 2, NAN, 3, roots) == 0,
"Nan C"
);
REPORTER_ASSERT(reporter,
SkCubics::RootsReal(1, 2, 3, NAN, roots) == 0,
"Nan D"
);
{
skiatest::ReporterContext subtest(reporter, "oss-fuzz:55419 C and D are large");
int numRoots = SkCubics::RootsReal(
-2, 0,
sk_bit_cast<double>(0xd5422020202020ff), //-5.074559e+102
sk_bit_cast<double>(0x600fff202020ff20), // 5.362551e+154
roots);
REPORTER_ASSERT(reporter, numRoots == 0, "No finite roots expected, got %d", numRoots);
}
{
skiatest::ReporterContext subtest(reporter, "oss-fuzz:55829 A is zero and B is NAN");
int numRoots = SkCubics::RootsReal(
0,
sk_bit_cast<double>(0xffffffffffff2020), //-nan
sk_bit_cast<double>(0x20202020202020ff), // 6.013470e-154
sk_bit_cast<double>(0xff20202020202020), //-2.211661e+304
roots);
REPORTER_ASSERT(reporter, numRoots == 0, "No finite roots expected, got %d", numRoots);
}
}
static void testCubicValidT(skiatest::Reporter* reporter, std::string name,
double A, double B, double C, double D,
SkSpan<const double> expectedRoots) {
skiatest::ReporterContext subtest(reporter, name);
// Validate test case
REPORTER_ASSERT(reporter, expectedRoots.size() <= 3,
"Invalid test case, up to 3 roots allowed");
for (size_t i = 0; i < expectedRoots.size(); i++) {
double x = expectedRoots[i];
REPORTER_ASSERT(reporter, x >= 0 && x <= 1,
"Invalid test case root %zu. Roots must be in [0, 1]", i);
// A*x^3 + B*x^2 + C*x + D should equal 0
double y = A * x * x * x + B * x * x + C * x + D;
REPORTER_ASSERT(reporter, sk_double_nearly_zero(y),
"Invalid test case root %zu. %.16f != 0", i, y);
if (i > 0) {
REPORTER_ASSERT(reporter, expectedRoots[i-1] <= expectedRoots[i],
"Invalid test case root %zu. Roots should be sorted in ascending order", i);
}
}
{
skiatest::ReporterContext subsubtest(reporter, "Pathops Implementation");
double roots[3] = {0, 0, 0};
int rootCount = SkDCubic::RootsValidT(A, B, C, D, roots);
REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
"Wrong number of roots returned %zu != %d",
expectedRoots.size(), rootCount);
// We don't care which order the roots are returned from the algorithm.
// For determinism, we will sort them (and ensure the provided solutions are also sorted).
std::sort(std::begin(roots), std::begin(roots) + rootCount);
for (int i = 0; i < rootCount; i++) {
if (sk_double_nearly_zero(expectedRoots[i])) {
REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
"0 != %.16f at index %d", roots[i], i);
} else {
REPORTER_ASSERT(reporter,
sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
"%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
}
}
}
{
skiatest::ReporterContext subsubtest(reporter, "SkCubics Analytic Implementation");
double roots[3] = {0, 0, 0};
int rootCount = SkCubics::RootsValidT(A, B, C, D, roots);
REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
"Wrong number of roots returned %zu != %d",
expectedRoots.size(), rootCount);
// We don't care which order the roots are returned from the algorithm.
// For determinism, we will sort them (and ensure the provided solutions are also sorted).
std::sort(std::begin(roots), std::begin(roots) + rootCount);
for (int i = 0; i < rootCount; i++) {
if (sk_double_nearly_zero(expectedRoots[i])) {
REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
"0 != %.16f at index %d", roots[i], i);
} else {
REPORTER_ASSERT(reporter,
sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
"%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
}
}
}
{
skiatest::ReporterContext subsubtest(reporter, "SkCubics Binary Search Implementation");
double roots[3] = {0, 0, 0};
int rootCount = SkCubics::BinarySearchRootsValidT(A, B, C, D, roots);
REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
"Wrong number of roots returned %zu != %d", expectedRoots.size(),
rootCount);
// We don't care which order the roots are returned from the algorithm.
// For determinism, we will sort them (and ensure the provided solutions are also sorted).
std::sort(std::begin(roots), std::begin(roots) + rootCount);
for (int i = 0; i < rootCount; i++) {
double delta = std::abs(roots[i] - expectedRoots[i]);
REPORTER_ASSERT(reporter,
// Binary search is not absolutely accurate all the time, but
// it should be accurate enough reliably
delta < 0.000001,
"%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
}
}
}
DEF_TEST(CubicRootsValidT, reporter) {
// All answers are given with 16 significant digits (max for a double) or as an integer
// when the answer is exact.
testCubicValidT(reporter, "three roots 24x^3 - 46x^2 + 29x - 6",
24, -46, 29, -6,
{0.5,
0.6666666666666667,
0.75});
testCubicValidT(reporter, "three roots total, two in range 54x^3 - 117x^2 + 45x + 0",
54, -117, 45, 0,
{0.0,
0.5,
// 5/3 is the other root, but not in [0, 1]
});
testCubicValidT(reporter, "one root = 1 10x^3 - 20x^2 - 30x + 40",
10, -20, -30, 40,
{1.0});
testCubicValidT(reporter, "one root = 0 2x^3 - 3x^2 - 4x + 0",
2, -3, -4, 0,
{0.0});
testCubicValidT(reporter, "three roots total, two in range -2x^3 - 3x^2 + 4x + 0",
-2, -3, 4, 0,
{ 0.0,
0.8507810593582122,
// 0.8507810593582121716220544 from Wolfram Alpha
});
// x(x-1) = x^2 - x
testCubicValidT(reporter, "Two roots at exactly 0 and 1",
0, 1, -1, 0,
{0.0, 1.0});
testCubicValidT(reporter, "Single point has one root",
0, 0, 0, 0,
{0.0});
}
DEF_TEST(CubicRootsValidT_ClampToZeroAndOne, reporter) {
{
// (x + 0.00001)(x - 1.00005), but the answers will be 0 and 1
double A = 0.;
double B = 1.;
double C = -1.00004;
double D = -0.0000100005;
double roots[3] = {0, 0, 0};
int rootCount = SkDCubic::RootsValidT(A, B, C, D, roots);
REPORTER_ASSERT(reporter, rootCount == 2);
std::sort(std::begin(roots), std::begin(roots) + rootCount);
REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[0]), "%.16f != 0", roots[0]);
REPORTER_ASSERT(reporter, sk_doubles_nearly_equal_ulps(roots[1], 1), "%.16f != 1", roots[1]);
}
{
// Three very small roots, all of them are nearly equal zero
// (1 - 10000000000x)(1 - 20000000000x)(1 - 30000000000x)
// -6000000000000000000000000000000 x^3 + 1100000000000000000000 x^2 - 60000000000 x + 1
double A = -6.0e30;
double B = 1.1e21;
double C = -6.0e10;
double D = 1;
double roots[3] = {0, 0, 0};
int rootCount = SkDCubic::RootsValidT(A, B, C, D, roots);
REPORTER_ASSERT(reporter, rootCount == 1);
std::sort(std::begin(roots), std::begin(roots) + rootCount);
REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[0]), "%.16f != 0", roots[0]);
}
}