Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #include "NumericalAlgorithms/RootFinding/QuadraticEquation.hpp"
5 :
6 : #include <gsl/gsl_poly.h>
7 : #include <limits>
8 :
9 : #include "DataStructures/DataVector.hpp"
10 : #include "Utilities/EqualWithinRoundoff.hpp"
11 : #include "Utilities/ErrorHandling/Assert.hpp"
12 : #include "Utilities/ErrorHandling/Error.hpp"
13 : #include "Utilities/GenerateInstantiations.hpp"
14 :
15 1 : double positive_root(const double a, const double b, const double c) noexcept {
16 : const auto roots = real_roots(a, b, c);
17 : ASSERT(roots[0] <= 0.0 and roots[1] >= 0.0,
18 : "There are two positive roots, " << roots[0] << " and " << roots[1]
19 : << ", with a=" << a << " b=" << b << " c=" << c);
20 : return roots[1];
21 : }
22 :
23 1 : std::array<double, 2> real_roots(const double a, const double b,
24 : const double c) noexcept {
25 : double x0 = std::numeric_limits<double>::signaling_NaN();
26 : double x1 = std::numeric_limits<double>::signaling_NaN();
27 : // clang-tidy: value stored ... never read (true if in Release Build)
28 : // NOLINTNEXTLINE
29 : const int num_real_roots = gsl_poly_solve_quadratic(a, b, c, &x0, &x1);
30 : ASSERT(num_real_roots == 2,
31 : "There are only " << num_real_roots << " real roots with a=" << a
32 : << " b=" << b << " c=" << c);
33 : return {{x0, x1}};
34 : }
35 :
36 : namespace detail {
37 : template <typename T>
38 : struct root_between_values_impl;
39 : enum class RootToChoose { min, max };
40 : } // namespace detail
41 :
42 : template <typename T>
43 1 : T smallest_root_greater_than_value_within_roundoff(
44 : const T& a, const T& b, const T& c, const double value) noexcept {
45 : return detail::root_between_values_impl<T>::function(
46 : a, b, c, value, std::numeric_limits<double>::max(),
47 : detail::RootToChoose::min);
48 : }
49 :
50 : template <typename T>
51 1 : T largest_root_between_values_within_roundoff(const T& a, const T& b,
52 : const T& c,
53 : const double min_value,
54 : const double max_value) noexcept {
55 : return detail::root_between_values_impl<T>::function(
56 : a, b, c, min_value, max_value, detail::RootToChoose::max);
57 : }
58 :
59 : namespace detail {
60 : template <>
61 : struct root_between_values_impl<double> {
62 : static double function(const double a, const double b, const double c,
63 : const double min_value, const double max_value,
64 : const RootToChoose min_or_max) noexcept {
65 : // Roots are returned in increasing order.
66 : const auto roots = real_roots(a, b, c);
67 :
68 : const std::array<bool, 2> roots_out_of_bounds_low{
69 : {roots[0] < min_value and
70 : not equal_within_roundoff(roots[0], min_value),
71 : roots[1] < min_value and
72 : not equal_within_roundoff(roots[1], min_value)}};
73 : const std::array<bool, 2> roots_out_of_bounds_high{
74 : {roots[0] > max_value and
75 : not equal_within_roundoff(roots[0], max_value),
76 : roots[1] > max_value and
77 : not equal_within_roundoff(roots[1], max_value)}};
78 :
79 : double return_value = std::numeric_limits<double>::signaling_NaN();
80 : const auto error_message = [&a, &b, &c,
81 : &roots](const std::string& message) noexcept {
82 : ERROR(message << " Roots are " << roots[0] << " and " << roots[1]
83 : << ", with a=" << a << " b=" << b << " c=" << c);
84 : };
85 :
86 : if (min_or_max == RootToChoose::min) {
87 : // Check roots[0] first because it is the smallest
88 : if (roots_out_of_bounds_low[0]) {
89 : if (roots_out_of_bounds_low[1]) {
90 : error_message("No root >= (within roundoff) min_value.");
91 : }
92 : if (roots_out_of_bounds_high[1]) {
93 : error_message(
94 : "No root between min_value and max_value (within roundoff).");
95 : }
96 : return_value = roots[1];
97 : } else {
98 : if (roots_out_of_bounds_high[0]) {
99 : error_message("No root <= (within roundoff) max_value.");
100 : }
101 : return_value = roots[0];
102 : }
103 : } else {
104 : // Check roots[1] first because it is the largest
105 : if (roots_out_of_bounds_high[1]) {
106 : if (roots_out_of_bounds_high[0]) {
107 : error_message("No root <= (within roundoff) max_value.");
108 : }
109 : if (roots_out_of_bounds_low[0]) {
110 : error_message(
111 : "No root between min_value and max_value (within "
112 : "roundoff).");
113 : }
114 : return_value = roots[0];
115 : } else {
116 : if (roots_out_of_bounds_low[1]) {
117 : error_message("No root >= (within roundoff) min_value.");
118 : }
119 : return_value = roots[1];
120 : }
121 : }
122 : return return_value;
123 : }
124 : };
125 :
126 : template <>
127 : struct root_between_values_impl<DataVector> {
128 : static DataVector function(const DataVector& a, const DataVector& b,
129 : const DataVector& c, const double min_value,
130 : const double max_value,
131 : const RootToChoose min_or_max) noexcept {
132 : ASSERT(a.size() == b.size(),
133 : "Size mismatch a vs b: " << a.size() << " " << b.size());
134 : ASSERT(a.size() == c.size(),
135 : "Size mismatch a vs c: " << a.size() << " " << c.size());
136 : DataVector result(a.size());
137 : for (size_t i = 0; i < a.size(); ++i) {
138 : result[i] = root_between_values_impl<double>::function(
139 : a[i], b[i], c[i], min_value, max_value, min_or_max);
140 : }
141 : return result;
142 : }
143 : };
144 : } // namespace detail
145 :
146 : // Explicit instantiations
147 : /// \cond
148 : #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
149 :
150 : #define INSTANTIATE(_, data) \
151 : template DTYPE(data) smallest_root_greater_than_value_within_roundoff( \
152 : const DTYPE(data) & a, const DTYPE(data) & b, const DTYPE(data) & c, \
153 : double value) noexcept; \
154 : template DTYPE(data) largest_root_between_values_within_roundoff( \
155 : const DTYPE(data) & a, const DTYPE(data) & b, const DTYPE(data) & c, \
156 : double min_value, double max_value) noexcept;
157 :
158 : GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector))
159 :
160 : #undef DTYPE
161 : #undef INSTANTIATE
162 : /// \endcond
|
