[libc-commits] [libc] 46aa659 - [libc][math] Improve exp2f performance.

[Tue Ly via libc-commits]

Author: Tue Ly
Date: 2023-06-20T09:34:20-04:00
New Revision: 46aa659a32e78481d90950eda9d464fc22eaad31

URL: https://github.com/llvm/llvm-project/commit/46aa659a32e78481d90950eda9d464fc22eaad31
DIFF: https://github.com/llvm/llvm-project/commit/46aa659a32e78481d90950eda9d464fc22eaad31.diff

LOG: [libc][math] Improve exp2f performance.

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

Reviewed By: michaelrj

Differential Revision: https://reviews.llvm.org/D153134

Added: 
    

Modified: 
    libc/src/math/generic/CMakeLists.txt
    libc/src/math/generic/exp2f.cpp
    utils/bazel/llvm-project-overlay/libc/BUILD.bazel

Removed: 
    


################################################################################
diff  --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 098b69817fb20..7772c6d826905 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -561,6 +561,7 @@ add_entrypoint_object(
     ../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

diff  --git a/libc/src/math/generic/exp2f.cpp b/libc/src/math/generic/exp2f.cpp
index f2c005d5e00a0..b967a4bb45371 100644
--- a/libc/src/math/generic/exp2f.cpp
+++ b/libc/src/math/generic/exp2f.cpp
@@ -10,11 +10,13 @@
 #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>
 
@@ -33,13 +35,38 @@ LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
   uint32_t x_u = xbits.uintval();
   uint32_t x_abs = x_u & 0x7fff'ffffU;
 
-  // |x| < 2^-25
-  if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
-    return 1.0f + x;
-  }
+  // 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()) {
@@ -73,17 +100,6 @@ LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
     }
   }
 
-  // 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
@@ -101,11 +117,20 @@ LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
   // 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 =

diff  --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index eb7a3c87fc0f7..2d9d352ce5204 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -1200,6 +1200,7 @@ libc_math_function(
         ":__support_fputil_polyeval",
         ":__support_fputil_rounding_mode",
         ":__support_macros_optimization",
+        ":__support_macros_properties_cpu_features",
         ":explogxf",
     ],
 )