Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <limits>
8 : #include <type_traits>
9 :
10 : #include "Utilities/Algorithm.hpp"
11 : #include "Utilities/ConstantExpressions.hpp"
12 : #include "Utilities/ErrorHandling/Assert.hpp"
13 : #include "Utilities/ForceInline.hpp"
14 : #include "Utilities/TypeTraits/IsIterable.hpp"
15 : #include "Utilities/TypeTraits/IsMaplike.hpp"
16 :
17 : namespace EqualWithinRoundoffImpls {
18 : /*!
19 : * \brief Specialize this class to add support for the `equal_within_roundoff`
20 : * function.
21 : *
22 : * Ensure the `Lhs` and `Rhs` are symmetric. A specialization must implement a
23 : * static `apply` function with this signature:
24 : *
25 : * ```cpp
26 : * static bool apply(const Lhs& lhs, const Rhs& rhs, const double eps,
27 : * const double scale);
28 : * ```
29 : *
30 : * It can be helpful to invoke the `equal_within_roundoff` function for floating
31 : * points from within your specialization.
32 : */
33 : template <typename Lhs, typename Rhs, typename = std::nullptr_t>
34 1 : struct EqualWithinRoundoffImpl;
35 : } // namespace EqualWithinRoundoffImpls
36 :
37 : /*!
38 : * \ingroup UtilitiesGroup
39 : * \brief Checks if two values `lhs` and `rhs` are equal within roundoff, by
40 : * comparing `abs(lhs - rhs) < (max(abs(lhs), abs(rhs)) + scale) * eps`.
41 : *
42 : * The two values can be floating-point numbers, or any types for which
43 : * `EqualWithinRoundoffImpls::EqualWithinRoundoffImpl` has been specialized. For
44 : * example, a default implementation exists for the case where `lhs`, `rhs`, or
45 : * both, are iterable, and compares the values point-wise.
46 : */
47 : template <typename Lhs, typename Rhs>
48 1 : constexpr SPECTRE_ALWAYS_INLINE bool equal_within_roundoff(
49 : const Lhs& lhs, const Rhs& rhs,
50 : const double eps = std::numeric_limits<double>::epsilon() * 100.0,
51 : const double scale = 1.0) {
52 : return EqualWithinRoundoffImpls::EqualWithinRoundoffImpl<Lhs, Rhs>::apply(
53 : lhs, rhs, eps, scale);
54 : }
55 :
56 : /// Specializations of `EqualWithinRoundoffImpl` for custom types, to add
57 : /// support for the `equal_within_roundoff` function.
58 : namespace EqualWithinRoundoffImpls {
59 :
60 : // Compare two floating points
61 : template <typename Floating>
62 : struct EqualWithinRoundoffImpl<Floating, Floating,
63 : Requires<std::is_floating_point_v<Floating>>> {
64 : static constexpr SPECTRE_ALWAYS_INLINE bool apply(const Floating& lhs,
65 : const Floating& rhs,
66 : const double eps,
67 : const double scale) {
68 : return ce_fabs(lhs - rhs) <=
69 : (std::max(ce_fabs(lhs), ce_fabs(rhs)) + scale) * eps;
70 : }
71 : };
72 :
73 : // Compare a complex number to a floating point, interpreting the latter as a
74 : // real number
75 : template <typename Floating>
76 : struct EqualWithinRoundoffImpl<std::complex<Floating>, Floating,
77 : Requires<std::is_floating_point_v<Floating>>> {
78 : static SPECTRE_ALWAYS_INLINE bool apply(const std::complex<Floating>& lhs,
79 : const Floating& rhs, const double eps,
80 : const double scale) {
81 : return equal_within_roundoff(lhs.real(), rhs, eps, scale) and
82 : equal_within_roundoff(lhs.imag(), 0., eps, scale);
83 : }
84 : };
85 :
86 : // Compare a floating point to a complex number, interpreting the former as a
87 : // real number
88 : template <typename Floating>
89 : struct EqualWithinRoundoffImpl<Floating, std::complex<Floating>,
90 : Requires<std::is_floating_point_v<Floating>>> {
91 : static SPECTRE_ALWAYS_INLINE bool apply(const std::complex<Floating>& lhs,
92 : const Floating& rhs, const double eps,
93 : const double scale) {
94 : return equal_within_roundoff(rhs, lhs, eps, scale);
95 : }
96 : };
97 :
98 : // Compare two complex numbers
99 : template <typename Floating>
100 : struct EqualWithinRoundoffImpl<std::complex<Floating>, std::complex<Floating>,
101 : Requires<std::is_floating_point_v<Floating>>> {
102 : static SPECTRE_ALWAYS_INLINE bool apply(const std::complex<Floating>& lhs,
103 : const std::complex<Floating>& rhs,
104 : const double eps,
105 : const double scale) {
106 : return equal_within_roundoff(lhs.real(), rhs.real(), eps, scale) and
107 : equal_within_roundoff(lhs.imag(), rhs.imag(), eps, scale);
108 : }
109 : };
110 :
111 : // Compare an iterable to a floating point
112 : template <typename Lhs, typename Rhs>
113 : struct EqualWithinRoundoffImpl<
114 : Lhs, Rhs,
115 : Requires<tt::is_iterable_v<Lhs> and not tt::is_maplike_v<Lhs> and
116 : std::is_floating_point_v<Rhs>>> {
117 : static SPECTRE_ALWAYS_INLINE bool apply(const Lhs& lhs, const Rhs& rhs,
118 : const double eps,
119 : const double scale) {
120 : return alg::all_of(lhs, [&rhs, &eps, &scale](const auto& lhs_element) {
121 : return equal_within_roundoff(lhs_element, rhs, eps, scale);
122 : });
123 : }
124 : };
125 :
126 : // Compare a floating point to an iterable
127 : template <typename Lhs, typename Rhs>
128 : struct EqualWithinRoundoffImpl<
129 : Lhs, Rhs,
130 : Requires<tt::is_iterable_v<Rhs> and not tt::is_maplike_v<Rhs> and
131 : std::is_floating_point_v<Lhs>>> {
132 : static SPECTRE_ALWAYS_INLINE bool apply(const Lhs& lhs, const Rhs& rhs,
133 : const double eps,
134 : const double scale) {
135 : return equal_within_roundoff(rhs, lhs, eps, scale);
136 : }
137 : };
138 :
139 : // Compare two iterables
140 : template <typename Lhs, typename Rhs>
141 0 : struct EqualWithinRoundoffImpl<
142 : Lhs, Rhs,
143 : Requires<tt::is_iterable_v<Lhs> and not tt::is_maplike_v<Lhs> and
144 : tt::is_iterable_v<Rhs> and not tt::is_maplike_v<Rhs>>> {
145 0 : static bool apply(const Lhs& lhs, const Rhs& rhs, const double eps,
146 : const double scale) {
147 : auto lhs_it = lhs.begin();
148 : auto rhs_it = rhs.begin();
149 : while (lhs_it != lhs.end() and rhs_it != rhs.end()) {
150 : if (not equal_within_roundoff(*lhs_it, *rhs_it, eps, scale)) {
151 : return false;
152 : }
153 : ++lhs_it;
154 : ++rhs_it;
155 : }
156 : ASSERT(lhs_it == lhs.end() and rhs_it == rhs.end(),
157 : "Can't compare lhs and rhs because they have different lengths.");
158 : return true;
159 : }
160 : };
161 :
162 : } // namespace EqualWithinRoundoffImpls
|
