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             : template <size_t N, typename DataType>
     114           1 : DataType smoothstep(const double lower_edge, const double upper_edge,
     115             :                     const DataType& arg) {
     116             :   ASSERT(lower_edge < upper_edge,
     117             :          "Requires lower_edge < upper_edge, but lower_edge="
     118             :              << lower_edge << " and upper_edge=" << upper_edge);
     119             :   constexpr auto coeffs = []() {
     120             :     std::array<double, 2 * N + 2> result{};
     121             :     for (size_t k = 0; k <= N; ++k) {
     122             :       result[N + 1 + k] = (k % 2 == 0 ? 1. : -1.) * binomial(2 * N + 1, N - k) *
     123             :                           binomial(N + k, k);
     124             :     }
     125             :     return result;
     126             :   }();
     127             :   using std::clamp;
     128             :   const DataType x = clamp(
     129             :       static_cast<DataType>((arg - lower_edge) / (upper_edge - lower_edge)), 0.,
     130             :       1.);
     131             :   return evaluate_polynomial(coeffs, x);
     132             : }
     133             : 
     134             : /*!
     135             :  * \ingroup UtilitiesGroup
     136             :  * \brief Derivative of the $N$-th order smoothstep function
     137             :  *
     138             :  * The general form of the derivative of the smoothstep function is:
     139             :  *
     140             :  * \begin{equation}
     141             :  * \tilde{S}_N'(x) = (2N + 1) \binom{2N}{N} (x - x^2)^N
     142             :  * \end{equation}
     143             :  *
     144             :  * Since the smoothstep function is $C^N$ continuous at the edges, this
     145             :  * function is $C^{N-1}$ continuous at the edges.
     146             :  */
     147             : template <size_t N, typename DataType>
     148           1 : DataType smoothstep_deriv(const double lower_edge, const double upper_edge,
     149             :                           const DataType& arg) {
     150             :   ASSERT(lower_edge < upper_edge,
     151             :          "Requires lower_edge < upper_edge, but lower_edge="
     152             :              << lower_edge << " and upper_edge=" << upper_edge);
     153             :   using std::clamp;
     154             :   const DataType x = clamp(
     155             :       static_cast<DataType>((arg - lower_edge) / (upper_edge - lower_edge)), 0.,
     156             :       1.);
     157             :   constexpr auto coeff = (2 * N + 1) * binomial(2 * N, N);
     158             :   return coeff * pow<N>(x - square(x));
     159             : }
     160             : 
     161             : /// \ingroup UtilitiesGroup
     162             : /// \brief Defines the inverse square-root (\f$1/\sqrt{x}\f$) for arithmetic
     163             : /// and complex types
     164             : template <typename T, Requires<std::is_arithmetic<T>::value or
     165             :                                tt::is_a_v<std::complex, T>> = nullptr>
     166           1 : auto invsqrt(const T& arg) {
     167             :   return static_cast<T>(1.0) / sqrt(arg);
     168             : }
     169             : 
     170             : /// \ingroup UtilitiesGroup
     171             : /// \brief Defines the inverse cube-root (\f$1/\sqrt[3]{x}\f$) for arithmetic
     172             : /// types
     173             : template <typename T, Requires<std::is_arithmetic<T>::value> = nullptr>
     174           1 : auto invcbrt(const T& arg) {
     175             :   return static_cast<T>(1.0) / cbrt(arg);
     176             : }
     177             : 
     178             : namespace sgn_detail {
     179             : template <typename T>
     180             : constexpr T sgn(const T& val, std::true_type /*is_signed*/) {
     181             :   return static_cast<T>(static_cast<T>(0) < val) -
     182             :          static_cast<T>(val < static_cast<T>(0));
     183             : }
     184             : 
     185             : template <typename T>
     186             : constexpr T sgn(const T& val, std::false_type /*is_signed*/) {
     187             :   return static_cast<T>(static_cast<T>(0) < val);
     188             : }
     189             : }  // namespace sgn_detail
     190             : 
     191             : /// \ingroup UtilitiesGroup
     192             : /// \brief Compute the sign function of `val` defined as `1` if `val > 0`, `0`
     193             : /// if `val == 0`, and `-1` if `val < 0`.
     194             : template <typename T>
     195           1 : constexpr T sgn(const T& val) {
     196             :   return sgn_detail::sgn(val, std::is_signed<T>{});
     197             : }
     198             : 
     199             : /// \ingroup UtilitiesGroup
     200             : /// \brief Raises a double to the integer power n.
     201           1 : inline double integer_pow(const double x, const int e) {
     202             :   ASSERT(e >= 0, "Negative powers are not implemented");
     203             :   int ecount = e;
     204             :   int bitcount = 1;
     205             :   while (ecount >>= 1) {
     206             :     ++bitcount;
     207             :   }
     208             :   double result = 1.;
     209             :   while (bitcount) {
     210             :     result *= result;
     211             :     if ((e >> --bitcount) & 0x1) {
     212             :       result *= x;
     213             :     }
     214             :   }
     215             :   return result;
     216             : }