libc/src/__support/math/expf_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
/// Utilities for float-only implementations of expf.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_UTILS_H

#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace float_eval {

namespace cr_expf_rn_internal {

using fputil::FloatFloat;

LIBC_ALWAYS_INLINE FloatFloat horner_step(float x, const FloatFloat &poly,
                                          const FloatFloat &coeff) {
  float hi = fputil::multiply_add(x, poly.hi, coeff.hi);
  float lo = fputil::multiply_add(x, poly.hi, coeff.hi - hi);
  lo = fputil::multiply_add(x, poly.lo, lo + coeff.lo);
  return {lo, hi};
}

// Approximates exp(u) on [-log(2)/2, log(2)/2] with the degree-10
// polynomial P(u) = 1 + u + u^2 * Q(u) using Horner evaluation.
//
// Coefficients generated by Sollya with:
//   > prec = 165!;
//   > display = hexadecimal!;
//   > p = fpminimax(exp(x), [|2, ..., 10|],
//                  [|49, 49, 49, 49, 49, 24, 24, 24, 24|],
//                  [-log(2)/2; log(2)/2], relative,
//                  floating, 1 + x);
//
// Float-float coefficients are split into:
//   > hi = round(coeff(p, i), SG, RN);
//   > lo = round(coeff(p, i) - hi, SG, RN);
LIBC_INLINE FloatFloat exp_eval(const FloatFloat &u) {
  // Coefficients of u^2, ..., u^6, in {lo, hi} order.
  constexpr FloatFloat COEFFS[] = {
      {-0x1.dp-46f, 0x1p-1f},
      {-0x1.555edcp-28f, 0x1.555556p-3f},
      {-0x1.54c28ep-30f, 0x1.555556p-5f},
      {-0x1.b504e8p-32f, 0x1.111112p-7f},
      {0x1.158eb2p-35f, 0x1.6c16cp-10f},
  };

  // The float-only tail contributes to degrees 7 through 10 of P.
  float tail = fputil::polyeval(u.hi, 0x1.a0198ep-13f, 0x1.a01c0ap-16f,
                                0x1.72e878p-19f, 0x1.2694cap-22f);
  FloatFloat poly{0.0f, tail};
  for (int i = 4; i >= 0; --i)
    poly = horner_step(u.hi, poly, COEFFS[i]);

  // Form expm1(u.hi) before adding 1, retaining its low component.
  poly = horner_step(u.hi, poly, FloatFloat{0.0f, 1.0f});
  FloatFloat expm1_u_hi = fputil::quick_mult(u.hi, poly);
  FloatFloat exp_u = fputil::exact_add(1.0f, expm1_u_hi.hi);

  // exp(u.hi + u.lo) = exp(u.hi) + u.lo * exp(u.hi) + O(u.lo^2).
  //
  // Combine the retained low part of expm1(u.hi) with the correction
  // u.lo * exp(u.hi), approximating exp(u.hi) by exp_u.hi.
  float corr = fputil::multiply_add(u.lo, exp_u.hi, expm1_u_hi.lo);

  // Rounding low sums, exp_u.lo + corr, to nearest can create a midpoint for
  // exp_u.hi + RN(exp_u.lo + corr).
  //
  // Round inexact low sums to odd whenever their rounded value is a power of
  // two. This covers all such midpoint cases in this evaluation and preserves
  // which side of the midpoint the exact sum lies on.
  //
  // FastTwoSum here is exact in RN: for nonzero corr, exp_u.lo is an integer
  // multiple of ulp(corr), as required by Theorem 1 of "FastTwoSum revisited".
  FloatFloat lo_sum = fputil::exact_add(exp_u.lo, corr);

  uint32_t lo_sum_u = fputil::FPBits<float>(lo_sum.hi).uintval();
  bool should_round_to_odd =
      (lo_sum_u & fputil::FPBits<float>::FRACTION_MASK) == 0 &&
      lo_sum.lo != 0.0f;
  uint32_t rounding_adj =
      ((lo_sum.hi > 0.0f) == (lo_sum.lo > 0.0f)) ? 1U : uint32_t(-1);
  lo_sum_u += should_round_to_odd ? rounding_adj : 0U;
  exp_u.lo = fputil::FPBits<float>(lo_sum_u).get_val();

  return exp_u;
}

// Computes 2^k * (exp_u.hi + exp_u.lo).
LIBC_INLINE float scale_result(const FloatFloat &exp_u, int k,
                               bool produces_subnormal) {
  using FPBits = fputil::FPBits<float>;

  // Scaling the rounded sum by 2^k is exact for normal results: only the
  // exponent changes. For subnormal results, the fixed spacing 2^-149 is
  // coarser than the scaled spacing of the rounded sum. Scaling can then
  // discard significand bits and round again.
  //
  // Choose B = 2^(-126 - k), so ulp(B) * 2^k = 2^-149. For inputs producing
  // subnormals, B + exp_u.hi + exp_u.lo lies in B's binade [B, 2B). Rounding
  // this biased sum therefore uses the desired spacing while staying normal.
  // Use B = 0 for normal results to share the same arithmetic path.
  uint32_t bias_u = produces_subnormal ? (1U - static_cast<uint32_t>(k))
                                             << FPBits::FRACTION_LEN
                                       : 0U;
  float bias = FPBits(bias_u).get_val();

  // For inputs producing subnormals, B >= exp_u.hi, so FastTwoSum retains
  // the rounding residual exactly in RN. It is also exact when B = 0,
  // reducing the final addition to RN(exp_u.hi + exp_u.lo).
  FloatFloat sum_biased = fputil::exact_add(bias, exp_u.hi);
  float sum_rounded = sum_biased.hi + (sum_biased.lo + exp_u.lo);
  uint32_t res_u = FPBits(sum_rounded).uintval();

  // For subnormal results, sum_rounded lies in B's binade [B, 2B):
  // sum_rounded = B + M * ulp(B), where M is the integer encoded by
  // its fraction field. Clearing the exponent makes the same fraction
  // bits represent M * 2^-149 = (sum_rounded - B) * 2^k. This removes
  // the bias and applies the scale without further rounding.
  uint32_t res_mask = produces_subnormal ? FPBits::FRACTION_MASK : ~0U;
  uint32_t exponent_adj = produces_subnormal ? 0U
                                             : static_cast<uint32_t>(k)
                                                   << FPBits::FRACTION_LEN;
  return FPBits((res_u & res_mask) + exponent_adj).get_val();
}

} // namespace cr_expf_rn_internal

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

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_FLOAT_UTILS_H