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 implementation 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/exp2f_float_utils.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace float_eval {

LIBC_INLINE float expf(float x) {
  using FPBits = fputil::FPBits<float>;
  FPBits xbits(x);

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

  // When |x| >= 89, |x| < 2^-25, or x is nan
  if (LIBC_UNLIKELY(x_abs >= 0x42b2'0000U || 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 >= 89 or nan
    if (xbits.is_pos() && (xbits.uintval() >= 0x42b2'0000)) {
      // 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), 12, 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.62ep-1f;
  constexpr float LOG_2_LO = 0x1.0bfbe8p-15f;

  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);
}

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

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_EVAL_H