libc/src/__support/math/expf_integer_eval.h

//===----------------------------------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
///
/// \file
/// Integer-only implementation of expf.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_INTEGER_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_INTEGER_EVAL_H

#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/frac64.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/math/check/exp_exceptions.h"
#include "src/__support/math/exp2f_integer_utils.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace static_rounding {

// Constants generated by Sollya with:
// > display = hexadecimal;
// > round(log2(exp(1)) * 2^63, 64, RN);
// LSB(LOG2_E) = 2^-63
LIBC_INLINE_VAR constexpr Frac64 LOG2_E = Frac64(0xb8aa'3b29'5c17'f0bc);

LIBC_ALWAYS_INLINE uint32_t expf_bits(uint32_t x_u,
                                      [[maybe_unused]] int rounding) {
  using FPBits = fputil::FPBits<float>;
  using FPBounds = math::check::exp_internal::Bounds<float>;

  bool is_neg = (x_u >> 31) != 0;
  uint32_t x_abs = x_u & 0x7fff'ffffU;
  uint32_t x_e = x_abs >> FPBits::FRACTION_LEN;

  // When |x| < 2^-25 (x_e <= 101) or |x| >= 64 (x_e >= 133):
  if (LIBC_UNLIKELY(x_e - 102U >= 31U)) {
    // |x| < 2^-25
    // exp(x) = 1 + 2^-23   if x > 0 and FE_UPWARD
    //        = 1 - 2^-24   if x < 0 and FE_DOWNWARD or FE_TOWARDZERO
    //        = 1           otherwise.
    if (x_e <= 101U) {
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      if (x_abs != 0) {
        if (rounding == FE_UPWARD && !is_neg)
          return 0x3f80'0001U; // 0x1.000002p+0f = 1.0f + 2^-23f

        if ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) && is_neg)
          return 0x3f7f'ffffU; // 0x1.fffffep-1f = 1.0f - 2^-24f
      }
#endif                     // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      return 0x3f80'0000U; // 1.0f
    }

    // exp(nan) = nan
    if (x_abs > 0x7f80'0000U)    // |x| > +inf (x is NaN)
      return x_u | 0x0040'0000U; // silence signaling NaN

    // x >= log(2^127 * (2 - 2^-23)):
    // exp(x) = +inf     for FE_TONEAREST, FE_UPWARD
    //        = FLT_MAX  for FE_TOWARDZERO, FE_DOWNWARD
    if (!is_neg && x_u >= FPBounds::UPPER_BITS) {
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      if ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) &&
          x_u < 0x7f80'0000U) // x < +inf
        return 0x7f7f'ffffU;  // 0x1.fffffep+127f (FLT_MAX)
#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      return 0x7f80'0000U; // +inf
    }

    // x <= -0x1.9fe36ap+6f < log(2^-150)
    // exp(x) = 2^-149    for FE_UPWARD
    //        = 0         otherwise
    if (x_u >= 0xc2cf'f1b5U) {
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      if (rounding == FE_UPWARD && x_u < 0xff80'0000U) // x > -inf
        return 0x0000'0001U; // 0x1.0p-149f (FLT_TRUE_MIN)
#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      return 0x0000'0000U; // 0.0f
    }
  }

  // Range reduction:
  //   k = floor(x * log2(e))
  //   x * log2(e) = k + u, with u in [0, 1)
  //   e^x = 2^(k + u) = 2^k * 2^u
  uint64_t u_bits = exp_b_scale(x_u, x_e, LOG2_E, 133U);
  return exp2f_eval_bits(u_bits, is_neg, rounding);
}

// Statically rounded, no-except implementation of expf using integer-only
// arithmetic.
LIBC_ALWAYS_INLINE float expf(float x, [[maybe_unused]] int rounding) {
  return cpp::bit_cast<float>(expf_bits(cpp::bit_cast<uint32_t>(x), rounding));
}

} // namespace static_rounding

namespace integer_eval {

LIBC_ALWAYS_INLINE float expf(float x) {
  return static_rounding::expf(x, FE_TONEAREST);
}

} // namespace integer_eval
} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_INTEGER_EVAL_H