[libc][Experimental] Use nearest integer instructions to speed up exp*f.

[lntue]

Summary

Performance before:

$ ./perf.sh expf
LIBC-location: /home/lnt/experiment/llvm/llvm-project/build/projects/libc/lib/libllvmlibc.a
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH reciprocal throughput   : 9.469
System LIBC reciprocal throughput : 7.341
LIBC reciprocal throughput        : 11.764

$ ./perf.sh exp2f
LIBC-location: /home/lnt/experiment/llvm/llvm-project/build/projects/libc/lib/libllvmlibc.a
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH reciprocal throughput   : 9.659
System LIBC reciprocal throughput : 6.424
LIBC reciprocal throughput        : 20.633

$ ./perf.sh expm1f
LIBC-location: /home/lnt/experiment/llvm/llvm-project/build/projects/libc/lib/libllvmlibc.a
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH reciprocal throughput   : 9.942
System LIBC reciprocal throughput : 43.645
LIBC reciprocal throughput        : 11.431

Performance after:

$ ./perf.sh expf
LIBC-location: /home/lnt/experiment/llvm-project/build/projects/libc/lib/libllvmlibc.a
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH reciprocal throughput   : 9.514
System LIBC reciprocal throughput : 7.356
LIBC reciprocal throughput        : 8.840

$ ./perf.sh exp2f
LIBC-location: /home/lnt/experiment/llvm-project/build/projects/libc/lib/libllvmlibc.a
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH reciprocal throughput   : 9.672
System LIBC reciprocal throughput : 6.387
LIBC reciprocal throughput        : 8.674

$ ./perf.sh expm1f
LIBC-location: /home/lnt/experiment/llvm-project/build/projects/libc/lib/libllvmlibc.a
GNU libc version: 2.31
GNU libc release: stable
CORE-MATH reciprocal throughput   : 9.945
System LIBC reciprocal throughput : 43.685
LIBC reciprocal throughput        : 9.138

libc/src/__support/FPUtil/nearest_integer.h

if (unlikely(diff > 0.5))
       return r + 1.0;
     if (unlikely(diff < -0.5))
       return r - 1.0;
     return r;
   }
   return x;
 }
 
+static inline float nearest_integer(float x) {
+  if (x < 0x1p24f && x > -0x1p24f) {
+    float r = x < 0 ? (x - 0x1.0p23f) + 0x1.0p23f : (x + 0x1.0p23f) - 0x1.0p23f;
+    float diff = x - r;
+    // The expression above is correct for the default rounding mode, round-to-
+    // nearest, tie-to-even.  For other rounding modes, it might be off by 1,
+    // which is corrected below.
+    if (unlikely(diff > 0.5f))
+      return r + 1.0f;
+    if (unlikely(diff < -0.5))
+      return r - 1.0f;
+    return r;
+  }
+  return x;
+}
+
 } // namespace fputil
 } // namespace __llvm_libc
 
 #endif
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_H

libc/src/__support/FPUtil/x86_64/nearest_integer.h

 
 static inline double nearest_integer(double x) {
   __m128d xmm = _mm_set_sd(x); // NOLINT
   __m128d ymm =
       _mm_round_sd(xmm, xmm, _MM_ROUND_NEAREST | _MM_FROUND_NO_EXC); // NOLINT
   return ymm[0];
 }
 
+static inline float nearest_integer(float x) {
+  __m128 xmm = _mm_set_ss(x); // NOLINT
+  __m128 ymm =
+      _mm_round_ss(xmm, xmm, _MM_ROUND_NEAREST | _MM_FROUND_NO_EXC); // NOLINT
+  return ymm[0];
+}
+
 } // namespace fputil
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_NEAREST_INTEGER_H

libc/src/math/generic/CMakeLists.txt

