48 const double b2 =
b *
b;
49 const double ac =
a * c;
53 const double roughDiscriminant = b2 - ac;
84 if (3 *
std::abs(roughDiscriminant) >= b2 + ac) {
85 return roughDiscriminant;
90 const double b2RoundingError =
std::fma(
b,
b, -b2);
91 const double acRoundingError =
std::fma(
a, c, -ac);
94 const double discriminant = (b2 - ac) + (b2RoundingError - acRoundingError);
105 root = std::numeric_limits<double>::infinity();
107 root = std::numeric_limits<double>::quiet_NaN();
117 if (discriminant == 0) {
118 return {discriminant,
B /
A,
B /
A};
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};
128 return {discriminant, NAN, NAN};
142 auto [discriminant, root0, root1] =
Roots(
A, -0.5 *
B,
C);
151 solution[
roots++] = r0;
154 solution[
roots++] = r1;
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)
static bool close_to_linear(double A, double B)
static int solve_linear(const double M, const double B, double solution[2])
static double zero_if_tiny(double x)
static double Discriminant(double A, double B, double C)
static RootResult Roots(double A, double B, double C)
static int RootsReal(double A, double B, double C, double solution[2])
static double EvalAt(double A, double B, double C, double t)
SINT bool isfinite(const Vec< N, T > &v)
SIN Vec< N, float > abs(const Vec< N, float > &x)
SIN Vec< N, float > sqrt(const Vec< N, float > &x)