Flutter Engine
The Flutter Engine
SkQuads.cpp
Go to the documentation of this file.
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
12
13#include <cmath>
14#include <limits>
15
16// Solve 0 = M * x + B. If M is 0, there are no solutions, unless B is also 0,
17// in which case there are infinite solutions, so we just return 1 of them.
18static int solve_linear(const double M, const double B, double solution[2]) {
20 solution[0] = 0;
22 return 1;
23 }
24 return 0;
25 }
26 solution[0] = -B / M;
27 if (!std::isfinite(solution[0])) {
28 return 0;
29 }
30 return 1;
31}
32
33// When B >> A, then the x^2 component doesn't contribute much to the output, so the second root
34// will be very large, but have massive round off error. Because of the round off error, the
35// second root will not evaluate to zero when substituted back into the quadratic equation. In
36// the situation when B >> A, then just treat the quadratic as a linear equation.
37static bool close_to_linear(double A, double B) {
38 if (A != 0) {
39 // Return if B is much bigger than A.
40 return std::abs(B / A) >= 1.0e+16;
41 }
42
43 // Otherwise A is zero, and the quadratic is linear.
44 return true;
45}
46
47double SkQuads::Discriminant(const double a, const double b, const double c) {
48 const double b2 = b * b;
49 const double ac = a * c;
50
51 // Calculate the rough discriminate which may suffer from a loss in precision due to b2 and
52 // ac being too close.
53 const double roughDiscriminant = b2 - ac;
54
55 // We would like the calculated discriminant to have a relative error of 2-bits or less. For
56 // doubles, this means the relative error is <= E = 3*2^-53. This gives a relative error
57 // bounds of:
58 //
59 // |D - D~| / |D| <= E,
60 //
61 // where D = B*B - AC, and D~ is the floating point approximation of D.
62 // Define the following equations
63 // B2 = B*B,
64 // B2~ = B2(1 + eB2), where eB2 is the floating point round off,
65 // AC = A*C,
66 // AC~ = AC(1 + eAC), where eAC is the floating point round off, and
67 // D~ = B2~ - AC~.
68 // We can now rewrite the above bounds as
69 //
70 // |B2 - AC - (B2~ - AC~)| / |B2 - AC| = |B2 - AC - B2~ + AC~| / |B2 - AC| <= E.
71 //
72 // Substituting B2~ and AC~, and canceling terms gives
73 //
74 // |eAC * AC - eB2 * B2| / |B2 - AC| <= max(|eAC|, |eBC|) * (|AC| + |B2|) / |B2 - AC|.
75 //
76 // We know that B2 is always positive, if AC is negative, then there is no cancellation
77 // problem, and max(|eAC|, |eBC|) <= 2^-53, thus
78 //
79 // 2^-53 * (AC + B2) / |B2 - AC| <= 3 * 2^-53. Leading to
80 // AC + B2 <= 3 * |B2 - AC|.
81 //
82 // If 3 * |B2 - AC| >= AC + B2 holds, then the roughDiscriminant has 2-bits of rounding error
83 // or less and can be used.
84 if (3 * std::abs(roughDiscriminant) >= b2 + ac) {
85 return roughDiscriminant;
86 }
87
88 // Use the extra internal precision afforded by fma to calculate the rounding error for
89 // b^2 and ac.
90 const double b2RoundingError = std::fma(b, b, -b2);
91 const double acRoundingError = std::fma(a, c, -ac);
92
93 // Add the total rounding error back into the discriminant guess.
94 const double discriminant = (b2 - ac) + (b2RoundingError - acRoundingError);
95 return discriminant;
96}
97
98SkQuads::RootResult SkQuads::Roots(double A, double B, double C) {
99 const double discriminant = Discriminant(A, B, C);
100
101 if (A == 0) {
102 double root;
103 if (B == 0) {
104 if (C == 0) {
105 root = std::numeric_limits<double>::infinity();
106 } else {
107 root = std::numeric_limits<double>::quiet_NaN();
108 }
109 } else {
110 // Solve -2*B*x + C == 0; x = c/(2*b).
111 root = C / (2 * B);
112 }
113 return {discriminant, root, root};
114 }
115
116 SkASSERT(A != 0);
117 if (discriminant == 0) {
118 return {discriminant, B / A, B / A};
119 }
120
121 if (discriminant > 0) {
122 const double D = sqrt(discriminant);
123 const double R = B > 0 ? B + D : B - D;
124 return {discriminant, R / A, C / R};
125 }
126
127 // The discriminant is negative or is not finite.
128 return {discriminant, NAN, NAN};
129}
130
131static double zero_if_tiny(double x) {
132 return sk_double_nearly_zero(x) ? 0 : x;
133}
134
135int SkQuads::RootsReal(const double A, const double B, const double C, double solution[2]) {
136
137 if (close_to_linear(A, B)) {
138 return solve_linear(B, C, solution);
139 }
140
141 SkASSERT(A != 0);
142 auto [discriminant, root0, root1] = Roots(A, -0.5 * B, C);
143
144 // Handle invariants to mesh with existing code from here on.
145 if (!std::isfinite(discriminant) || discriminant < 0) {
146 return 0;
147 }
148
149 int roots = 0;
150 if (const double r0 = zero_if_tiny(root0); std::isfinite(r0)) {
151 solution[roots++] = r0;
152 }
153 if (const double r1 = zero_if_tiny(root1); std::isfinite(r1)) {
154 solution[roots++] = r1;
155 }
156
157 if (roots == 2 && sk_doubles_nearly_equal_ulps(solution[0], solution[1])) {
158 roots = 1;
159 }
160
161 return roots;
162}
163
164double SkQuads::EvalAt(double A, double B, double C, double t) {
165 // Use fused-multiply-add to reduce the amount of round-off error between terms.
166 return std::fma(std::fma(A, t, B), t, C);
167}
#define SkASSERT(cond)
Definition: SkAssert.h:116
bool sk_double_nearly_zero(double a)
bool sk_doubles_nearly_equal_ulps(double a, double b, uint8_t maxUlpsDiff=16)
static skvx::float4 fma(const skvx::float4 &f, float m, const skvx::float4 &a)
Definition: SkGeometry.cpp:597
static bool close_to_linear(double A, double B)
Definition: SkQuads.cpp:37
static int solve_linear(const double M, const double B, double solution[2])
Definition: SkQuads.cpp:18
static double zero_if_tiny(double x)
Definition: SkQuads.cpp:131
static double Discriminant(double A, double B, double C)
Definition: SkQuads.cpp:47
static RootResult Roots(double A, double B, double C)
Definition: SkQuads.cpp:98
static int RootsReal(double A, double B, double C, double solution[2])
Definition: SkQuads.cpp:135
static double EvalAt(double A, double B, double C, double t)
Definition: SkQuads.cpp:164
static bool b
struct MyStruct a[10]
#define R(r)
#define B
double x
string root
Definition: scale_cpu.py:20
SINT bool isfinite(const Vec< N, T > &v)
Definition: SkVx.h:1003
SIN Vec< N, float > abs(const Vec< N, float > &x)
Definition: SkVx.h:707
SIN Vec< N, float > sqrt(const Vec< N, float > &x)
Definition: SkVx.h:706
#define M(PROC, DITHER)