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 Hermite interpolation of polynomial degree `2 * N + 1`.
80 : *
81 : * The smoothstep function is
82 : *
83 : * \begin{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 : * \end{align*}
91 : *
92 : * where \f$x_0\f$ is `lower_edge` and \f$x_1\f$ is `upper_edge`. The general
93 : * form of the polynomial \f$\tilde{S}_N(x)\f$ is:
94 : *
95 : * \begin{equation}
96 : * \tilde{S}_N(x) = x^{N+1} \sum_{k=0}^{N} \binom{N+k}{k}
97 : * \binom{2N+1}{N-k} (-x)^k
98 : * \end{equation}
99 : *
100 : * The first few polynomials are:
101 : *
102 : * \begin{align*}
103 : * \tilde{S}_0(x) &= x \\
104 : * \tilde{S}_1(x) &= 3x^2 - 2x^3 \\
105 : * \tilde{S}_2(x) &= 10x^3 - 15x^4 + 6x^5 \\
106 : * \tilde{S}_3(x) &= 35x^4 - 84x^5 + 70x^6 - 20x^7
107 : * \text{.}
108 : * \end{align*}
109 : *
110 : * This function is $C^N$ continuous at the edges, i.e. the first $N$
111 : * derivatives are continuous at the edges.
112 : *
113 : * If `lower_edge` and `upper_edge` are equal, this function reduces to the
114 : * step function $\Theta(x - x_0)$.
115 : */
116 : template <size_t N, typename DataType>
117 1 : DataType smoothstep(const double lower_edge, const double upper_edge,
118 : const DataType& arg) {
119 : if (lower_edge == upper_edge) {
120 : return step_function(arg - lower_edge);
121 : }
122 : ASSERT(lower_edge < upper_edge,
123 : "Requires lower_edge < upper_edge, but lower_edge="
124 : << lower_edge << " and upper_edge=" << upper_edge);
125 : constexpr auto coeffs = []() {
126 : std::array<double, 2 * N + 2> result{};
127 : for (size_t k = 0; k <= N; ++k) {
128 : result[N + 1 + k] = (k % 2 == 0 ? 1. : -1.) * binomial(2 * N + 1, N - k) *
129 : binomial(N + k, k);
130 : }
131 : return result;
132 : }();
133 : using std::clamp;
134 : const DataType x = clamp(
135 : static_cast<DataType>((arg - lower_edge) / (upper_edge - lower_edge)), 0.,
136 : 1.);
137 : return evaluate_polynomial(coeffs, x);
138 : }
139 :
140 : /*!
141 : * \ingroup UtilitiesGroup
142 : * \brief Derivative of the $N$-th order smoothstep function
143 : *
144 : * The general form of the derivative of the smoothstep function is:
145 : *
146 : * \begin{equation}
147 : * \tilde{S}_N'(x) = (2N + 1) \binom{2N}{N} (x - x^2)^N
148 : * \end{equation}
149 : *
150 : * Since the smoothstep function is $C^N$ continuous at the edges, this
151 : * function is $C^{N-1}$ continuous at the edges.
152 : *
153 : * If `lower_edge` and `upper_edge` are equal, this function reduces to the
154 : * derivative of the step function, which is zero everywhere except at the edge
155 : * where it is not defined. We define it to be zero in this case as well.
156 : */
157 : template <size_t N, typename DataType>
158 1 : DataType smoothstep_deriv(const double lower_edge, const double upper_edge,
159 : const DataType& arg) {
160 : if (lower_edge == upper_edge) {
161 : return make_with_value<DataType>(arg, 0.);
162 : }
163 : ASSERT(lower_edge < upper_edge,
164 : "Requires lower_edge < upper_edge, but lower_edge="
165 : << lower_edge << " and upper_edge=" << upper_edge);
166 : using std::clamp;
167 : const DataType x = clamp(
168 : static_cast<DataType>((arg - lower_edge) / (upper_edge - lower_edge)), 0.,
169 : 1.);
170 : constexpr auto coeff = (2 * N + 1) * binomial(2 * N, N);
171 : return coeff * pow<N>(x - square(x)) / (upper_edge - lower_edge);
172 : }
173 :
174 : /// \ingroup UtilitiesGroup
175 : /// \brief Defines the inverse square-root (\f$1/\sqrt{x}\f$) for arithmetic
176 : /// and complex types
177 : template <typename T, Requires<std::is_arithmetic<T>::value or
178 : tt::is_a_v<std::complex, T>> = nullptr>
179 1 : auto invsqrt(const T& arg) {
180 : return static_cast<T>(1.0) / sqrt(arg);
181 : }
182 :
183 : /// \ingroup UtilitiesGroup
184 : /// \brief Defines the inverse cube-root (\f$1/\sqrt[3]{x}\f$) for arithmetic
185 : /// types
186 : template <typename T, Requires<std::is_arithmetic<T>::value> = nullptr>
187 1 : auto invcbrt(const T& arg) {
188 : return static_cast<T>(1.0) / cbrt(arg);
189 : }
190 :
191 : namespace sgn_detail {
192 : template <typename T>
193 : constexpr T sgn(const T& val, std::true_type /*is_signed*/) {
194 : return static_cast<T>(static_cast<T>(0) < val) -
195 : static_cast<T>(val < static_cast<T>(0));
196 : }
197 :
198 : template <typename T>
199 : constexpr T sgn(const T& val, std::false_type /*is_signed*/) {
200 : return static_cast<T>(static_cast<T>(0) < val);
201 : }
202 : } // namespace sgn_detail
203 :
204 : /// \ingroup UtilitiesGroup
205 : /// \brief Compute the sign function of `val` defined as `1` if `val > 0`, `0`
206 : /// if `val == 0`, and `-1` if `val < 0`.
207 : template <typename T>
208 1 : constexpr T sgn(const T& val) {
209 : return sgn_detail::sgn(val, std::is_signed<T>{});
210 : }
211 :
212 : /// \ingroup UtilitiesGroup
213 : /// \brief Raises a double to the integer power n.
214 1 : inline double integer_pow(const double x, const int e) {
215 : ASSERT(e >= 0, "Negative powers are not implemented");
216 : int ecount = e;
217 : int bitcount = 1;
218 : while (ecount >>= 1) {
219 : ++bitcount;
220 : }
221 : double result = 1.;
222 : while (bitcount) {
223 : result *= result;
224 : if ((e >> --bitcount) & 0x1) {
225 : result *= x;
226 : }
227 : }
228 : return result;
229 : }
|
