[libc][math] Improve exp2f performance.

[lntue]

Summary

Re-organize special cases and add a special case when |x| < 2^-5.

[lntue]

Improve performance for no-FMA, no-rounding instructions targets.

[michaelrj]

LGTM

libc/src/math/generic/CMakeLists.txt

 add_entrypoint_object(
   exp2f
   SRCS
     exp2f.cpp
   HDRS
     ../exp2f.h
   DEPENDS
     .explogxf
+    libc.src.__support.FPUtil.except_value_utils
     libc.src.__support.FPUtil.fenv_impl
     libc.src.__support.FPUtil.fp_bits
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.macros.optimization
     libc.include.errno

libc/src/math/generic/exp2f.cpp

 //===-- Single-precision 2^x function -------------------------------------===//
 //
 // 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
 //
 //===----------------------------------------------------------------------===//
 
 #include "src/math/exp2f.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/except_value_utils.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/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/properties/cpu_features.h"
 
 #include <errno.h>
 
 #include "explogxf.h"
 
 namespace __llvm_libc {
 
 constexpr uint32_t EXVAL1 = 0x3b42'9d37U;
 constexpr uint32_t EXVAL2 = 0xbcf3'a937U;
 constexpr uint32_t EXVAL_MASK = EXVAL1 & EXVAL2;
 
 LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
   uint32_t x_u = xbits.uintval();
   uint32_t x_abs = x_u & 0x7fff'ffffU;
 
+  // When |x| >= 128, or x is nan, or |x| <= 2^-5
+  if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U || x_abs <= 0x3d00'0000U)) {
+    // |x| <= 2^-5
+    if (x_abs <= 0x3d00'0000) {
       // |x| < 2^-25
       if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
         return 1.0f + x;
       }
 
-  // // When |x| >= 128, or x is nan
-  if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U)) {
+      // Check exceptional values.
+      if (LIBC_UNLIKELY((x_u & EXVAL_MASK) == EXVAL_MASK)) {
+        if (LIBC_UNLIKELY(x_u == EXVAL1)) { // x = 0x1.853a6ep-9f
+          return fputil::round_result_slightly_down(0x1.00870ap+0f);
+        } else if (LIBC_UNLIKELY(x_u == EXVAL2)) { // x = -0x1.e7526ep-6f
+          return fputil::round_result_slightly_down(0x1.f58d62p-1f);
+        }
+      }
+
+      // Minimax polynomial generated by Sollya with:
+      // > P = fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-5, 2^-5]);
+      constexpr double COEFFS[] = {
+          0x1.62e42fefa39f3p-1, 0x1.ebfbdff82c57bp-3,  0x1.c6b08d6f2d7aap-5,
+          0x1.3b2ab6fc92f5dp-7, 0x1.5d897cfe27125p-10, 0x1.43090e61e6af1p-13};
+      double xd = static_cast<double>(x);
+      double xsq = xd * xd;
+      double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]);
+      double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]);
+      double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]);
+      double p = fputil::polyeval(xsq, c0, c1, c2);
+      double r = fputil::multiply_add(p, xd, 1.0);
+      return r;
+    }
 
     // x >= 128
     if (!xbits.get_sign()) {
       // x is finite
       if (x_u < 0x7f80'0000U) {
         int rounding = fputil::quick_get_round();
         if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
           return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
if (x_u >= 0xc316'0000U) {
       if (x != 0.0f) {
         fputil::set_errno_if_required(ERANGE);
         fputil::raise_except_if_required(FE_UNDERFLOW);
       }
       return 0.0f;
     }
   }
 
-  // Check exceptional values.
-  if (LIBC_UNLIKELY((x_u & EXVAL_MASK) == EXVAL_MASK)) {
-    if (LIBC_UNLIKELY(x_u == EXVAL1)) { // x = 0x1.853a6ep-9f
-      if (fputil::fenv_is_round_to_nearest())
-        return 0x1.00870ap+0f;
-    } else if (LIBC_UNLIKELY(x_u == EXVAL2)) { // x = -0x1.e7526ep-6f
-      if (fputil::fenv_is_round_to_nearest())
-        return 0x1.f58d62p-1f;
-    }
-  }
-
   // For -150 < x < 128, to compute 2^x, we perform the following range
   // reduction: find hi, mid, lo such that:
   //   x = hi + mid + lo, in which
   //     hi is an integer,
   //     0 <= mid * 2^5 < 32 is an integer
   //     -2^(-6) <= lo <= 2^-6.
   // In particular,
   //   hi + mid = round(x * 2^5) * 2^(-5).
   // Then,
   //   2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
   // 2^mid is stored in the lookup table of 32 elements.
   // 2^lo is computed using a degree-5 minimax polynomial
   // generated by Sollya.
   // We perform 2^hi * 2^mid by simply add hi to the exponent field
   // of 2^mid.
 
   // kf = (hi + mid) * 2^5 = round(x * 2^5)
-  float kf = fputil::nearest_integer(x * 32.0f);
+  float kf;
+  int k;
+#ifdef LIBC_TARGET_CPU_HAS_NEAREST_INT
+  kf = fputil::nearest_integer(x * 32.0f);
+  k = static_cast<int>(kf);
+#else
+  constexpr float HALF[2] = {0.5f, -0.5f};
+  k = static_cast<int>(fputil::multiply_add(x, 32.0f, HALF[x < 0.0f]));
+  kf = static_cast<float>(k);
+#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT
+
   // dx = lo = x - (hi + mid) = x - kf * 2^(-5)
   double dx = fputil::multiply_add(-0x1.0p-5f, kf, x);
 
-  int k = static_cast<int>(kf);
   // hi = floor(kf * 2^(-4))
   // exp_hi = shift hi to the exponent field of double precision.
   int64_t exp_hi =
       static_cast<int64_t>(static_cast<uint64_t>(k >> ExpBase::MID_BITS)
                            << fputil::FloatProperties<double>::MANTISSA_WIDTH);
   // mh = 2^hi * 2^mid
   // mh_bits = bit field of mh
   int64_t mh_bits = ExpBase::EXP_2_MID[k & ExpBase::MID_MASK] + exp_hi;

utils/bazel/llvm-project-overlay/libc/BUILD.bazel

libc_math_function(
     name = "exp2f",
     additional_deps = [
         ":__support_fputil_fma",
         ":__support_fputil_multiply_add",
         ":__support_fputil_nearest_integer",
         ":__support_fputil_polyeval",
         ":__support_fputil_rounding_mode",
         ":__support_macros_optimization",
+        ":__support_macros_properties_cpu_features",
         ":explogxf",
     ],
 )
 
 libc_math_function(
     name = "logf",
     additional_deps = [
         ":__support_fputil_fma",