fast_transcendental.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_fast_transcendental_h
#define dealii_fast_transcendental_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/numbers.h>
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/vectorization.h>
 
#include <array>
#include <cstdint>
#include <cstring>
#include <type_traits>
 
#ifdef _MSC_VER
#  include <intrin.h>
#elif defined(__ARM_NEON)
#  include <arm_neon.h>
#elif defined(__x86_64__)
#  include <x86intrin.h>
#endif
 
DEAL_II_NAMESPACE_OPEN
 
namespace fast_transcendental
{
  namespace internal
  {
    /*
     * These arrays contain the polynomial coefficients used to approximate
     * the correction function K(y_f) that appears in the vectorized exponential
     * approximation. Each specialization provides the coefficients for a
     * specific polynomial degree.
     *
     * The coefficients have been determined via a least-squares fit using 1e8
     * uniformly spaced sample points over the interval [0, 1].
     *
     * The coefficients are ordered from the lowest to the highest polynomial
     * term, i.e., exp_pol_coeff[i] corresponds to the coefficient of
     * (y_f)^i.
     *
     * The primary template triggers a static assertion to ensure that an
     * unsupported polynomial degree results in a clear, compile-time error.
     */
    template <int degree, typename Number>
    constexpr std::array<Number, degree + 1> exp_pol_coeff =
      []() -> std::array<Number, degree + 1> {
      static_assert(degree > 2 && degree < 12,
                    "Unsupported polynomial degree for exp_pol_coeff.");
      return std::array<Number, degree + 1>{};
    }();
 
    template <typename Number>
    constexpr std::array<Number, 4> exp_pol_coeff<3, Number> = {
      {1.8809206555322363e-04,
       3.0316116404944521e-01,
       -2.2412622972531135e-01,
       -7.9019886953785187e-02}};
 
    template <typename Number>
    constexpr std::array<Number, 5> exp_pol_coeff<4, Number> = {
      {-7.2868375805796964e-06,
       3.0706874220004554e-01,
       -2.4171033150558635e-01,
       -5.1666839694436362e-02,
       -1.3676523629674565e-02}};
 
    template <typename Number>
    constexpr std::array<Number, 6> exp_pol_coeff<5, Number> = {
      {2.3053723876408911e-07,
       3.0684322094757277e-01,
       -2.4013168272249241e-01,
       -5.5876569798473615e-02,
       -8.9405772567122607e-03,
       -1.8943785491846432e-03}};
 
    template <typename Number>
    constexpr std::array<Number, 7> exp_pol_coeff<6, Number> = {
      {-6.1646222105268035e-09,
       3.0685316242622318e-01,
       -2.4023109751048086e-01,
       -5.5478910644035311e-02,
       -9.6861881733345725e-03,
       -1.2382409419015978e-03,
       -2.1871253576104385e-04}};
 
    template <typename Number>
    constexpr std::array<Number, 8> exp_pol_coeff<7, Number> = {
      {1.4276280848886066e-10,
       3.0685280921264796e-01,
       -2.4022632912712596e-01,
       -5.5505401662914296e-02,
       -9.6133378710895074e-03,
       -1.3431453773612783e-03,
       -1.4294822119902167e-04,
       -2.1646947017882258e-05}};
 
    template <typename Number>
    constexpr std::array<Number, 9> exp_pol_coeff<8, Number> = {
      {-2.9146986289710858e-12,
       3.0685281970141953e-01,
       -2.4022651268062989e-01,
       -5.5504055603868258e-02,
       -9.6183855925465120e-03,
       -1.3326461166925901e-03,
       -1.5519735866800178e-04,
       -1.4147475092332136e-05,
       -1.8748679813608566e-06}};
 
    template <typename Number>
    constexpr std::array<Number, 10> exp_pol_coeff<9, Number> = {
      {5.4337475594757437e-14,
       3.0685281943421089e-01,
       -2.4022650680203564e-01,
       -5.5504110470759059e-02,
       -9.6181181164443925e-03,
       -1.3333950497836748e-03,
       -1.5394913684739714e-04,
       -1.5370223000068715e-05,
       -1.2252831550509094e-06,
       -1.4435218373241001e-07}};
 
    template <typename Number>
    constexpr std::array<Number, 11> exp_pol_coeff<10, Number> = {
      {7.6072080297756739e-12,
       3.0685281862960340e-01,
       -2.4022648563631835e-01,
       -5.5504349821691862e-02,
       -9.6166785909533817e-03,
       -1.3384970969245304e-03,
       -1.4276337493295606e-04,
       -3.0711077326624553e-05,
       1.1584454591775664e-05,
       -6.0984990187293866e-06,
       1.1810109241254786e-06}};
 
