1 /*
2 * Copyright 2023 Google LLC
3 *
4 * Use of this source code is governed by a BSD-style license that can be
5 * found in the LICENSE file.
6 */
7
8 #include "include/core/SkSpan.h"
9 #include "include/core/SkTypes.h"
10 #include "include/private/base/SkFloatingPoint.h"
11 #include "src/base/SkCubics.h"
12 #include "src/base/SkUtils.h"
13 #include "src/pathops/SkPathOpsCubic.h"
14 #include "tests/Test.h"
15
16 #include <algorithm>
17 #include <cmath>
18 #include <cstddef>
19 #include <iterator>
20 #include <string>
21
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)22 static void testCubicRootsReal(skiatest::Reporter* reporter, std::string name,
23 double A, double B, double C, double D,
24 SkSpan<const double> expectedRoots,
25 bool skipPathops = false,
26 bool skipRootValidation = false) {
27 skiatest::ReporterContext subtest(reporter, name);
28 // Validate test case
29 REPORTER_ASSERT(reporter, expectedRoots.size() <= 3,
30 "Invalid test case, up to 3 roots allowed");
31
32 for (size_t i = 0; i < expectedRoots.size(); i++) {
33 double x = expectedRoots[i];
34 // A*x^3 + B*x^2 + C*x + D should equal 0 (unless floating point error causes issues)
35 double y = A * x * x * x + B * x * x + C * x + D;
36 if (!skipRootValidation) {
37 REPORTER_ASSERT(reporter, sk_double_nearly_zero(y),
38 "Invalid test case root %zu. %.16f != 0", i, y);
39 }
40
41 if (i > 0) {
42 REPORTER_ASSERT(reporter, expectedRoots[i-1] <= expectedRoots[i],
43 "Invalid test case root %zu. Roots should be sorted in ascending order", i);
44 }
45 }
46
47 // The old pathops implementation sometimes gives incorrect solutions. We can opt
48 // our tests out of checking that older implementation if that causes issues.
49 if (!skipPathops) {
50 skiatest::ReporterContext subsubtest(reporter, "Pathops Implementation");
51 double roots[3] = {0, 0, 0};
52 int rootCount = SkDCubic::RootsReal(A, B, C, D, roots);
53 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
54 "Wrong number of roots returned %zu != %d", expectedRoots.size(),
55 rootCount);
56
57 // We don't care which order the roots are returned from the algorithm.
58 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
59 std::sort(std::begin(roots), std::begin(roots) + rootCount);
60 for (int i = 0; i < rootCount; i++) {
61 if (sk_double_nearly_zero(expectedRoots[i])) {
62 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
63 "0 != %.16f at index %d", roots[i], i);
64 } else {
65 REPORTER_ASSERT(reporter,
66 sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
67 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
68 }
69 }
70 }
71 {
72 skiatest::ReporterContext subsubtest(reporter, "SkCubics Analytic Implementation");
73 double roots[3] = {0, 0, 0};
74 int rootCount = SkCubics::RootsReal(A, B, C, D, roots);
75 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
76 "Wrong number of roots returned %zu != %d", expectedRoots.size(),
77 rootCount);
78
79 // We don't care which order the roots are returned from the algorithm.
80 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
81 std::sort(std::begin(roots), std::begin(roots) + rootCount);
82 for (int i = 0; i < rootCount; i++) {
83 if (sk_double_nearly_zero(expectedRoots[i])) {
84 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
85 "0 != %.16f at index %d", roots[i], i);
86 } else {
87 REPORTER_ASSERT(reporter,
88 sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
89 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
90 }
91 }
92 }
93 }
94
DEF_TEST(CubicRootsReal_ActualCubics,reporter)95 DEF_TEST(CubicRootsReal_ActualCubics, reporter) {
96 // All answers are given with 16 significant digits (max for a double) or as an integer
97 // when the answer is exact.
98 testCubicRootsReal(reporter, "one root 1x^3 + 2x^2 + 3x + 4",
99 1, 2, 3, 4,
100 {-1.650629191439388});
101 //-1.650629191439388218880801 from Wolfram Alpha
102
103 // (3x-5)(6x-10)(x+4) = 18x^3 + 12x^2 - 190x + 200
104 testCubicRootsReal(reporter, "touches y axis 18x^3 + 12x^2 - 190x + 200",
105 18, 12, -190, 200,
106 {-4.,
107 1.666666666666667, // 5/3
108 });
109
110 testCubicRootsReal(reporter, "three roots 10x^3 - 20x^2 - 30x + 40",
111 10, -20, -30, 40,
112 {-1.561552812808830,
113 //-1.561552812808830274910705 from Wolfram Alpha
114 1.,
115 2.561552812808830,
116 // 2.561552812808830274910705 from Wolfram Alpha
117 });
118
119 testCubicRootsReal(reporter, "three roots -10x^3 + 200x^2 + 300x - 400",
120 -10, 200, 300, -400,
121 {-2.179884793243323,
122 //-2.179884793243323422232630 from Wolfram Alpha
123 0.8607083693981839,
124 // 0.8607083693981838897123320 from Wolfram Alpha
125 21.31917642384514,
126 //21.31917642384513953252030 from Wolfram Alpha
127 });
128
129 testCubicRootsReal(reporter, "one root -x^3 + 0x^2 + 5x - 7",
130 -1, 0, 5, -7,
131 {-2.747346540307211,
132 //-2.747346540307210849971436 from Wolfram Alpha
133 });
134
135 testCubicRootsReal(reporter, "one root 2x^3 - 3x^2 + 0x + 3",
136 2, -3, 0, 3,
137 {-0.806443932358772,
138 //-0.8064439323587723772036250 from Wolfram Alpha
139 });
140
141 testCubicRootsReal(reporter, "one root x^3 + 0x^2 + 0x - 9",
142 1, 0, 0, -9,
143 {2.080083823051904,
144 //2.0800838230519041145300568 from Wolfram Alpha
145 });
146
147 testCubicRootsReal(reporter, "three roots 2x^3 - 3x^2 - 4x + 0",
148 2, -3, -4, 0,
149 {-0.8507810593582122,
150 //-0.8507810593582121716220544 from Wolfram Alpha
151 0.,
152 2.350781059358212
153 //2.350781059358212171622054 from Wolfram Alpha
154 });
155
156 testCubicRootsReal(reporter, "R^2 and Q^3 are near zero",
157 -0.33790159225463867, -0.81997990608215332,
158 -0.66327774524688721, -0.17884063720703125,
159 {-0.7995944894729731});
160
161 // The following three cases fallback to treating the cubic as a quadratic.
162 // Otherwise, floating point error mangles the solutions near +- 1
163 // This means we don't find all the roots, but usually we only care about roots
164 // in the range [0, 1], so that is ok.
165 testCubicRootsReal(reporter, "oss-fuzz:55625 Two roots near zero, one big root",
166 sk_bit_cast<double>(0xbf1a8de580000000), // -0.00010129655
167 sk_bit_cast<double>(0x4106c0c680000000), // 186392.8125
168 0.0,
169 sk_bit_cast<double>(0xc104c0ce80000000), // -170009.8125
170 { -0.9550418733785169, // Wolfram Alpha puts the root at X = 0.955042
171 0.9550418733785169, // (~2e7 error)
172 // 1.84007e9 is the other root, which we do not find.
173 },
174 true /* == skipPathops */, true /* == skipRootValidation */);
175
176 testCubicRootsReal(reporter, "oss-fuzz:55625 Two roots near zero, one big root, near linear",
177 sk_bit_cast<double>(0x3c04040400000000), // -1.3563156-19
178 sk_bit_cast<double>(0x4106c0c680000000), // 186392.8125
179 0.0,
180 sk_bit_cast<double>(0xc104c0ce80000000), // -170009.8125
181 { -0.9550418733785169,
182 0.9550418733785169,
183 // 1.84007e9 is the other root, which we do not find.
184 },
185 true /* == skipPathops */);
186
187 testCubicRootsReal(reporter, "oss-fuzz:55680 A nearly zero, C is zero",
188 sk_bit_cast<double>(0x3eb0000000000000), // 9.5367431640625000e-07
189 sk_bit_cast<double>(0x409278a560000000), // 1182.1614990234375
190 0.0,
191 sk_bit_cast<double>(0xc092706160000000), // -1180.0950927734375
192 { -0.9991256228290017,
193 // -0.9991256232316570469050229 according to Wolfram Alpha (~1e-09 error)
194 0.9991256228290017,
195 // 0.9991256224263463476403026 according to Wolfram Alpha (~1e-09 error)
196 // 1.239586176×10^9 is the other root, which we do not find.
197 },
198 true, true /* == skipRootValidation */);
199 }
200
DEF_TEST(CubicRootsReal_Quadratics,reporter)201 DEF_TEST(CubicRootsReal_Quadratics, reporter) {
202 testCubicRootsReal(reporter, "two roots -2x^2 + 3x + 4",
203 0, -2, 3, 4,
204 {-0.8507810593582122,
205 //-0.8507810593582121716220544 from Wolfram Alpha
206 2.350781059358212,
207 // 2.350781059358212171622054 from Wolfram Alpha
208 });
209
210 testCubicRootsReal(reporter, "touches y axis -x^2 + 3x + 4",
211 0, -2, 3, 4,
212 {-0.8507810593582122,
213 //-0.8507810593582121716220544 from Wolfram Alpha
214 2.350781059358212,
215 // 2.350781059358212171622054 from Wolfram Alpha
216 });
217
218 testCubicRootsReal(reporter, "no roots x^2 + 2x + 7",
219 0, 1, 2, 7,
220 {});
221
222 // similar to oss-fuzz:55680
223 testCubicRootsReal(reporter, "two roots one small one big (and ignored)",
224 0, -0.01, 200000000000000, -120000000000000,
225 { 0.6 },
226 true /* == skipPathops */);
227 }
228
DEF_TEST(CubicRootsReal_Linear,reporter)229 DEF_TEST(CubicRootsReal_Linear, reporter) {
230 testCubicRootsReal(reporter, "positive slope 3x + 4",
231 0, 0, 3, 4,
232 {-1.333333333333333});
233
234 testCubicRootsReal(reporter, "negative slope -2x - 8",
235 0, 0, -2, -8,
236 {-4.});
237 }
238
DEF_TEST(CubicRootsReal_Constant,reporter)239 DEF_TEST(CubicRootsReal_Constant, reporter) {
240 testCubicRootsReal(reporter, "No intersections y = 4",
241 0, 0, 0, 4,
242 {});
243
244 testCubicRootsReal(reporter, "Infinite solutions y = 0",
245 0, 0, 0, 0,
246 {0.});
247 }
248
DEF_TEST(CubicRootsReal_NonFiniteNumbers,reporter)249 DEF_TEST(CubicRootsReal_NonFiniteNumbers, reporter) {
250 // The Pathops implementation does not check for infinities nor nans in all cases.
251 double roots[3] = {0, 0, 0};
252 REPORTER_ASSERT(reporter,
253 SkCubics::RootsReal(NAN, 1, 2, 3, roots) == 0,
254 "Nan A"
255 );
256 REPORTER_ASSERT(reporter,
257 SkCubics::RootsReal(1, NAN, 2, 3, roots) == 0,
258 "Nan B"
259 );
260 REPORTER_ASSERT(reporter,
261 SkCubics::RootsReal(1, 2, NAN, 3, roots) == 0,
262 "Nan C"
263 );
264 REPORTER_ASSERT(reporter,
265 SkCubics::RootsReal(1, 2, 3, NAN, roots) == 0,
266 "Nan D"
267 );
268
269 {
270 skiatest::ReporterContext subtest(reporter, "oss-fuzz:55419 C and D are large");
271 int numRoots = SkCubics::RootsReal(
272 -2, 0,
273 sk_bit_cast<double>(0xd5422020202020ff), //-5.074559e+102
274 sk_bit_cast<double>(0x600fff202020ff20), // 5.362551e+154
275 roots);
276 REPORTER_ASSERT(reporter, numRoots == 0, "No finite roots expected, got %d", numRoots);
277 }
278 {
279 skiatest::ReporterContext subtest(reporter, "oss-fuzz:55829 A is zero and B is NAN");
280 int numRoots = SkCubics::RootsReal(
281 0,
282 sk_bit_cast<double>(0xffffffffffff2020), //-nan
283 sk_bit_cast<double>(0x20202020202020ff), // 6.013470e-154
284 sk_bit_cast<double>(0xff20202020202020), //-2.211661e+304
285 roots);
286 REPORTER_ASSERT(reporter, numRoots == 0, "No finite roots expected, got %d", numRoots);
287 }
288 }
289
testCubicValidT(skiatest::Reporter * reporter,std::string name,double A,double B,double C,double D,SkSpan<const double> expectedRoots)290 static void testCubicValidT(skiatest::Reporter* reporter, std::string name,
291 double A, double B, double C, double D,
292 SkSpan<const double> expectedRoots) {
293 skiatest::ReporterContext subtest(reporter, name);
294 // Validate test case
295 REPORTER_ASSERT(reporter, expectedRoots.size() <= 3,
296 "Invalid test case, up to 3 roots allowed");
297
298 for (size_t i = 0; i < expectedRoots.size(); i++) {
299 double x = expectedRoots[i];
300 REPORTER_ASSERT(reporter, x >= 0 && x <= 1,
301 "Invalid test case root %zu. Roots must be in [0, 1]", i);
302
303 // A*x^3 + B*x^2 + C*x + D should equal 0
304 double y = A * x * x * x + B * x * x + C * x + D;
305 REPORTER_ASSERT(reporter, sk_double_nearly_zero(y),
306 "Invalid test case root %zu. %.16f != 0", i, y);
307
308 if (i > 0) {
309 REPORTER_ASSERT(reporter, expectedRoots[i-1] <= expectedRoots[i],
310 "Invalid test case root %zu. Roots should be sorted in ascending order", i);
311 }
312 }
313
314 {
315 skiatest::ReporterContext subsubtest(reporter, "Pathops Implementation");
316 double roots[3] = {0, 0, 0};
317 int rootCount = SkDCubic::RootsValidT(A, B, C, D, roots);
318 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
319 "Wrong number of roots returned %zu != %d",
320 expectedRoots.size(), rootCount);
321
322 // We don't care which order the roots are returned from the algorithm.
323 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
324 std::sort(std::begin(roots), std::begin(roots) + rootCount);
325 for (int i = 0; i < rootCount; i++) {
326 if (sk_double_nearly_zero(expectedRoots[i])) {
327 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
328 "0 != %.16f at index %d", roots[i], i);
329 } else {
330 REPORTER_ASSERT(reporter,
331 sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
332 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
333 }
334 }
335 }
336 {
337 skiatest::ReporterContext subsubtest(reporter, "SkCubics Analytic Implementation");
338 double roots[3] = {0, 0, 0};
339 int rootCount = SkCubics::RootsValidT(A, B, C, D, roots);
340 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
341 "Wrong number of roots returned %zu != %d",
342 expectedRoots.size(), rootCount);
343
344 // We don't care which order the roots are returned from the algorithm.
345 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
346 std::sort(std::begin(roots), std::begin(roots) + rootCount);
347 for (int i = 0; i < rootCount; i++) {
348 if (sk_double_nearly_zero(expectedRoots[i])) {
349 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
350 "0 != %.16f at index %d", roots[i], i);
351 } else {
352 REPORTER_ASSERT(reporter,
353 sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
354 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
355 }
356 }
357 }
358 {
359 skiatest::ReporterContext subsubtest(reporter, "SkCubics Binary Search Implementation");
360 double roots[3] = {0, 0, 0};
361 int rootCount = SkCubics::BinarySearchRootsValidT(A, B, C, D, roots);
362 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
363 "Wrong number of roots returned %zu != %d", expectedRoots.size(),
364 rootCount);
365
366 // We don't care which order the roots are returned from the algorithm.
367 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
368 std::sort(std::begin(roots), std::begin(roots) + rootCount);
369 for (int i = 0; i < rootCount; i++) {
370 double delta = std::abs(roots[i] - expectedRoots[i]);
371 REPORTER_ASSERT(reporter,
372 // Binary search is not absolutely accurate all the time, but
373 // it should be accurate enough reliably
374 delta < 0.000001,
375 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
376 }
377 }
378 }
379
DEF_TEST(CubicRootsValidT,reporter)380 DEF_TEST(CubicRootsValidT, reporter) {
381 // All answers are given with 16 significant digits (max for a double) or as an integer
382 // when the answer is exact.
383 testCubicValidT(reporter, "three roots 24x^3 - 46x^2 + 29x - 6",
384 24, -46, 29, -6,
385 {0.5,
386 0.6666666666666667,
387 0.75});
388
389 testCubicValidT(reporter, "three roots total, two in range 54x^3 - 117x^2 + 45x + 0",
390 54, -117, 45, 0,
391 {0.0,
392 0.5,
393 // 5/3 is the other root, but not in [0, 1]
394 });
395
396 testCubicValidT(reporter, "one root = 1 10x^3 - 20x^2 - 30x + 40",
397 10, -20, -30, 40,
398 {1.0});
399
400 testCubicValidT(reporter, "one root = 0 2x^3 - 3x^2 - 4x + 0",
401 2, -3, -4, 0,
402 {0.0});
403
404 testCubicValidT(reporter, "three roots total, two in range -2x^3 - 3x^2 + 4x + 0",
405 -2, -3, 4, 0,
406 { 0.0,
407 0.8507810593582122,
408 // 0.8507810593582121716220544 from Wolfram Alpha
409 });
410
411 // x(x-1) = x^2 - x
412 testCubicValidT(reporter, "Two roots at exactly 0 and 1",
413 0, 1, -1, 0,
414 {0.0, 1.0});
415
416 testCubicValidT(reporter, "Single point has one root",
417 0, 0, 0, 0,
418 {0.0});
419 }
420
DEF_TEST(CubicRootsValidT_ClampToZeroAndOne,reporter)421 DEF_TEST(CubicRootsValidT_ClampToZeroAndOne, reporter) {
422 {
423 // (x + 0.00001)(x - 1.00005), but the answers will be 0 and 1
424 double A = 0.;
425 double B = 1.;
426 double C = -1.00004;
427 double D = -0.0000100005;
428 double roots[3] = {0, 0, 0};
429 int rootCount = SkDCubic::RootsValidT(A, B, C, D, roots);
430
431 REPORTER_ASSERT(reporter, rootCount == 2);
432 std::sort(std::begin(roots), std::begin(roots) + rootCount);
433 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[0]), "%.16f != 0", roots[0]);
434 REPORTER_ASSERT(reporter, sk_doubles_nearly_equal_ulps(roots[1], 1), "%.16f != 1", roots[1]);
435 }
436 {
437 // Three very small roots, all of them are nearly equal zero
438 // (1 - 10000000000x)(1 - 20000000000x)(1 - 30000000000x)
439 // -6000000000000000000000000000000 x^3 + 1100000000000000000000 x^2 - 60000000000 x + 1
440 double A = -6.0e30;
441 double B = 1.1e21;
442 double C = -6.0e10;
443 double D = 1;
444 double roots[3] = {0, 0, 0};
445 int rootCount = SkDCubic::RootsValidT(A, B, C, D, roots);
446
447 REPORTER_ASSERT(reporter, rootCount == 1);
448 std::sort(std::begin(roots), std::begin(roots) + rootCount);
449 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[0]), "%.16f != 0", roots[0]);
450 }
451 }
452