libc/src/__support/math/exp2f_integer_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 integer-only single-precision exponential
/// functions.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_INTEGER_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_INTEGER_UTILS_H

#include "hdr/fenv_macros.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/frac64.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/macros/properties/types.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace static_rounding {

// Degree-11 minimax polynomial approximating 2^u - 1 on [0, 1], generated by
// Sollya with:
// > P = fpminimax(2^x, 11, [|1, 64...|], [0, 1], absolute, fixed);
// > dirtyinfnorm(2^x - P(x), [0, 1]);
//   0x1.6238...p-58
// Store the fractional part of the coefficients below with LSB = 2^-64:
LIBC_INLINE_VAR constexpr Frac64 EXP2F_COEFFS[11] = {
    Frac64(0xb172'17f7'd1cf'b7cf), // u
    Frac64(0x3d7f'7bff'057d'4a5e), // u^2
    Frac64(0x0e35'846b'8363'9484), // u^3
    Frac64(0x0276'556d'ec97'dcd4), // u^4
    Frac64(0x0057'61ff'dc04'c7ff), // u^5
    Frac64(0x000a'1847'b6e7'92ec), // u^6
    Frac64(0x0000'ffe8'14e5'7033), // u^7
    Frac64(0x0000'1628'b6e9'70c8), // u^8
    Frac64(0x0000'01b8'8ce7'4088), // u^9
    Frac64(0x0000'001c'18d5'cb29), // u^10
    Frac64(0x0000'0002'b43f'4490), // u^11
};

// Compute |x| * log2(b) in unsigned 64-bit fixed-point format with binary
// point at bit 54 (LSB = 2^-54), for x_e in [max_e - 31, max_e].
// The binary point of log2_b is at bit (63 - (133 - max_e)).
LIBC_ALWAYS_INLINE uint64_t exp_b_scale(uint32_t x_u, uint32_t x_e,
                                        Frac64 log2_b, uint32_t max_e) {
#if defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
  // When 128-bit integers are available (64-bit targets) and small code size is
  // not requested, variable 64-bit shifts and 64x64 -> 128 multiplications are
  // single native instructions. Shift the 24-bit significand of |x| into a
  // 64-bit word and take the upper 64 bits of its 128-bit product with log2_b.
  uint64_t m = static_cast<uint64_t>((x_u & 0x007f'ffffU) | 0x0080'0000U)
               << (x_e - (max_e - 38U));
  return (Frac64(m) * log2_b).val[0];
#else  // !defined(LIBC_TYPES_HAS_INT128) || defined(LIBC_MATH_HAS_SMALL_TABLES)
  // On 32-bit targets without native 128-bit integers (or when optimizing for
  // size), avoid 64-bit variable shifts (which lower to branches on 32-bit
  // targets) and 64x64 -> 128 multiplications:
  // 1. Place the 24-bit significand of |x| at bits [8..31] (leading 1 at
  //    bit 31) and compute the upper 64 bits of the 32x64 -> 96 product with
  //    (log2_b >> 2) using only two 32x32 -> 64 multiplications.
  // 2. Since x_e is in [max_e - 31, max_e], s = max_e - x_e is in [0, 31], so
  //    shifting the 64-bit product right by s using two 32-bit shifts places
  //    the binary point of |x| * log2(b) at bit 54 without branches or UB at
  //    s = 0.
  uint64_t m = (x_u << 8) | 0x8000'0000U;
  uint64_t log2_b_scaled = log2_b.val[0] >> 2;
  uint64_t prod = m * (log2_b_scaled >> 32) +
                  ((m * static_cast<uint32_t>(log2_b_scaled)) >> 32);
  uint32_t s = max_e - x_e;
  uint32_t prod_lo = static_cast<uint32_t>(prod);
  uint32_t prod_hi = static_cast<uint32_t>(prod >> 32);
  uint32_t u_lo = (prod_lo >> s) | ((prod_hi << 1) << (31U - s));
  uint32_t u_hi = prod_hi >> s;
  return (static_cast<uint64_t>(u_hi) << 32) | u_lo;
#endif // defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
}

// Evaluate 2^(x * log2(b)) from u_bits = |x| * log2(b) (with binary point at
// bit 54) and return the IEEE 754 single-precision bit representation.
LIBC_ALWAYS_INLINE uint32_t exp2f_eval_bits(uint64_t u_bits, bool is_neg,
                                            [[maybe_unused]] int rounding) {
  // Negate in two's complement when x < 0 so that floor(x * log2(b)) is the
  // arithmetic right shift by 54, and the fractional part u in [0, 1) is the
  // low 54 bits.
  // Setting the LSB (| 1U) prevents InstCombine's demanded-bits analysis from
  // narrowing the free 32-bit low-word extraction mask 0xffff'ffffU to
  // 0xffff'fc00U inside quick_mul_hi on 32-bit targets. Even when the
  // fractional part is 0, u = Frac64(1) satisfies (a * 1) >> 64 == 0, so p
  // remains 0 in all rounding modes.
  int64_t u_signed =
      is_neg ? -static_cast<int64_t>(u_bits) : static_cast<int64_t>(u_bits);
  int k =
      (static_cast<int32_t>(static_cast<uint64_t>(u_signed) >> 32) >> 22) + 126;
  Frac64 u((static_cast<uint64_t>(u_signed) << 10) | 1U);

#if defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
  // When 128-bit integers are available (64-bit targets) and small code size is
  // not requested, each Frac64 multiplication is a single native 64x64 -> 128
  // multiply-high instruction, and out-of-order 64-bit CPUs have ample
  // registers and multiple integer multiplier pipelines:
  // 1. Evaluate p ~ 2^u - 1 on [0, 1) using an unrolled 2-way Estrin/Horner
  //    scheme (evaluating even and odd powers in parallel in u^2) to halve the
  //    critical-path latency of the polynomial evaluation.
  // 2. Take a fast return path for the common case of normal outputs in
  //    round-to-nearest mode (k >= 0 && rounding == FE_TONEAREST), bypassing
  //    the subnormal and directed-rounding logic below.
  Frac64 u2 = u * u;
  Frac64 p_even =
      fputil::polyeval(u2, EXP2F_COEFFS[0], EXP2F_COEFFS[2], EXP2F_COEFFS[4],
                       EXP2F_COEFFS[6], EXP2F_COEFFS[8], EXP2F_COEFFS[10]);
  Frac64 p_odd =
      fputil::polyeval(u2, EXP2F_COEFFS[1], EXP2F_COEFFS[3], EXP2F_COEFFS[5],
                       EXP2F_COEFFS[7], EXP2F_COEFFS[9]);
  Frac64 p = fputil::multiply_add(p_odd, u, p_even) * u;

  if (LIBC_LIKELY(k >= 0 && rounding == FE_TONEAREST))
    return (((static_cast<uint32_t>(p.val[0] >> 40) + 1U) >> 1) +
            (static_cast<uint32_t>(k + 1) << 23));
#else  // !defined(LIBC_TYPES_HAS_INT128) || defined(LIBC_MATH_HAS_SMALL_TABLES)
  // On 32-bit targets without native 128-bit integers (such as 32-bit ARM
  // microcontrollers without hardware double-precision FPU), or when
  // LIBC_MATH_HAS_SMALL_TABLES is enabled, unrolling or using Estrin's scheme
  // would increase register pressure and bloat code size (by ~470-810 bytes on
  // 32-bit targets and ~70-360 bytes on 64-bit targets). Instead, evaluate
  // p ~ 2^u - 1 using a compact, non-unrolled Horner loop over EXP2F_COEFFS and
  // share the unified normal/subnormal rounding epilogue below.
  Frac64 p(0);
  LIBC_LOOP_NOUNROLL
  for (int i = 10; i >= 0; --i)
    p = (p + EXP2F_COEFFS[i]) * u;
#endif // defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)

  // Combine 2^u = 1 + p with 2^(k - 126):
  // Place the explicit leading 1 of 2^u at bit 31 and the top 31 bits of p in
  // bits [0..30].
  // 1. Normal result (k >= 0, biased exponent k + 1 in [1, 254]):
  //    Shift right by 7 to place the leading 1 at bit 24 (which after >> 1
  //    contributes +1 to the biased exponent field, added to k << 23).
  // 2. Subnormal result (k < 0, i.e., k in [-24, -1]):
  //    Shift right by 7 - k in [8, 31] to place both the leading 1 and the
  //    fractional bits directly at their subnormal bit positions.
  uint32_t mant = static_cast<uint32_t>(p.val[0] >> 33) | 0x8000'0000U;
  uint32_t shift = 7;
  uint32_t exp_bits = static_cast<uint32_t>(k) << 23;

  if (LIBC_UNLIKELY(k < 0)) {
    shift = static_cast<uint32_t>(7 - k);
    exp_bits = 0;
  }

#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
  if (LIBC_UNLIKELY(rounding != FE_TONEAREST)) {
    uint32_t round_up = 0;
    if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
      uint32_t mask = ~0U >> (31U - shift);
      round_up = static_cast<uint32_t>(
          ((mant & mask) | static_cast<uint32_t>(p.val[0])) != 0);
    }

    return ((mant >> shift) >> 1) + round_up + exp_bits;
  }
#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY

  // Round to nearest.
  return (((mant >> shift) + 1U) >> 1) + exp_bits;
}

} // namespace static_rounding
} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_INTEGER_UTILS_H