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 #include "src/base/SkQuads.h"
8
9 #include "include/private/base/SkFloatingPoint.h"
10
11 #include <cmath>
12
13 // Solve 0 = M * x + B. If M is 0, there are no solutions, unless B is also 0,
14 // in which case there are infinite solutions, so we just return 1 of them.
solve_linear(const double M,const double B,double solution[2])15 static int solve_linear(const double M, const double B, double solution[2]) {
16 if (sk_double_nearly_zero(M)) {
17 solution[0] = 0;
18 if (sk_double_nearly_zero(B)) {
19 return 1;
20 }
21 return 0;
22 }
23 solution[0] = -B / M;
24 if (!std::isfinite(solution[0])) {
25 return 0;
26 }
27 return 1;
28 }
29
30 // When the A coefficient of a quadratic is close to 0, there can be floating point error
31 // that arises from computing a very large root. In those cases, we would rather be
32 // precise about the one smaller root, so we have this arbitrary cutoff for when A is
33 // really small or small compared to B.
close_to_linear(double A,double B)34 static bool close_to_linear(double A, double B) {
35 if (sk_double_nearly_zero(B)) {
36 return sk_double_nearly_zero(A);
37 }
38 // This is a different threshold (tighter) than the close_to_a_quadratic in SkCubics.cpp
39 // because the SkQuads::RootsReal gives better answers for longer as A/B -> 0.
40 return std::abs(A / B) < 1.0e-16;
41 }
42
RootsReal(const double A,const double B,const double C,double solution[2])43 int SkQuads::RootsReal(const double A, const double B, const double C, double solution[2]) {
44 if (close_to_linear(A, B)) {
45 return solve_linear(B, C, solution);
46 }
47 // If A is zero (e.g. B was nan and thus close_to_linear was false), we will
48 // temporarily have infinities rolling about, but will catch that when checking
49 // p2 - q.
50 const double p = sk_ieee_double_divide(B, 2 * A);
51 const double q = sk_ieee_double_divide(C, A);
52 /* normal form: x^2 + px + q = 0 */
53 const double p2 = p * p;
54 if (!std::isfinite(p2 - q) ||
55 (!sk_double_nearly_zero(p2 - q) && p2 < q)) {
56 return 0;
57 }
58 double sqrt_D = 0;
59 if (p2 > q) {
60 sqrt_D = sqrt(p2 - q);
61 }
62 solution[0] = sqrt_D - p;
63 solution[1] = -sqrt_D - p;
64 if (sk_double_nearly_zero(sqrt_D) ||
65 sk_doubles_nearly_equal_ulps(solution[0], solution[1])) {
66 return 1;
67 }
68 return 2;
69 }
70