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 "src/base/SkQuads.h"
9
10 #include "include/core/SkSpan.h"
11 #include "include/core/SkTypes.h"
12 #include "include/private/base/SkFloatingPoint.h"
13 #include "src/pathops/SkPathOpsQuad.h"
14 #include "tests/Test.h"
15
16 #include <algorithm>
17 #include <cfloat>
18 #include <cmath>
19 #include <cstddef>
20 #include <cstdint>
21 #include <iterator>
22 #include <limits>
23 #include <string>
24
testQuadRootsReal(skiatest::Reporter * reporter,const std::string & name,double A,double B,double C,SkSpan<const double> expectedRoots)25 static void testQuadRootsReal(skiatest::Reporter* reporter, const std::string& name,
26 double A, double B, double C,
27 SkSpan<const double> expectedRoots) {
28 skiatest::ReporterContext subtest(reporter, name);
29 // Validate test case
30 REPORTER_ASSERT(reporter, expectedRoots.size() <= 2,
31 "Invalid test case, up to 2 roots allowed");
32
33 for (size_t i = 0; i < expectedRoots.size(); i++) {
34 double x = expectedRoots[i];
35 // A*x^2 + B*x + C should equal 0
36 double y = A * x * x + B * x + C;
37 REPORTER_ASSERT(reporter, sk_double_nearly_zero(y),
38 "Invalid test case root %zu. %.16f != 0", i, y);
39
40 if (i > 0) {
41 REPORTER_ASSERT(reporter, expectedRoots[i-1] <= expectedRoots[i],
42 "Invalid test case root %zu. Roots should be sorted in ascending order", i);
43 }
44 }
45
46 {
47 skiatest::ReporterContext subsubtest(reporter, "Pathops Implementation");
48 double roots[2] = {0, 0};
49 int rootCount = SkDQuad::RootsReal(A, B, C, roots);
50 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
51 "Wrong number of roots returned %zu != %d", expectedRoots.size(),
52 rootCount);
53
54 // We don't care which order the roots are returned from the algorithm.
55 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
56 std::sort(std::begin(roots), std::begin(roots) + rootCount);
57 for (int i = 0; i < rootCount; i++) {
58 if (sk_double_nearly_zero(expectedRoots[i])) {
59 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
60 "0 != %.16f at index %d", roots[i], i);
61 } else {
62 REPORTER_ASSERT(reporter,
63 sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
64 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
65 }
66 }
67 }
68 {
69 skiatest::ReporterContext subsubtest(reporter, "SkQuads Implementation");
70 double roots[2] = {0, 0};
71 int rootCount = SkQuads::RootsReal(A, B, C, roots);
72 REPORTER_ASSERT(reporter, expectedRoots.size() == size_t(rootCount),
73 "Wrong number of roots returned %zu != %d", expectedRoots.size(),
74 rootCount);
75
76 // We don't care which order the roots are returned from the algorithm.
77 // For determinism, we will sort them (and ensure the provided solutions are also sorted).
78 std::sort(std::begin(roots), std::begin(roots) + rootCount);
79 for (int i = 0; i < rootCount; i++) {
80 if (sk_double_nearly_zero(expectedRoots[i])) {
81 REPORTER_ASSERT(reporter, sk_double_nearly_zero(roots[i]),
82 "0 != %.16f at index %d", roots[i], i);
83 } else {
84 REPORTER_ASSERT(reporter,
85 sk_doubles_nearly_equal_ulps(expectedRoots[i], roots[i], 64),
86 "%.16f != %.16f at index %d", expectedRoots[i], roots[i], i);
87 }
88 }
89 }
90 }
91
DEF_TEST(QuadRootsReal_ActualQuadratics,reporter)92 DEF_TEST(QuadRootsReal_ActualQuadratics, reporter) {
93 // All answers are given with 16 significant digits (max for a double) or as an integer
94 // when the answer is exact.
95 testQuadRootsReal(reporter, "two roots 3x^2 - 20x - 40",
96 3, -20, -40,
97 {-1.610798991397109,
98 //-1.610798991397108632474265 from Wolfram Alpha
99 8.277465658063775,
100 // 8.277465658063775299140932 from Wolfram Alpha
101 });
102
103 // (2x - 4)(x + 17)
104 testQuadRootsReal(reporter, "two roots 2x^2 + 30x - 68",
105 2, 30, -68,
106 {-17, 2});
107
108 testQuadRootsReal(reporter, "two roots x^2 - 5",
109 1, 0, -5,
110 {-2.236067977499790,
111 //-2.236067977499789696409174 from Wolfram Alpha
112 2.236067977499790,
113 // 2.236067977499789696409174 from Wolfram Alpha
114 });
115
116 testQuadRootsReal(reporter, "one root x^2 - 2x + 1",
117 1, -2, 1,
118 {1});
119
120 testQuadRootsReal(reporter, "no roots 5x^2 + 6x + 7",
121 5, 6, 7,
122 {});
123
124 testQuadRootsReal(reporter, "no roots 4x^2 + 1",
125 4, 0, 1,
126 {});
127
128 testQuadRootsReal(reporter, "one root is zero, another is big",
129 14, -13, 0,
130 {0,
131 0.9285714285714286
132 //0.9285714285714285714285714 from Wolfram Alpha
133 });
134
135 // Values from a failing test case observed during testing.
136 testQuadRootsReal(reporter, "one root is zero, another is small",
137 0.2929016490705016, 0.0000030451558069, 0,
138 {-0.00001039651301576329, 0});
139
140 testQuadRootsReal(reporter, "b and c are zero, a is positive 4x^2",
141 4, 0, 0,
142 {0});
143
144 testQuadRootsReal(reporter, "b and c are zero, a is negative -4x^2",
145 -4, 0, 0,
146 {0});
147
148 testQuadRootsReal(reporter, "a and b are huge, c is zero",
149 4.3719914983870202e+291, 1.0269509510194551e+152, 0,
150 // One solution is 0, the other is so close to zero it returns
151 // true for sk_double_nearly_zero, so it is collapsed into one.
152 {0});
153
154 testQuadRootsReal(reporter, "Very small A B, very large C",
155 0x1p-1055, 0x1.3000006p-1044, -0x1.c000008p+1009,
156 // The roots are not in the range of doubles.
157 {});
158 }
159
DEF_TEST(QuadRootsReal_Linear,reporter)160 DEF_TEST(QuadRootsReal_Linear, reporter) {
161 testQuadRootsReal(reporter, "positive slope 5x + 6",
162 0, 5, 6,
163 {-1.2});
164
165 testQuadRootsReal(reporter, "negative slope -3x - 9",
166 0, -3, -9,
167 {-3.});
168 }
169
DEF_TEST(QuadRootsReal_Constant,reporter)170 DEF_TEST(QuadRootsReal_Constant, reporter) {
171 testQuadRootsReal(reporter, "No intersections y = -10",
172 0, 0, -10,
173 {});
174
175 testQuadRootsReal(reporter, "Infinite solutions y = 0",
176 0, 0, 0,
177 {0.});
178 }
179
DEF_TEST(QuadRootsReal_NonFiniteNumbers,reporter)180 DEF_TEST(QuadRootsReal_NonFiniteNumbers, reporter) {
181 // The Pathops implementation does not check for infinities nor nans in all cases.
182 double roots[2];
183 REPORTER_ASSERT(reporter,
184 SkQuads::RootsReal(DBL_MAX, 0, DBL_MAX, roots) == 0,
185 "Discriminant is negative infinity"
186 );
187 REPORTER_ASSERT(reporter,
188 SkQuads::RootsReal(DBL_MAX, DBL_MAX, DBL_MAX, roots) == 0,
189 "Double Overflow"
190 );
191
192 REPORTER_ASSERT(reporter,
193 SkQuads::RootsReal(1, NAN, -3, roots) == 0,
194 "Nan quadratic"
195 );
196 REPORTER_ASSERT(reporter,
197 SkQuads::RootsReal(0, NAN, 3, roots) == 0,
198 "Nan linear"
199 );
200 REPORTER_ASSERT(reporter,
201 SkQuads::RootsReal(0, 0, NAN, roots) == 0,
202 "Nan constant"
203 );
204 }
205
206 // Test the discriminant using
207 // Use quadratics of the form F_n * x^2 - 2 * F_(n-1) * x + F_(n-2).
208 // This has a discriminant of F_(n-1)^2 - F_n * F_(n-2) = 1 if n is even else -1.
DEF_TEST(QuadDiscriminant_Fibonacci,reporter)209 DEF_TEST(QuadDiscriminant_Fibonacci, reporter) {
210 // n, n-1, n-2
211 int64_t F[] = {1, 1, 0};
212 // F_79 just fits in the 53 significant bits of a double.
213 for (int i = 2; i < 79; ++i) {
214 F[0] = F[1] + F[2];
215
216 const int expectedDiscriminant = i % 2 == 0 ? 1 : -1;
217 REPORTER_ASSERT(reporter, SkQuads::Discriminant(F[0], F[1], F[2]) == expectedDiscriminant);
218
219 F[2] = F[1];
220 F[1] = F[0];
221 }
222 }
223
DEF_TEST(QuadRoots_Basic,reporter)224 DEF_TEST(QuadRoots_Basic, reporter) {
225 {
226 // (x - 1) (x - 1) normal quadratic form A = 1, B = 2, C =1.
227 auto [discriminant, r0, r1] = SkQuads::Roots(1, -0.5 * -2, 1);
228 REPORTER_ASSERT(reporter, discriminant == 0);
229 REPORTER_ASSERT(reporter, r0 == 1 && r1 == 1);
230 }
231
232 {
233 // (x + 2) (x + 2) normal quadratic form A = 1, B = 4, C = 4.
234 auto [discriminant, r0, r1] = SkQuads::Roots(1, -0.5 * 4, 4);
235 REPORTER_ASSERT(reporter, discriminant == 0);
236 REPORTER_ASSERT(reporter, r0 == -2 && r1 == -2);
237 }
238 }
239
240 // Test the roots using
241 // Use quadratics of the form F_n * x^2 - 2 * F_(n-1) * x + F_(n-2).
242 // The roots are (F_(n–1) ± 1)/F_n if n is even otherwise there are no roots.
DEF_TEST(QuadRoots_Fibonacci,reporter)243 DEF_TEST(QuadRoots_Fibonacci, reporter) {
244 // n, n-1, n-2
245 int64_t F[] = {1, 1, 0};
246 // F_79 just fits in the 53 significant bits of a double.
247 for (int i = 2; i < 79; ++i) {
248 F[0] = F[1] + F[2];
249
250 const int expectedDiscriminant = i % 2 == 0 ? 1 : -1;
251 auto [discriminant, r0, r1] = SkQuads::Roots(F[0], F[1], F[2]);
252 REPORTER_ASSERT(reporter, discriminant == expectedDiscriminant);
253
254 // There are only real roots when i is even.
255 if (i % 2 == 0) {
256 const double expectedLittle = ((double)F[1] - 1) / F[0];
257 const double expectedBig = ((double)F[1] + 1) / F[0];
258 if (r0 <= r1) {
259 REPORTER_ASSERT(reporter, r0 == expectedLittle);
260 REPORTER_ASSERT(reporter, r1 == expectedBig);
261 } else {
262 REPORTER_ASSERT(reporter, r1 == expectedLittle);
263 REPORTER_ASSERT(reporter, r0 == expectedBig);
264 }
265 } else {
266 REPORTER_ASSERT(reporter, std::isnan(r0));
267 REPORTER_ASSERT(reporter, std::isnan(r1));
268 }
269
270 F[2] = F[1];
271 F[1] = F[0];
272 }
273 }
274
275 // These are test cases used in the paper "The Ins and Outs of Solving Quadratic Equations with
276 // Floating-Point Arithmetic" located at
277 // https://github.com/goualard-f/QuadraticEquation.jl/blob/main/test/tests.jl
278
279 struct TestCase {
280 const double A;
281 const double B;
282 const double C;
283 const double answerLo;
284 const double answerHi;
285 };
286
DEF_TEST(QuadRoots_Hard,reporter)287 DEF_TEST(QuadRoots_Hard, reporter) {
288 const double nan = std::numeric_limits<double>::quiet_NaN();
289 const double infinity = std::numeric_limits<double>::infinity();
290
291 auto specialEqual = [] (double actual, double test) {
292 if (std::isnan(actual)) {
293 return std::isnan(test);
294 }
295
296 if (std::isinf(actual)) {
297 return std::isinf(test);
298 }
299
300 // Comparison function from the paper "The Ins and Outs ...."
301 const double errorFactor = std::sqrt(std::numeric_limits<double>::epsilon());
302 return std::abs(test - actual) <= errorFactor * std::max(std::abs(test), std::abs(actual));
303 };
304
305 auto p2 = [](double a) {
306 return std::exp2(a);
307 };
308
309 TestCase cases[] = {
310 // no real solutions
311 {2, 0, 3, nan, nan},
312 {1, 1, 1, nan, nan},
313 {2.0 * p2(600), 0, 2.0 * p2(600), nan, nan},
314 {-2.0 * p2(600), 0, -2.0 * p2(600), nan, nan},
315
316 // degenerate cases
317 {0, 0, 0, infinity, infinity},
318 {0, 1, 0, 0, 0},
319 {0, 1, 2, -2, -2},
320 {0, 3, 4, -4.0/3.0, -4.0/3.0},
321 {0, p2(600), -p2(600), 1, 1},
322 {0, p2(600), p2(600), -1, -1},
323 {0, p2(-600), p2(600), -infinity, -infinity},
324 {0, p2(600), p2(-600), 0, 0},
325 {0, 2, -1.0e-323, 5.0e-324, 5.0e-324},
326 {3, 0, 0, 0, 0},
327 {p2(600), 0, 0, 0, 0},
328 {2, 0, -3, -sqrt(3.0/2.0), sqrt(3.0/2.0)},
329 // {p2(600), 0, -p2(600), -1, 1}, determinant is infinity
330 {3, 2, 0, -2.0/3.0, 0},
331 // {p2(600), p2(700), 0, -p2(100), 0},
332 {p2(-600), p2(700), 0, -infinity, 0},
333 {p2(600), p2(-700), 0, 0, 0},
334
335 // two solutions tests
336 {1, -1, -1, -0.6180339887498948, 1.618033988749895},
337 {1, 1 + p2(-52), 0.25 + p2(-53), (-1 - p2(-51)) / 2.0, -0.5},
338 {1, p2(-511) + p2(-563), std::exp2(-1024), -7.458340888372987e-155,-7.458340574027429e-155},
339 {1, p2(27), 0.75, -134217728.0, -5.587935447692871e-09},
340 {1, -1e9, 1, 1e-09, 1000000000.0},
341 // {1.3407807929942596e154, -1.3407807929942596e154, -1.3407807929942596e154, -0.6180339887498948, 1.618033988749895},
342 {p2(600), 0.5, -p2(-600), -3.086568504549085e-181, 1.8816085719976428e-181},
343 // {p2(600), 0.5, -p2(600), -1.0, 1.0},
344 // {8.0, p2(800),-p2(500), -8.335018041099818e+239, 4.909093465297727e-91},
345 {1, p2(26), -0.125, -67108864.0, 1.862645149230957e-09},
346 // {p2(-1073), -p2(-1073), -p2(-1073), -0.6180339887498948,1.618033988749895},
347 {p2(600), -p2(-600), -p2(-600), -2.409919865102884e-181, 2.409919865102884e-181},
348
349 // Tests in Nivergelt paper
350 {-158114166017, 316227766017, -158113600000, 0.99999642020057874, 1},
351 {-312499999999.0, 707106781186.0, -400000000000.0, 1.131369396027, 1.131372303775},
352 {-67, 134, -65, 0.82722631488372798, 1.17277368511627202},
353 {0.247260273973, 0.994520547945, -0.138627953316, -4.157030027041105, 0.1348693622211607},
354 {1, -2300000, 2.0e11, 90518.994979145, 2209481.005020854},
355 {1.5*p2(-1026), 0, -p2(1022), -1.4678102981723264e308, 1.4678102981723264e308},
356
357 // one solution tests
358 {1.5*p2(-1026), 0, -p2(1022), -1.4678102981723264e308, 1.4678102981723264e308},
359 };
360
361 for (auto testCase : cases) {
362 double A = testCase.A,
363 B = testCase.B,
364 C = testCase.C,
365 answerLo = testCase.answerLo,
366 answerHi = testCase.answerHi;
367 if (SkIsFinite(answerLo, answerHi)) {
368 SkASSERT(answerLo <= answerHi);
369 }
370 auto [discriminate, r0, r1] = SkQuads::Roots(A, -0.5*B, C);
371 double rLo = std::min(r0, r1),
372 rHi = std::max(r0, r1);
373 REPORTER_ASSERT(reporter, specialEqual(rLo, answerLo));
374 REPORTER_ASSERT(reporter, specialEqual(rHi, answerHi));
375 }
376 }
377
378