libc/src/__support/math/expf_float_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
/// Float-only implementations of expf.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_EVAL_H

#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/macros/properties/cpu_features.h"
#include "src/__support/math/check/exp_exceptions.h"
#include "src/__support/math/exp2f_float_utils.h"
#include "src/__support/math/expf_float_utils.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace float_eval {

LIBC_INLINE float expf(float x) {
  using FPBits = fputil::FPBits<float>;
  using FPBounds = math::check::exp_internal::Bounds<float>;

  FPBits xbits(x);

  uint32_t x_u = xbits.uintval();
  uint32_t x_abs = x_u & 0x7fff'ffffU;

  // When |x| >= log(2^127 * (2 - 2^-23)), |x| < 2^-25, or x is nan
  if (LIBC_UNLIKELY(x_abs >= FPBounds::UPPER_BITS || x_abs <= 0x3280'0000U)) {
    // |x| < 2^-25
    if (xbits.get_biased_exponent() <= 101) {
      return 1.0f + x;
    }

    // When x < log(2^-150) or nan
    if (xbits.uintval() >= 0xc2cf'f1b5U) {
      // exp(-Inf) = 0
      if (xbits.is_inf())
        return 0.0f;
      // exp(nan) = nan
      if (xbits.is_nan())
        return x;
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      if (fputil::fenv_is_round_up())
        return FPBits::min_subnormal().get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
      fputil::set_errno_if_required(ERANGE);
      fputil::raise_except_if_required(FE_UNDERFLOW);
      return 0.0f;
    }
    // x >= log(2^127 * (2 - 2^-23)) or nan
    if (xbits.is_pos() && (xbits.uintval() >= FPBounds::UPPER_BITS)) {
      // x is finite
      if (xbits.uintval() < 0x7f80'0000U) {
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
        int rounding = fputil::quick_get_round();
        if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
          return FPBits::max_normal().get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY

        fputil::set_errno_if_required(ERANGE);
        fputil::raise_except_if_required(FE_OVERFLOW);
      }
      // x is +inf or nan
      return x + FPBits::inf().get_val();
    }
  }

  // Range reduction:
  //   k = round(x * log2(e))
  //   x * log2(e) = k + u, with |u| <= 0.5
  //   e^x = 2^(k + u) = 2^k * 2^u
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
  // Constants generated by Sollya with:
  // > display = hexadecimal;
  // > hi = round(log2(exp(1)), SG, RN);
  // > lo = round(log2(exp(1)) - hi, SG, RN);
  constexpr fputil::FloatFloat LOG2_E = {/*lo=*/0x1.4ae0cp-26f,
                                         /*hi=*/0x1.715476p+0f};

  float kf = fputil::nearest_integer(x * LOG2_E.hi);
  int k = static_cast<int>(kf);

  float u_hi = fputil::multiply_add(x, LOG2_E.hi, -kf);
  float u = fputil::multiply_add(x, LOG2_E.lo, u_hi);
#else  // !LIBC_TARGET_CPU_HAS_FMA_FLOAT
  // Cody-Waite reduction for non-FMA targets:
  // Constants generated by Sollya with:
  // > display = hexadecimal;
  // > LOG2_E = round(log2(exp(1)), SG, RN);
  // > LOG_2_HI = round(log(2), 16, RN);
  // > LOG_2_LO = round(log(2) - LOG_2_HI, SG, RN);
  constexpr float LOG2_E = 0x1.715476p+0f;
  constexpr float LOG_2_HI = 0x1.62e4p-1f;
  constexpr float LOG_2_LO = 0x1.7f7d1cp-20f;

  float kf = fputil::nearest_integer(x * LOG2_E);
  int k = static_cast<int>(kf);

  float v_hi = fputil::multiply_add(-kf, LOG_2_HI, x);
  float v = fputil::multiply_add(-kf, LOG_2_LO, v_hi);

  // Convert reduced argument to base-2: u = v * log2(e)
  float u = v * LOG2_E;
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT

  return exp2f_eval(u, k);
}

// Designed for GPU performance. Correctly rounded only in round-to-nearest,
// ties-to-even mode. Requires single-precision FMA. Does not support errno
// or floating-point exception reporting.
LIBC_INLINE float cr_expf_rn(float x) {
  using FPBits = fputil::FPBits<float>;
  using FPBounds = math::check::exp_internal::Bounds<float>;

  FPBits x_bits(x);

  uint32_t x_u = x_bits.uintval();
  uint32_t x_abs_u = x_u & 0x7fff'ffffU;

  // |x| <= 2^-25
  bool is_tiny = x_abs_u <= 0x3300'0000U;

  // x >= log(2^127 * (2 - 2^-23)) or x is nan
  bool should_overflow = x_u >= FPBounds::UPPER_BITS && x_u < 0x8000'0000U;

  // x <= log(2^-150) or x is nan
  bool should_underflow = x_u >= 0xc2cf'f1b5U;

  // x is nan with its sign bit set
  bool is_neg_nan = should_underflow && x_bits.is_nan();

  // log(2^-150) < x < log(2^-126)
  bool produces_subnormal = x_u >= FPBounds::LOWER_BITS && !should_underflow;

  // Use x_safe for range reduction to keep the conversion to int defined
  float x_safe = (should_underflow || should_overflow) ? 0.0f : x;

  // Range reduction:
  //   k = round(x * log2(e))
  //   x * log2(e) = k + u * log2(e)
  //   e^x = 2^k * e^u
  //
  // Constants generated by Sollya with:
  //   > display = hexadecimal!;
  //   > LOG2_E = round(log2(exp(1)), SG, RN);
  //   > LOG_2_HI = round(log(2), SG, RN);
  //   > LOG_2_LO.hi = round(log(2) - LOG_2_HI, SG, RN);
  //   > LOG_2_LO.lo = round(log(2) - LOG_2_HI - LOG_2_LO.hi, SG, RN);
  constexpr float LOG2_E = 0x1.715476p+0f;
  constexpr float LOG_2_HI = 0x1.62e43p-1f;
  constexpr fputil::FloatFloat LOG_2_LO = {/*lo=*/-0x1.950d88p-54f,
                                           /*hi=*/-0x1.05c61p-29f};

  float kf = fputil::nearest_integer(x_safe * LOG2_E);
  int k = static_cast<int>(kf);

  float u_hi = fputil::multiply_add(-kf, LOG_2_HI, x_safe);
  fputil::FloatFloat u_lo = fputil::quick_mult(-kf, LOG_2_LO);

  fputil::FloatFloat u = fputil::exact_add(u_hi, u_lo.hi);
  u.lo += u_lo.lo;

  // Polynomial approximation: e^u
  fputil::FloatFloat exp_u = cr_expf_rn_internal::exp_eval(u);

  // Reconstruction: 2^k * e^u
  float res = cr_expf_rn_internal::scale_result(exp_u, k, produces_subnormal);

  // Special cases
  constexpr float INF = FPBits::inf().get_val();

  res = is_tiny ? 1.0f : res;
  res = should_overflow ? x + INF : res;
  res = should_underflow ? (is_neg_nan ? x : 0.0f) : res;

  return res;
}

} // namespace float_eval
} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_EVAL_H