libc/src/__support/math/exp2f_float_utils.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
/// Shared utilities and constants for single-precision exponential functions.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_FLOAT_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_FLOAT_UTILS_H

#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace float_eval {

// Compute y * 2^m for y in [0.5, 2.0) and m in [-150, 128].
//
// We divide into three cases depending on m:
// 1. -126 <= m <= 127:
//    2^m is a normal float (exponent field m + 127 in [1, 254]).
//    We can directly compute y * 2^m with a single multiplication.
//    The condition is checked using unsigned comparison: (m + 126) <= 253.
//
// 2. m > 127 (i.e., m = 128):
//    2^m cannot be represented as a finite normal float (exponent field
//    128 + 127 = 255 is inf). Instead, we compute:
//      y * 2^m = (2 * y) * 2^(m - 1).
//    Since y < 2.0, 2 * y is exact (simply increments the exponent of y),
//    and 2^(m - 1) = 2^127 is a normal float.
//
// 3. m < -126 (i.e., m in [-150, -127]):
//    The result is in the subnormal range (or underflows to 0).
//    To avoid double-rounding errors, we scale in two steps:
//      y * 2^m = (y * 2^(m + 25)) * 2^-25.
//    Since m >= -150, m + 25 >= -125, so 2^(m + 25) is a normal float.
//    The intermediate product (y * 2^(m + 25)) is normal and computed with
//    full precision. Then the final multiplication by 2^-25 rounds directly
//    to the target subnormal precision in a single rounding step.
LIBC_INLINE float scale_exp2f(float y, int m) {
  // -126 <= m <= 127: 2^m is a normal float.
  if (LIBC_LIKELY(static_cast<uint32_t>(m + 126) <= 253)) {
    uint32_t scale_bits = static_cast<uint32_t>(m + 127) << 23;
    float scale = cpp::bit_cast<float>(scale_bits);
    return y * scale;
  }

  // m = 128: 2^m cannot be represented as a normal float.
  // Compute (2 * y) * 2^(m - 1), where 2 * y is exact and 2^(m - 1) is normal.
  if (m > 127) {
    uint32_t scale_bits = static_cast<uint32_t>(m - 1 + 127) << 23;
    float scale = cpp::bit_cast<float>(scale_bits);
    return (y * 2.0f) * scale;
  }

  // Subnormal result: m in [-150, -127].
  // Compute (y * 2^(m + 25)) * 2^-25 so that the first product is normal and
  // the second multiplication rounds directly to the subnormal result.
  uint32_t scale_bits = static_cast<uint32_t>(m + 25 + 127) << 23;
  float scale = cpp::bit_cast<float>(scale_bits);
  return (y * scale) * 0x1.0p-25f;
}

// Degree-6 polynomial approximating 2^u - 1 = u * P(u) on [-0.5, 0.5] and
// exponent scaling to evaluate 2^x = 2^m * 2^u.
// Generated by Sollya with:
// > display = hexadecimal;
// > P = fpminimax((2^x - 1)/x, 5, [|SG...|], [-0.5, 0.5]);
// > dirtyinfnorm(1 + x * P - 2^x, [-0.5, 0.5]);
//   0x1.68cb...p-28
// > dirtyinfnorm((1 + x * P - 2^x) / 2^x, [-0.5, 0.5]);
//   0x1.5109...p-28
//
// Instead of using pure Estrin's scheme on P(u), we evaluate the degree-4
// part p = c1 + c2 * u + c3 * u^2 + c4 * u^3 + c5 * u^4 with Estrin's
// scheme, and then evaluate y = 1.0f + u * (c0 + u * p) with Horner's
// scheme.
// In pure Estrin's scheme, p0 = c0 + c1 * u, q0 = p0 + u^2 * p1, and
// p = q0 + u^4 * p2 all round at the exponent of c0 (~2^-1), accumulating
// rounding errors on non-FMA targets. Evaluating c1..c5 first keeps only a
// single rounding at c0's scale.
LIBC_INLINE float exp2f_eval(float u, int m) {
  constexpr float COEFFS[6] = {0x1.62e430p-1f,  0x1.ebfbdep-3f,
                               0x1.c6af9ep-5f,  0x1.3b2c58p-7f,
                               0x1.5ef4a2p-10f, 0x1.427918p-13f};

  float u2 = u * u;
  float u4 = u2 * u2;

  float p0 = fputil::multiply_add(u, COEFFS[2], COEFFS[1]);
  float p1 = fputil::multiply_add(u, COEFFS[4], COEFFS[3]);
  float q0 = fputil::multiply_add(u2, p1, p0);
  float p = fputil::multiply_add(u4, COEFFS[5], q0);

  float y = fputil::polyeval(u, 1.0f, COEFFS[0], p);

  return scale_exp2f(y, m);
}

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

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_FLOAT_UTILS_H