    template <typename Number>
    constexpr std::array<Number, 12> exp_pol_coeff<11, Number> = {
      {3.6083394569535254e-13,
       3.0685281939395403e-01,
       -2.4022650549710606e-01,
       -5.5504128716495800e-02,
       -9.6179813320517119e-03,
       -1.3340079120703319e-03,
       -1.5221190054222076e-04,
       -1.8563222423500779e-05,
       2.5694826992914108e-06,
       -2.9577981133179920e-06,
       1.1827521620524028e-06,
       -2.1525063398151332e-07}};
 
    template <typename Number>
    constexpr Number exp_coeff_a = numbers::signaling_nan<Number>();
    template <>
    inline constexpr double exp_coeff_a<double> = 4503599627370496;
    template <>
    inline constexpr float exp_coeff_a<float> = 8388608;
 
    template <typename Number>
    inline constexpr Number exp_coeff_b = numbers::signaling_nan<Number>();
    template <>
    inline constexpr double exp_coeff_b<double> = 4.60718241880001741e+18;
    template <>
    inline constexpr float exp_coeff_b<float> = 1065353216;
 
    template <typename Number>
    inline constexpr Number max_exponent = numbers::signaling_nan<Number>();
    template <>
    inline constexpr double max_exponent<double> = 9.218868437227405312e+18;
    template <>
    inline constexpr float max_exponent<float> = 2139095040;
 
    template <typename Number>
    inline constexpr Number exp_max_abs_x = numbers::signaling_nan<Number>();
    template <>
    inline constexpr double exp_max_abs_x<double> = 709.196208642;
    template <>
    inline constexpr float exp_max_abs_x<float> = 87.49823353f;
 
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 8>
    type_cast(const VectorizedArray<double, 8> &in)
    {
      VectorizedArray<double, 8> out{};
      const __m512i              integer = _mm512_cvtpd_epi64(in.data);
      out.data                           = _mm512_castsi512_pd(integer);
      return out;
    }
 
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 16>
    type_cast(const VectorizedArray<float, 16> &in)
    {
      VectorizedArray<float, 16> out{};
      const __m512i              integer = _mm512_cvtps_epi32(in.data);
      out.data                           = _mm512_castsi512_ps(integer);
      return out;
    }
#endif
 
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 4>
    type_cast(const VectorizedArray<double, 4> &in)
    {
      VectorizedArray<double, 4> out;
      int64_t                    int_values[4];
      for (int i = 0; i < 4; i++)
        {
          int_values[i] = static_cast<int64_t>(in[i]);
        }
      out.data = _mm256_castsi256_pd(
        _mm256_loadu_si256(reinterpret_cast<__m256i *>(int_values)));
      return out;
    }
 
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 8>
    type_cast(const VectorizedArray<float, 8> &in)
    {
      VectorizedArray<float, 8> out;
      int32_t                   int_values[8];
      for (int i = 0; i < 8; i++)
        {
          int_values[i] = static_cast<int32_t>(in[i]);
        }
      out.data = _mm256_castsi256_ps(
        _mm256_loadu_si256(reinterpret_cast<__m256i *>(int_values)));
      return out;
    }
#endif
 
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 2>
    type_cast(const VectorizedArray<double, 2> &in)
    {
      VectorizedArray<double, 2> out;
      int64_t                    int_values[2];
      for (int i = 0; i < 2; i++)
        {
          int_values[i] = static_cast<int64_t>(in[i]);
        }
      out.data = _mm_castsi128_pd(
        _mm_loadu_si128(reinterpret_cast<__m128i *>(int_values)));
      return out;
    }
 
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 4>
    type_cast(const VectorizedArray<float, 4> &in)
    {
      VectorizedArray<float, 4> out;
      int32_t                   int_values[4];
      for (int i = 0; i < 4; i++)
        {
          int_values[i] = static_cast<int32_t>(in[i]);
        }
      out.data = _mm_castsi128_ps(
        _mm_loadu_si128(reinterpret_cast<__m128i *>(int_values)));
      return out;
    }
#endif
 
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 2>
    type_cast(const VectorizedArray<double, 2> &in)
    {
      VectorizedArray<double, 2> out;
      const int64x2_t            int_values = vcvtq_s64_f64(in.data);
      out.data                              = vreinterpretq_f64_s64(int_values);
      return out;
    }
 
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 4>
    type_cast(const VectorizedArray<float, 4> &in)
    {
      VectorizedArray<float, 4> out;
      const int32x4_t           int_values = vcvtq_s32_f32(in.data);
      out.data                             = vreinterpretq_f32_s32(int_values);
      return out;
    }
#endif
 
