27 if (!std::isfinite(solution[0])) {
40 return std::abs(
B /
A) >= 1.0e+16;
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();
113 return {discriminant, root, root};
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);
145 if (!std::isfinite(discriminant) || discriminant < 0) {
150 if (
const double r0 =
zero_if_tiny(root0); std::isfinite(r0)) {
151 solution[roots++] = r0;
153 if (
const double r1 =
zero_if_tiny(root1); std::isfinite(r1)) {
154 solution[roots++] = r1;
166 return std::fma(std::fma(
A, t,
B), t,
C);
bool sk_double_nearly_zero(double a)
bool sk_doubles_nearly_equal_ulps(double a, double b, uint8_t maxUlpsDiff=16)
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)