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