Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cmath>
8 : #include <cstddef>
9 : #include <numeric>
10 : #include <type_traits>
11 : #include <vector>
12 :
13 : #include "Utilities/ConstantExpressions.hpp"
14 : #include "Utilities/ErrorHandling/Assert.hpp"
15 : #include "Utilities/ForceInline.hpp"
16 : #include "Utilities/MakeWithValue.hpp"
17 : #include "Utilities/Requires.hpp"
18 : #include "Utilities/TypeTraits.hpp"
19 : #include "Utilities/TypeTraits/IsA.hpp"
20 : #include "Utilities/TypeTraits/IsInteger.hpp"
21 :
22 : // using for overload resolution with blaze
23 : // clang-tidy doesn't want these in the global namespace
24 : using std::conj; // NOLINT
25 : using std::imag; // NOLINT
26 : using std::real; // NOLINT
27 :
28 : /*!
29 : * \ingroup UtilitiesGroup
30 : * \brief Returns the number of digits in an integer number
31 : */
32 : template <typename T>
33 1 : SPECTRE_ALWAYS_INLINE T number_of_digits(const T number) {
34 : static_assert(tt::is_integer_v<std::decay_t<T>>,
35 : "Must call number_of_digits with an integer number");
36 : return number == 0 ? 1
37 : : static_cast<decltype(number)>(
38 : std::ceil(std::log10(std::abs(number) + 1)));
39 : }
40 :
41 : /*!
42 : * \ingroup UtilitiesGroup
43 : * \brief Evaluate a polynomial \f$\sum_{p=0}^N c_p x^p\f$ with Horner's rule
44 : *
45 : * \param coeffs The polynomial coefficients \f$c_p\f$ ordered from constant to
46 : * largest power
47 : * \param x The polynomial variable \f$x\f$
48 : *
49 : * \tparam CoeffsIterable The type of the polynomial coefficients \p coeffs. Can
50 : * be a `std::vector<double>` or `std::array<double>`, which means the
51 : * coefficients are constant for all values in \p x. Each coefficient can also
52 : * be a vector type of typically the same size as \p x, which means the
53 : * coefficients vary with the elements in \p x.
54 : * \tparam DataType The type of the polynomial variable \p x. Must support
55 : * `make_with_value<DataType, DataType>`, as well as (elementwise) addition with
56 : * `CoeffsIterable::value_type` and multiplication with `DataType`.
57 : */
58 : template <typename CoeffsIterable, typename DataType>
59 1 : DataType evaluate_polynomial(const CoeffsIterable& coeffs, const DataType& x) {
60 : return std::accumulate(coeffs.rbegin(), coeffs.rend(),
61 : make_with_value<DataType>(x, 0.),
62 : [&x](const DataType& state, const auto& element) {
63 : return state * x + element;
64 : });
65 : }
66 :
67 : /// \ingroup UtilitiesGroup
68 :
69 : /// \brief Defines the Heaviside step function \f$\Theta\f$ for arithmetic
70 : /// types. \f$\Theta(0) = 1\f$.
71 : template <typename T, Requires<std::is_arithmetic<T>::value> = nullptr>
72 1 : constexpr T step_function(const T& arg) {
73 : return static_cast<T>((arg >= static_cast<T>(0)) ? 1 : 0);
74 : }
75 :
76 : /*!
77 : * \ingroup UtilitiesGroup
78 : * \brief Smoothly interpolates from 0 to 1 between `lower_edge` and
79 : * `upper_edge` with a Hermite polynomial of degree `2 * N + 1`.
80 : *
81 : * The smoothstep function is
82 : *
83 : * \f{align*}
84 : * S_N(x) = \begin{cases}
85 : * 0 &\quad \text{for} \quad x\leq x_0 \\
86 : * \tilde{S}_N((x - x_0) / (x_1 - x_0))
87 : * &\quad \text{for} \quad x_0 \leq x\leq x_1 \\
88 : * 1 &\quad \text{for} \quad x_1\leq x \\
89 : * \end{cases}
90 : * \f}
91 : *
92 : * where \f$x_0\f$ is `lower_edge`, \f$x_1\f$ is `upper_edge`, and, up to
93 : * \f$N=3\f$,
94 : *
95 : * \f{align*}
96 : * \tilde{S}_0(x) &= x \\
97 : * \tilde{S}_1(x) &= 3x^2 - 2x^3 \\
98 : * \tilde{S}_2(x) &= 10x^3 - 15x^4 + 6x^5 \\
99 : * \tilde{S}_3(x) &= 35x^4 - 84x^5 + 70x^6 - 20x^7
100 : * \text{.}
101 : * \f}
102 : */
103 : template <size_t N, typename DataType>
104 1 : DataType smoothstep(const double lower_edge, const double upper_edge,
105 : const DataType& arg) {
106 : ASSERT(lower_edge < upper_edge,
107 : "Requires lower_edge < upper_edge, but lower_edge="
108 : << lower_edge << " and upper_edge=" << upper_edge);
109 : using std::clamp;
110 : return evaluate_polynomial(
111 : []() -> std::array<double, 2 * N + 2> {
112 : static_assert(N <= 3,
113 : "The smoothstep function is currently only implemented "
114 : "for N <= 3.");
115 : if constexpr (N == 0) {
116 : return {0., 1};
117 : } else if constexpr (N == 1) {
118 : return {0., 0., 3., -2};
119 : } else if constexpr (N == 2) {
120 : return {0., 0., 0., 10., -15., 6.};
121 : } else if constexpr (N == 3) {
122 : return {0., 0., 0., 0., 35., -84., 70., -20.};
123 : }
124 : }(),
125 : static_cast<DataType>(clamp(
126 : static_cast<DataType>((arg - lower_edge) / (upper_edge - lower_edge)),
127 : 0., 1.)));
128 : }
129 :
130 : /// \ingroup UtilitiesGroup
131 : /// \brief Defines the inverse square-root (\f$1/\sqrt{x}\f$) for arithmetic
132 : /// and complex types
133 : template <typename T, Requires<std::is_arithmetic<T>::value or
134 : tt::is_a_v<std::complex, T>> = nullptr>
135 1 : auto invsqrt(const T& arg) {
136 : return static_cast<T>(1.0) / sqrt(arg);
137 : }
138 :
139 : /// \ingroup UtilitiesGroup
140 : /// \brief Defines the inverse cube-root (\f$1/\sqrt[3]{x}\f$) for arithmetic
141 : /// types
142 : template <typename T, Requires<std::is_arithmetic<T>::value> = nullptr>
143 1 : auto invcbrt(const T& arg) {
144 : return static_cast<T>(1.0) / cbrt(arg);
145 : }
146 :
147 : namespace sgn_detail {
148 : template <typename T>
149 : constexpr T sgn(const T& val, std::true_type /*is_signed*/) {
150 : return static_cast<T>(static_cast<T>(0) < val) -
151 : static_cast<T>(val < static_cast<T>(0));
152 : }
153 :
154 : template <typename T>
155 : constexpr T sgn(const T& val, std::false_type /*is_signed*/) {
156 : return static_cast<T>(static_cast<T>(0) < val);
157 : }
158 : } // namespace sgn_detail
159 :
160 : /// \ingroup UtilitiesGroup
161 : /// \brief Compute the sign function of `val` defined as `1` if `val > 0`, `0`
162 : /// if `val == 0`, and `-1` if `val < 0`.
163 : template <typename T>
164 1 : constexpr T sgn(const T& val) {
165 : return sgn_detail::sgn(val, std::is_signed<T>{});
166 : }
167 :
168 : /// \ingroup UtilitiesGroup
169 : /// \brief Raises a double to the integer power n.
170 1 : inline double integer_pow(const double x, const int e) {
171 : ASSERT(e >= 0, "Negative powers are not implemented");
172 : int ecount = e;
173 : int bitcount = 1;
174 : while (ecount >>= 1) {
175 : ++bitcount;
176 : }
177 : double result = 1.;
178 : while (bitcount) {
179 : result *= result;
180 : if ((e >> --bitcount) & 0x1) {
181 : result *= x;
182 : }
183 : }
184 : return result;
185 : }
|