add_entrypoint_object(
   expf
   SRCS
     expf.cpp
   HDRS
     ../expf.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   exp2f
   SRCS
     exp2f.cpp
   HDRS
     ../exp2f.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   expm1f
   SRCS
     expm1f.cpp
   HDRS
     ../expm1f.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   copysign
add_entrypoint_object(
   HDRS
     ../fmodf.h
   DEPENDS
     libc.include.math
     libc.src.__support.FPUtil.generic.fmod
   COMPILE_OPTIONS
     -O3
 )
 No newline at end of file

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/BasicOperations.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FMA.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 namespace __llvm_libc {
 
 // Lookup table for 2^(m * 2^(-6)) with m = 0, ..., 63.
 // Table is generated with Sollya as follow:

 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;
 
   // Exceptional values.
-  switch (x_u) {
-  case 0x3b42'9d37U: // x = 0x1.853a6ep-9f
-    if (fputil::get_round() == FE_TONEAREST)
-      return 0x1.00870ap+0f;
-    break;
-  case 0x3c02'a9adU: // x = 0x1.05535ap-7f
-    if (fputil::get_round() == FE_TONEAREST)
-      return 0x1.016b46p+0f;
-    break;
-  case 0x3ca6'6e26U: { // x = 0x1.4cdc4cp-6f
-    int round_mode = fputil::get_round();
-    if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
-      return 0x1.03a16ap+0f;
-    return 0x1.03a168p+0f;
-  }
-  case 0x3d92'a282U: // x = 0x1.254504p-4f
-    if (fputil::get_round() == FE_UPWARD)
-      return 0x1.0d0688p+0f;
-    return 0x1.0d0686p+0f;
-  case 0xbcf3'a937U: // x = -0x1.e7526ep-6f
-    if (fputil::get_round() == FE_TONEAREST)
-      return 0x1.f58d62p-1f;
-    break;
-  case 0xb8d3'd026U: // x = -0x1.a7a04cp-14f
-    if (fputil::get_round() == FE_TONEAREST)
-      return 0x1.fff6d2p-1f;
-    break;
-  }
+  // Wait for https://reviews.llvm.org/D129918
+  // switch (x_u) {
+  // case 0x3b42'9d37U: // x = 0x1.853a6ep-9f
+  //   if (fputil::get_round() == FE_TONEAREST)
+  //     return 0x1.00870ap+0f;
+  //   break;
+  // case 0x3c02'a9adU: // x = 0x1.05535ap-7f
+  //   if (fputil::get_round() == FE_TONEAREST)
+  //     return 0x1.016b46p+0f;
+  //   break;
+  // case 0x3ca6'6e26U: { // x = 0x1.4cdc4cp-6f
+  //   int round_mode = fputil::get_round();
+  //   if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
+  //     return 0x1.03a16ap+0f;
+  //   return 0x1.03a168p+0f;
+  // }
+  // case 0x3d92'a282U: // x = 0x1.254504p-4f
+  //   if (fputil::get_round() == FE_UPWARD)
+  //     return 0x1.0d0688p+0f;
+  //   return 0x1.0d0686p+0f;
+  // case 0xbcf3'a937U: // x = -0x1.e7526ep-6f
+  //   if (fputil::get_round() == FE_TONEAREST)
+  //     return 0x1.f58d62p-1f;
+  //   break;
+  // case 0xb8d3'd026U: // x = -0x1.a7a04cp-14f
+  //   if (fputil::get_round() == FE_TONEAREST)
+  //     return 0x1.fff6d2p-1f;
+  //   break;
+  // }
 
   // // When |x| >= 128, |x| < 2^-25, or x is nan
   if (unlikely(x_abs >= 0x4300'0000U || x_abs <= 0x3280'0000U)) {
     // |x| < 2^-25
     if (x_abs <= 0x3280'0000U) {
       return 1.0f + x;
     }
     // x >= 128
LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
   // Multiply by 2^hi is simply adding hi to the exponent field.  We store
   // exp(mid) in the lookup tables EXP_M.  exp(lo) is computed using a degree-4
   // minimax polynomial generated by Sollya.
 
   // x_hi = round(hi + mid).
   // The default rounding mode for float-to-int conversion in C++ is
   // round-toward-zero. To make it round-to-nearest, we add (-1)^sign(x) * 0.5
   // before conversion.
-  int x_hi =
-      static_cast<int>(x * 0x1.0p+6f + (xbits.get_sign() ? -0.5f : 0.5f));
+  float kf = fputil::nearest_integer(x * 0x1.0p+6f);
+  int x_hi = static_cast<int>(kf);
   // For 2-complement integers, arithmetic right shift is the same as dividing
   // by a power of 2 and then round down (toward negative infinity).
   int e_hi = (x_hi >> 6) +
              static_cast<int>(fputil::FloatProperties<double>::EXPONENT_BIAS);
   fputil::FPBits<double> exp_hi(
       static_cast<uint64_t>(e_hi)
       << fputil::FloatProperties<double>::MANTISSA_WIDTH);
   // mid = x_hi & 0x0000'003fU;
   double exp_hi_mid = static_cast<double>(exp_hi) * EXP_M[x_hi & 0x3f];
   // Subtract (hi + mid) from x to get lo.
-  x -= static_cast<float>(x_hi) * 0x1.0p-6f;
-  double xd = static_cast<double>(x);
+  double xd = fputil::multiply_add(kf, -0x1.0p-6f, x);
   // Degree-4 minimax polynomial generated by Sollya with the following
   // commands:
   //   > display = hexadecimal;
   //   > Q = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-7, 2^-7]);
   //   > Q;
   double exp_lo =
       fputil::polyeval(xd, 0x1p0, 0x1.62e42fefa2417p-1, 0x1.ebfbdff82f809p-3,
                        0x1.c6b0b92131c47p-5, 0x1.3b2ab6fb568a3p-7);
   double result = exp_hi_mid * exp_lo;
   return static_cast<float>(result);
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/expf.cpp

 //===-- Single-precision e^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/expf.h"
 #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
 #include "src/__support/FPUtil/BasicOperations.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FMA.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, expf, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
LLVM_LIBC_FUNCTION(float, expf, (float x)) {
   // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2
   // respectively.  exp(lo) is computed using a degree-4 minimax polynomial
   // generated by Sollya.
 
   // x_hi = (hi + mid) * 2^7 = round(x * 2^7).
   // The default rounding mode for float-to-int conversion in C++ is
   // round-toward-zero. To make it round-to-nearest, we add (-1)^sign(x) * 0.5
   // before conversion.
-  int x_hi = static_cast<int>(x * 0x1.0p7f + (xbits.get_sign() ? -0.5f : 0.5f));
+  float kf = fputil::nearest_integer(x * 0x1.0p7f);
+  double xd = static_cast<double>(fputil::multiply_add(kf, -0x1.0p-7f, x));
+  int x_hi = static_cast<int>(kf);
   // Subtract (hi + mid) from x to get lo.
-  x -= static_cast<float>(x_hi) * 0x1.0p-7f;
-  double xd = static_cast<double>(x);
   x_hi += 104 << 7;
   // hi = x_hi >> 7
   double exp_hi = EXP_M1[x_hi >> 7];
   // mid * 2^7 = x_hi & 0x0000'007fU;
   double exp_mid = EXP_M2[x_hi & 0x7f];
   // Degree-4 minimax polynomial generated by Sollya with the following
   // commands:
   //   > display = hexadecimal;

libc/src/math/generic/expm1f.cpp

 //===-- Single-precision e^x - 1 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/expm1f.h"
 #include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
 #include "src/__support/FPUtil/BasicOperations.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FMA.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
#endif // LIBC_TARGET_HAS_FMA
   //   hi + mid = round(x * 2^7) * 2^(-7).
   // Then,
   //   expm1(x) = exp(hi + mid + lo) - 1 = exp(hi) * exp(mid) * exp(lo) - 1.
   // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2
   // respectively.  exp(lo) is computed using a degree-4 minimax polynomial
   // generated by Sollya.
 
   // x_hi = hi + mid.
-  int x_hi = static_cast<int>(x * 0x1.0p7f + (xbits.get_sign() ? -0.5f : 0.5f));
+  float kf = fputil::nearest_integer(x * 0x1.0p7f);
+  int x_hi = static_cast<int>(kf);
   // Subtract (hi + mid) from x to get lo.
-  x -= static_cast<float>(x_hi) * 0x1.0p-7f;
-  double xd = static_cast<double>(x);
+  double xd = static_cast<double>(fputil::multiply_add(kf, -0x1.0p-7f, x));
   x_hi += 104 << 7;
   // hi = x_hi >> 7
   double exp_hi = EXP_M1[x_hi >> 7];
   // lo = x_hi & 0x0000'007fU;
   double exp_mid = EXP_M2[x_hi & 0x7f];
   double exp_hi_mid = exp_hi * exp_mid;
   // Degree-4 minimax polynomial generated by Sollya with the following
   // commands: