libc/src/__support/math/exp10f_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 exp10f.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_FLOAT_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_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 exp10f(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| >= log10(2^128), or x is nan
  if (LIBC_UNLIKELY(x_abs >= 0x421a'209bU)) {
    // When x < log10(2^-150) or nan
    if (x_u > 0xc234'9e35U) {
      // 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 >= log10(2^128) or nan
    if (xbits.is_pos() && (x_u >= 0x421a'209bU)) {
      // x is finite
      if (x_u < 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();
    }
  }

  // |x| <= 2^-25
  // 10^x ~ 1 + log(10) * x
  if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
    return fputil::multiply_add(x, 0x1.26bb1cp+1f, 1.0f);
  }

  // Exact outputs when x = 1, 2, ..., 10.
  // Quick check mask: 0x800f'ffffU = ~(bits of 1.0f | ... | bits of 10.0f)
  if (LIBC_UNLIKELY((x_u & 0x800f'ffffU) == 0)) {
    switch (x_u) {
    case 0x3f800000U: // x = 1.0f
      return 10.0f;
    case 0x40000000U: // x = 2.0f
      return 100.0f;
    case 0x40400000U: // x = 3.0f
      return 1'000.0f;
    case 0x40800000U: // x = 4.0f
      return 10'000.0f;
    case 0x40a00000U: // x = 5.0f
      return 100'000.0f;
    case 0x40c00000U: // x = 6.0f
      return 1'000'000.0f;
    case 0x40e00000U: // x = 7.0f
      return 10'000'000.0f;
    case 0x41000000U: // x = 8.0f
      return 100'000'000.0f;
    case 0x41100000U: // x = 9.0f
      return 1'000'000'000.0f;
    case 0x41200000U: // x = 10.0f
      return 10'000'000'000.0f;
    }
  }

  // Range reduction:
  //   k = round(x * log2(10))
  //   x * log2(10) = k + u, with |u| <= 0.5
  //   10^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(10), SG, RN);
  // > lo = round(log2(10) - hi, SG, RN);
  constexpr fputil::FloatFloat LOG2_10 = {0x1.2f346ep-24f, 0x1.a934fp+1f};

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

  float u_hi = fputil::multiply_add(x, LOG2_10.hi, -kf);
  float u = fputil::multiply_add(x, LOG2_10.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_10 = round(log2(10), SG, RN);
  // > LOG10_2_HI = round(log10(2), 12, RN);
  // > LOG10_2_LO = round(log10(2) - LOG10_2_HI, SG, RN);
  constexpr float LOG2_10 = 0x1.a934fp+1f;
  constexpr float LOG10_2_HI = 0x1.344p-2f;
  constexpr float LOG10_2_LO = 0x1.3509f8p-18f;

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

  float v_hi = fputil::multiply_add(-kf, LOG10_2_HI, x);
  float v = fputil::multiply_add(-kf, LOG10_2_LO, v_hi);

  // Convert reduced argument to base-2: u = v * log2(10)
  float u = v * LOG2_10;
#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_EXP10F_FLOAT_EVAL_H