diff --git a/base/cpuid.cpp b/base/cpuid.cpp index 04a564bca2..158a0c376f 100644 --- a/base/cpuid.cpp +++ b/base/cpuid.cpp @@ -28,6 +28,7 @@ namespace { // This vector is not a static member variable of CPUIDFeatureFlag because we do // not want to include in the header. std::vector>& CPUIDFlags() { + std::basic_string_view s; static std::vector> result; return result; } diff --git a/numerics/root_finders.hpp b/numerics/root_finders.hpp index 92a3fba0dc..df23cdf136 100644 --- a/numerics/root_finders.hpp +++ b/numerics/root_finders.hpp @@ -65,14 +65,14 @@ Argument Brent( double eps = Sqrt(ScaleB(0.5, 1 - std::numeric_limits::digits))); // Returns the solutions of the quadratic equation: -// a2 * (x - origin)^2 + a1 * (x - origin) + a0 == 0 +// a₂ * (x - origin)² + a₁ * (x - origin) + a₀ == 0 // The result may have 0, 1 or 2 values and is sorted. template BoundedArray SolveQuadraticEquation( Argument const& origin, - Value const& a0, - Derivative const& a1, - Derivative, Argument> const& a2); + Value const& a₀, + Derivative const& a₁, + Derivative, Argument> const& a₂); } // namespace internal diff --git a/numerics/root_finders_body.hpp b/numerics/root_finders_body.hpp index c85fd2ef8e..23aa0a6797 100644 --- a/numerics/root_finders_body.hpp +++ b/numerics/root_finders_body.hpp @@ -343,46 +343,44 @@ Argument Brent(Function f, template BoundedArray SolveQuadraticEquation( Argument const& origin, - Value const& a0, - Derivative const& a1, - Derivative, Argument> const& a2) { + Value const& a₀, + Derivative const& a₁, + Derivative, Argument> const& a₂) { using Derivative1 = Derivative; using Discriminant = Square; - // This algorithm is after section 1.8 of Accuracy and Stability of Numerical - // Algorithms, Second Edition, Higham, ISBN 0-89871-521-0. + // This algorithm is after section 1.8 of [Hig02]. - static Discriminant const discriminant_zero{}; + constexpr Discriminant zero{}; // Use double precision for the discriminant because there can be - // cancellations. Higham mentions that it is necessary "to use extended - // precision (or some trick tantamount to the use of extended precision)." - DoublePrecision discriminant = TwoSum(a1 * a1, -4.0 * a0 * a2); + // cancellations. Higham mentions that it is necessary “to use extended + // precision (or some trick tantamount to the use of extended precision).” + DoublePrecision discriminant = TwoSum(a₁ * a₁, -4.0 * a₀ * a₂); - if (discriminant.value == discriminant_zero && - discriminant.error == discriminant_zero) { + if (discriminant.value == zero && discriminant.error == zero) { // One solution. - return {origin - 0.5 * a1 / a2}; - } else if (discriminant.value < discriminant_zero || - (discriminant.value == discriminant_zero && - discriminant.error < discriminant_zero)) { + return {origin - 0.5 * a₁ / a₂}; + } else if (discriminant.value < zero || + (discriminant.value == zero && discriminant.error < zero)) { // No solution. return {}; } else { - // Two solutions. Compute the numerator of the larger one. + // Two solutions. Compute the numerator of the larger one (in absolute + // value). Derivative1 numerator; - static Derivative1 derivative_zero{}; - if (a1 > derivative_zero) { - numerator = -a1 - Sqrt(discriminant.value + discriminant.error); + constexpr Derivative1 zero{}; + if (a₁ > zero) { + numerator = -a₁ - Sqrt(discriminant.value + discriminant.error); } else { - numerator = -a1 + Sqrt(discriminant.value + discriminant.error); + numerator = -a₁ + Sqrt(discriminant.value + discriminant.error); } - auto const solution1 = origin + numerator / (2.0 * a2); - auto const solution2 = origin + (2.0 * a0) / numerator; - if (solution1 < solution2) { - return {solution1, solution2}; + auto const x₁ = origin + numerator / (2.0 * a₂); + auto const x₂ = origin + (2.0 * a₀) / numerator; + if (x₁ < x₂) { + return {x₁, x₂}; } else { - return {solution2, solution1}; + return {x₂, x₁}; } } }