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             : }