| /* |
| * 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]); |
| } |
| } |
