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