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