4  * Use of this source code is governed by a BSD-style license that can be
26     solution[0] = -B / M;  in solve_linear()
49     const double ac = a * c;  in Discriminant()  local
52     // ac being too close.  in Discriminant()
53     const double roughDiscriminant = b2 - ac;  in Discriminant()
55     // We would like the calculated discriminant to have a relative error of 2-bits or less. For  in Discriminant()
56     // doubles, this means the relative error is <= E = 3*2^-53. This gives a relative error  in Discriminant()
59     //     |D - D~| / |D| <= E,  in Discriminant()
61     // where D = B*B - AC, and D~ is the floating point approximation of D.  in Discriminant()
65     //     AC = A*C,  in Discriminant()
66     //     AC~ = AC(1 + eAC), where eAC is the floating point round off, and  in Discriminant()
67     //     D~ = B2~ - AC~.  in Discriminant()
70     //     |B2 - AC - (B2~ - AC~)| / |B2 - AC| = |B2 - AC - B2~ + AC~| / |B2 - AC| <= E.  in Discriminant()
72     //  Substituting B2~ and AC~, and canceling terms gives  in Discriminant()
74     //     |eAC * AC - eB2 * B2| / |B2 - AC| <= max(|eAC|, |eBC|) * (|AC| + |B2|) / |B2 - AC|.  in Discriminant()
76     //  We know that B2 is always positive, if AC is negative, then there is no cancellation  in Discriminant()
77     //  problem, and max(|eAC|, |eBC|) <= 2^-53, thus  in Discriminant()
79     //     2^-53 * (AC + B2) / |B2 - AC| <= 3 * 2^-53. Leading to  in Discriminant()
80     //     AC + B2 <= 3 * |B2 - AC|.  in Discriminant()
82     // If 3 * |B2 - AC| >= AC + B2 holds, then the roughDiscriminant has 2-bits of rounding error  in Discriminant()
84     if (3 * std::abs(roughDiscriminant) >= b2 + ac) {  in Discriminant()
89     // b^2 and ac.  in Discriminant()
90     const double b2RoundingError = std::fma(b, b, -b2);  in Discriminant()
91     const double acRoundingError = std::fma(a, c, -ac);  in Discriminant()
94     const double discriminant = (b2 - ac) + (b2RoundingError - acRoundingError);  in Discriminant()
110             // Solve -2*B*x + C == 0; x = c/(2*b).  in Roots()
123         const double R = B > 0 ? B + D : B - D;  in Roots()
142     auto [discriminant, root0, root1] = Roots(A, -0.5 * B, C);  in RootsReal()
165     // Use fused-multiply-add to reduce the amount of round-off error between terms.  in EvalAt()