Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Define simple functions for constant expressions.
6 :
7 : #pragma once
8 :
9 : #include <algorithm>
10 : #include <array>
11 : #include <cassert>
12 : #include <cstddef>
13 : #include <cstdint>
14 : #include <functional>
15 : #include <initializer_list>
16 : #include <utility>
17 :
18 : #include "Utilities/ForceInline.hpp"
19 : #include "Utilities/MakeWithValue.hpp"
20 : #include "Utilities/Requires.hpp"
21 : #include "Utilities/TMPL.hpp"
22 : #include "Utilities/TypeTraits.hpp"
23 : #include "Utilities/TypeTraits/GetFundamentalType.hpp"
24 : #include "Utilities/TypeTraits/IsA.hpp"
25 : #include "Utilities/TypeTraits/IsInteger.hpp"
26 :
27 : // IWYU pragma: no_include "DataStructures/DataVector.hpp"
28 :
29 : /// \ingroup ConstantExpressionsGroup
30 : /// Compute 2 to the n for integral types.
31 : ///
32 : /// \param n the power of two to compute.
33 : /// \return 2^n
34 : template <typename T,
35 : Requires<tt::is_integer_v<T> and std::is_unsigned_v<T>> = nullptr>
36 1 : SPECTRE_ALWAYS_INLINE constexpr T two_to_the(T n) {
37 : return T(1) << n;
38 : }
39 :
40 : /// \ingroup ConstantExpressionsGroup
41 : /// Get the nth bit from the right, counting from zero.
42 : ///
43 : /// \param i the value to be acted on.
44 : /// \param N which place to extract the bit
45 : /// \return the value of the bit at that place.
46 1 : SPECTRE_ALWAYS_INLINE constexpr size_t get_nth_bit(const size_t i,
47 : const size_t N) {
48 : return (i >> N) % 2;
49 : }
50 :
51 : /// \ingroup ConstantExpressionsGroup
52 : /// \brief Compute the square of `x`
53 : template <typename T>
54 : // NOLINTNEXTLINE(readability-const-return-type)
55 1 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) square(const T& x) {
56 : return x * x;
57 : }
58 :
59 : /// \ingroup ConstantExpressionsGroup
60 : /// \brief Compute the cube of `x`
61 : template <typename T>
62 1 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) cube(const T& x) {
63 : return x * x * x;
64 : }
65 :
66 : /*!
67 : * \ingroup ConstantExpressionsGroup
68 : * \brief Compute the falling factorial of \f$(x)_{n}\f$
69 : *
70 : * See http://mathworld.wolfram.com/FallingFactorial.html
71 : * \note The largest representable factorial is 20!. It is up to the user to
72 : * ensure this is satisfied
73 : */
74 1 : constexpr uint64_t falling_factorial(const uint64_t x, const uint64_t n) {
75 : // clang-tidy: don't warn about STL internals, I can't fix them
76 : assert(n <= x); // NOLINT
77 : uint64_t r = 1;
78 : for (uint64_t k = 0; k < n; ++k) {
79 : r *= (x - k);
80 : }
81 : return r;
82 : }
83 :
84 : /*!
85 : * \ingroup ConstantExpressionsGroup
86 : * \brief Compute the factorial of \f$n!\f$
87 : */
88 1 : constexpr uint64_t factorial(const uint64_t n) {
89 : assert(n <= 20); // NOLINT
90 : return falling_factorial(n, n);
91 : }
92 :
93 : /// \ingroup ConstantExpressionsGroup
94 : namespace ConstantExpressions_detail {
95 :
96 : // Implementation functions for the pow template function below which computes
97 : // optimized powers where the exponent is a compile-time integer
98 :
99 : // base case power 0: Returns 1.0. This will need to be overloaded for types for
100 : // which the simple return is inappropriate.
101 : template <typename T>
102 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) pow_impl(
103 : const T& /*t*/, std::integral_constant<int, 0> /*meta*/,
104 : std::bool_constant<true> /*exponent_was_positive*/) {
105 : return static_cast<tt::get_fundamental_type_t<T>>(1.0);
106 : }
107 :
108 : // special case power 1: acts as a direct identity function for efficiency
109 : template <typename T>
110 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) pow_impl(
111 : const T& t, std::integral_constant<int, 1> /*meta*/,
112 : std::bool_constant<true> /*exponent_was_positive*/) {
113 : return t;
114 : }
115 :
116 : // general case for positive powers: return the power via recursive inline call,
117 : // which expands to a series of multiplication operations.
118 : template <int N, typename T>
119 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) pow_impl(
120 : const T& t, std::integral_constant<int, N> /*meta*/,
121 : std::bool_constant<true> /*exponent_was_positive*/) {
122 : return t * pow_impl(t, std::integral_constant<int, N - 1>{},
123 : std::bool_constant<true>{});
124 : }
125 :
126 : // general case for negative powers: return the multiplicative inverse of the
127 : // result from the recursive inline call, which expands to a series of
128 : // multiplication operations. This assumes that division is supported with
129 : // tt::get_fundamental_type_t<T> and T. If not, this utility will need further
130 : // specialization.
131 : template <int N, typename T>
132 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) pow_impl(
133 : const T& t, std::integral_constant<int, N> /*meta*/,
134 : std::bool_constant<false> /*exponent_was_positive*/) {
135 : return static_cast<tt::get_fundamental_type_t<T>>(1) /
136 : (pow_impl(t, std::integral_constant<int, -N>{},
137 : std::bool_constant<true>{}));
138 : }
139 : } // namespace ConstantExpressions_detail
140 :
141 : /// \ingroup ConstantExpressionsGroup
142 : /// \brief Compute t^N where N is an integer (positive or negative)
143 : ///
144 : /// \warning If passing an integer this will do integer arithmetic, e.g.
145 : /// `pow<-10>(2) == 0` evaluates to `true`
146 : ///
147 : /// \warning For optimization, it is assumed that the `pow<0>` of a vector type
148 : /// (e.g. `DataVector`) will not be used for initialization or directly indexed,
149 : /// so `pow<0>` returns simply `1.0`. In the case of use for initialization, a
150 : /// constructor should be used instead, and in the case of a direct index, the
151 : /// expression may be simplifyable to avoid the use of `pow<0>` altogether. If a
152 : /// more complete treatment of `pow<0>` is required, further overloads may be
153 : /// added to the `ConstantExpressions_detail` namespace.
154 : ///
155 : /// \tparam N the integer power being raised to in \f$t^N\f$
156 : /// \param t the value being exponentiated
157 : /// \return value \f$t^N\f$ determined via repeated multiplication
158 : template <int N, typename T>
159 1 : SPECTRE_ALWAYS_INLINE constexpr decltype(auto) pow(const T& t) {
160 : return ConstantExpressions_detail::pow_impl(
161 : t, std::integral_constant<int, N>{}, std::bool_constant<(N >= 0)>{});
162 : }
163 :
164 : /// \ingroup ConstantExpressionsGroup
165 : /// \brief Compute the absolute value of of its argument
166 : ///
167 : /// The argument must be comparable to an int and must be negatable.
168 : template <typename T, Requires<tt::is_integer_v<T> or
169 : std::is_floating_point_v<T>> = nullptr>
170 1 : SPECTRE_ALWAYS_INLINE constexpr T ce_abs(const T& x) {
171 : return x < 0 ? -x : x;
172 : }
173 :
174 : /// \cond
175 : template <>
176 : SPECTRE_ALWAYS_INLINE constexpr double ce_abs(const double& x) {
177 : return __builtin_fabs(x);
178 : }
179 :
180 : template <>
181 : SPECTRE_ALWAYS_INLINE constexpr float ce_abs(const float& x) {
182 : return __builtin_fabsf(x);
183 : }
184 : /// \endcond
185 :
186 : /// \ingroup ConstantExpressionsGroup
187 : /// \brief Compute the absolute value of its argument
188 1 : constexpr SPECTRE_ALWAYS_INLINE double ce_fabs(const double x) {
189 : return __builtin_fabs(x);
190 : }
191 :
192 0 : constexpr SPECTRE_ALWAYS_INLINE float ce_fabs(const float x) {
193 : return __builtin_fabsf(x);
194 : }
195 :
196 : namespace ConstantExpressions_detail {
197 : struct CompareByMagnitude {
198 : template <typename T>
199 : constexpr bool operator()(const T& a, const T& b) {
200 : return ce_abs(a) < ce_abs(b);
201 : }
202 : };
203 : } // namespace ConstantExpressions_detail
204 :
205 : /// \ingroup ConstantExpressionsGroup
206 : /// \brief Return the argument with the largest magnitude
207 : ///
208 : /// Magnitude is determined by constexpr_abs. In case of a tie,
209 : /// returns the leftmost of the tied values.
210 : /// @{
211 : template <typename T>
212 1 : constexpr const T& max_by_magnitude(const T& a, const T& b) {
213 : return std::max(a, b, ConstantExpressions_detail::CompareByMagnitude{});
214 : }
215 :
216 : template <typename T>
217 1 : constexpr T max_by_magnitude(std::initializer_list<T> ilist) {
218 : return std::max(std::move(ilist),
219 : ConstantExpressions_detail::CompareByMagnitude{});
220 : }
221 : /// @}
222 :
223 : /// \ingroup ConstantExpressionsGroup
224 : /// \brief Return the argument with the smallest magnitude
225 : ///
226 : /// Magnitude is determined by constexpr_abs. In case of a tie,
227 : /// returns the leftmost of the tied values.
228 : /// @{
229 : template <typename T>
230 1 : constexpr const T& min_by_magnitude(const T& a, const T& b) {
231 : return std::min(a, b, ConstantExpressions_detail::CompareByMagnitude{});
232 : }
233 :
234 : template <typename T>
235 1 : constexpr T min_by_magnitude(std::initializer_list<T> ilist) {
236 : return std::min(std::move(ilist),
237 : ConstantExpressions_detail::CompareByMagnitude{});
238 : }
239 : /// @}
240 :
241 : /// \ingroup ConstantExpressionsGroup
242 : /// \brief Returns `f(ic<0>{}) + f(ic<1>{}) + ... + f(ic<NumTerms-1>{})`
243 : /// where `ic<N>` stands for `std::integral_constant<size_t, N>`.
244 : /// This function allows the result types of each operation to be
245 : /// different, and so works efficiently with expression templates.
246 : /// \note When summing expression templates one must be careful of
247 : /// referring to temporaries in `f`.
248 : template <size_t NumTerms, typename Function, Requires<NumTerms == 1> = nullptr>
249 1 : constexpr decltype(auto) constexpr_sum(Function&& f) {
250 : return f(std::integral_constant<size_t, 0>{});
251 : }
252 :
253 : /// \cond HIDDEN_SYMBOLS
254 : template <size_t NumTerms, typename Function,
255 : Requires<(NumTerms > 1)> = nullptr>
256 : constexpr decltype(auto) constexpr_sum(Function&& f) {
257 : return constexpr_sum<NumTerms - 1>(f) +
258 : f(std::integral_constant<size_t, NumTerms - 1>{});
259 : }
260 : /// \endcond
261 :
262 : namespace detail {
263 : template <typename List, size_t... indices,
264 : Requires<not tt::is_a_v<tmpl::list, tmpl::front<List>>> = nullptr>
265 : inline constexpr std::array<std::decay_t<decltype(tmpl::front<List>::value)>,
266 : tmpl::size<List>::value>
267 : make_array_from_list_helper(
268 : std::integer_sequence<size_t, indices...> /*meta*/) {
269 : return std::array<std::decay_t<decltype(tmpl::front<List>::value)>,
270 : tmpl::size<List>::value>{
271 : {tmpl::at<List, tmpl::size_t<indices>>::value...}};
272 : }
273 : } // namespace detail
274 :
275 : /// \ingroup ConstantExpressionsGroup
276 : /// Make an array from a typelist that holds std::integral_constant's all of
277 : /// which have the same `value_type`
278 : ///
279 : /// \tparam List the typelist of std::integral_constant's
280 : /// \return array of integral values from the typelist
281 : template <typename List,
282 : Requires<not tt::is_a_v<tmpl::list, tmpl::front<List>>> = nullptr>
283 1 : inline constexpr auto make_array_from_list()
284 : -> std::array<std::decay_t<decltype(tmpl::front<List>::value)>,
285 : tmpl::size<List>::value> {
286 : return detail::make_array_from_list_helper<List>(
287 : std::make_integer_sequence<size_t, tmpl::size<List>::value>{});
288 : }
289 :
290 : template <typename TypeForZero,
291 : Requires<not tt::is_a_v<tmpl::list, TypeForZero>> = nullptr>
292 : inline constexpr std::array<std::decay_t<TypeForZero>, 0>
293 0 : make_array_from_list() {
294 : return std::array<std::decay_t<TypeForZero>, 0>{{}};
295 : }
296 :
297 : namespace detail {
298 : /// \cond
299 : template <typename List, size_t... indices,
300 : Requires<tt::is_a<tmpl::list, tmpl::front<List>>::value> = nullptr>
301 : inline constexpr std::array<
302 : std::decay_t<
303 : decltype(make_array_from_list<tmpl::at<List, tmpl::size_t<0>>>())>,
304 : tmpl::size<List>::value>
305 : make_array_from_list_helper(
306 : std::integer_sequence<size_t, indices...> /*unused*/) {
307 : return std::array<std::decay_t<decltype(make_array_from_list<
308 : tmpl::at<List, tmpl::size_t<0>>>())>,
309 : tmpl::size<List>::value>{
310 : {make_array_from_list_helper<tmpl::at<List, tmpl::size_t<indices>>>(
311 : std::make_integer_sequence<
312 : size_t,
313 : tmpl::size<tmpl::at<List, tmpl::size_t<indices>>>::value>{})...}};
314 : }
315 : /// \endcond
316 : } // namespace detail
317 :
318 : /// \ingroup ConstantExpressionsGroup
319 : ///
320 : /// Make an array of arrays from a typelist that holds typelists of
321 : /// std::integral_constant's all of which have the same `value_type`
322 : ///
323 : /// \tparam List the typelist of typelists of std::integral_constant's
324 : /// \return array of arrays of integral values from the typelists
325 : template <typename List,
326 : Requires<tt::is_a_v<tmpl::list, tmpl::front<List>>> = nullptr>
327 : inline constexpr auto make_array_from_list() {
328 : return detail::make_array_from_list_helper<List>(
329 : std::make_integer_sequence<size_t, tmpl::size<List>::value>{});
330 : }
331 :
332 : /// \ingroup ConstantExpressionsGroup
333 : /// \brief Compute the length of a const char* at compile time
334 1 : SPECTRE_ALWAYS_INLINE constexpr size_t cstring_length(const char* str) {
335 : // clang-tidy: do not use pointer arithmetic
336 : return *str != 0 ? 1 + cstring_length(str + 1) : 0; // NOLINT
337 : }
338 :
339 : /// \ingroup ConstantExpressionsGroup
340 : /// \brief Compute a hash of a const char* at compile time
341 1 : SPECTRE_ALWAYS_INLINE constexpr size_t cstring_hash(const char* str) {
342 : // clang-tidy: do not use pointer arithmetic
343 : return *str != 0
344 : ? (cstring_hash(str + 1) * 33) ^ // NOLINT
345 : static_cast<size_t>(*str)
346 : : 5381;
347 : }
348 :
349 : namespace ConstantExpression_detail {
350 : template <typename T, size_t Size, size_t... I, size_t... J>
351 : inline constexpr std::array<std::decay_t<T>, Size> replace_at_helper(
352 : const std::array<T, Size>& arr, const T& value, const size_t i,
353 : std::integer_sequence<size_t, I...> /*unused*/,
354 : std::integer_sequence<size_t, J...> /*unused*/) {
355 : // clang-tidy: Cannot use gsl::at because we want constexpr evaluation and
356 : // Parallel::abort violates this
357 : return std::array<std::decay_t<T>, Size>{
358 : {arr[I]..., value, arr[i + J]...}}; // NOLINT
359 : }
360 : } // namespace ConstantExpression_detail
361 :
362 : /// \ingroup ConstantExpressionsGroup
363 : /// Replace at compile time the `I`th entry in the array with `value`
364 : template <size_t I, typename T, size_t Size>
365 1 : inline constexpr std::array<std::decay_t<T>, Size> replace_at(
366 : const std::array<T, Size>& arr, T value) {
367 : return ConstantExpression_detail::replace_at_helper(
368 : arr, std::forward<T>(value), I + 1,
369 : std::make_integer_sequence<size_t, I>{},
370 : std::make_integer_sequence<size_t, Size - I - 1>{});
371 : }
372 :
373 : /// \ingroup ConstantExpressionsGroup
374 : /// Check at compile time if two `std::array`s are equal
375 : template <typename T, typename S, size_t size>
376 1 : inline constexpr bool array_equal(const std::array<T, size>& lhs,
377 : const std::array<S, size>& rhs,
378 : const size_t i = 0) {
379 : // clang-tidy: Cannot use gsl::at because we want constexpr evaluation and
380 : // Parallel::abort violates this
381 : return i < size ? (lhs[i] == rhs[i] // NOLINT
382 : and array_equal(lhs, rhs, i + 1))
383 : : true;
384 : }
|