    template <typename number>
    DEAL_II_ALWAYS_INLINE inline VectorizedArray<number, 1>
    type_cast(const VectorizedArray<number, 1> &in)
    {
      VectorizedArray<number, 1> out;
      const auto                 result = static_cast<
        std::conditional_t<std::is_same_v<number, double>, int64_t, int32_t>>(
        in.data);
      std::memcpy(&out.data, &result, sizeof(out.data));
      return out;
    }
 
    template <int degree, typename CoeffsType, typename ValueType>
    DEAL_II_ALWAYS_INLINE inline ValueType
    horner_scheme(const std::array<CoeffsType, degree + 1> &c,
                  const ValueType                          &x)
    {
      ValueType result = c.back();
      for (int i = c.size() - 2; i >= 0; --i)
        result = c[i] + x * result;
 
      return result;
    }
 
    template <int degree, typename CoeffsType, typename ValueType>
    DEAL_II_ALWAYS_INLINE inline ValueType
    estrin_scheme(const std::array<CoeffsType, degree + 1> &c,
                  const ValueType                          &x)
    {
      constexpr int n_precomputed_x_powers = degree / 2;
      std::array<ValueType, n_precomputed_x_powers> x_powers;
      x_powers[0] = x * x;
      for (unsigned i = 1; i < n_precomputed_x_powers; ++i)
        x_powers[i] = x_powers[i - 1] * x_powers[0];
 
      ValueType t1, t2;
      t1 = (c[0] + c[1] * x);
 
      for (unsigned i = 2; i < degree; i += 2)
        {
          t2 = (c[i] + c[i + 1] * x);
          t1 += t2 * x_powers[(i - 2) / 2];
        }
 
      if constexpr (degree % 2 == 0)
        t1 += c.back() * x_powers.back();
      return t1;
    }
  } // namespace internal
 
  enum class PolynomialEvalScheme : int
  {
    estrin,
    horner
  };
 
 
  template <int degree,
            typename Number,
            bool                 ensure_correct_treatment_of_zero     = true,
            bool                 ensure_correct_treatment_of_infinity = true,
            PolynomialEvalScheme polynomial_evaluation_scheme =
              PolynomialEvalScheme::horner>
  DEAL_II_ALWAYS_INLINE inline Number
  exp(Number x)
  {
    if constexpr (std::is_floating_point_v<Number>)
      {
        auto r = exp<degree,
                     ::VectorizedArray<Number, 1>,
                     ensure_correct_treatment_of_zero,
                     ensure_correct_treatment_of_infinity,
                     polynomial_evaluation_scheme>(
          ::VectorizedArray<Number, 1>(x));
        return r.data;
      }
    else
      {
        using floating_type = typename Number::value_type;
 
        static_assert(
          std::is_same_v<floating_type, double> ||
            std::is_same_v<floating_type, float>,
          "The fast transcendental approximation of the exp() function only supports the fundamental types float and double.");
 
        Number r = x * static_cast<floating_type>(numbers::LOG2E);
 
        Number fractional_exponent = r - r.get_floor();
 
        if constexpr (polynomial_evaluation_scheme ==
                      PolynomialEvalScheme::horner)
          r -= internal::horner_scheme<degree>(
            internal::exp_pol_coeff<degree, floating_type>,
            fractional_exponent);
        else if constexpr (polynomial_evaluation_scheme ==
                           PolynomialEvalScheme::estrin)
          r -= internal::estrin_scheme<degree>(
            internal::exp_pol_coeff<degree, floating_type>,
            fractional_exponent);
        else
          static_assert(
            polynomial_evaluation_scheme == PolynomialEvalScheme::estrin ||
              polynomial_evaluation_scheme == PolynomialEvalScheme::horner,
            "The provided scheme for the evaluation of the correction function is not supported!");
 
        r = internal::exp_coeff_a<floating_type> * r +
            internal::exp_coeff_b<floating_type>;
 
        // To ensure correct handling of zero and infinity, the input to the
        // static_cast operation is adjusted so that values falling outside the
        // valid range produce exactly zero or infinity.
        if constexpr (ensure_correct_treatment_of_zero)
          {
            // Denormal values are excluded since the approximation method does
            // not work for them.
            r = compare_and_apply_mask<SIMDComparison::less_than>(
              x,
              Number(-internal::exp_max_abs_x<floating_type>),
              Number(0.),
              r);
          }
        if constexpr (ensure_correct_treatment_of_infinity)
          {
            r = compare_and_apply_mask<SIMDComparison::greater_than>(
              x,
              Number(internal::exp_max_abs_x<floating_type>),
              Number(internal::max_exponent<floating_type>),
              r);
          }
 
        return internal::type_cast(r);
      }
  }
} // namespace fast_transcendental
 
DEAL_II_NAMESPACE_CLOSE
 
#endif