[libc] Improve the performance of exp2f.

[lntue]

Summary

Reduce the range-reduction table size from 128 entries down to 64 entries, and
reduce the polynomial's degree from 6 down to 4.

Currently we use a degree-6 minimax polynomial on an interval of length 2^-7
around 0 to compute exp2f. Based on the suggestion of @santoshn and the RLIBM
project (https://github.com/rutgers-apl/rlibm-prog/blob/main/libm/float/exp2.c)
it is possible to have a good polynomial of degree-4 on a subinterval of length
2^(-6) to approximate 2^x.

We did try to either reduce the degree of the polynomial down to 3 or increase
the interval size to 2^(-5), but in both cases the number of exceptional values
exploded. So we settle with using a degree-4 polynomial of the interval of
size 2^(-6) around 0.

[michaelrj]

LGTM from a code perspective

[sivachandra]

Fine with structuring but please curtail the CL description to 80 character width wherever possible.

[lntue]

In D122346#3403670, @sivachandra wrote:

Fine with structuring but please curtail the CL description to 80 character width wherever possible.

The commit's message is the same as the title of the review patch, which is under 80 char. I just added more info about the patch to the summary part which should not be carried to the commit message.

[sivachandra]

Depending on how you submit, you might end up with this message in the git log. I do not know exactly what arc land will do or if you plan to use it to submit, but it could pick up the message from here when submitting.

[lntue]

In D122346#3403767, @sivachandra wrote:

Depending on how you submit, you might end up with this message in the git log. I do not know exactly what arc land will do or if you plan to use it to submit, but it could pick up the message from here when submitting.

You are right, arc land does copy the summary from phabricators to commit messages. I updated the summary to stick with 80-char-width.

[zimmermann6]

I confirm the results are still correctly rounded for all four rounding modes. For rounding to nearest the reciprocal throughput/latency decreased from 56/103.5 cycles to 26.6/76.0 cycles on a Core i5-4590.
As a comparison, the core-math code runs in 21.2/63.2 cycles.

[santoshn]

Looks good to me. Thanks for incorporating RLIBM style range reduction for exp2f. A degree-4 polynomial is impressive.

libc/src/math/generic/exp2f.cpp

 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 namespace __llvm_libc {
 
-// Lookup table for 2^(m * 2^(-7)) with m = 0, ..., 127.
+// Lookup table for 2^(m * 2^(-6)) with m = 0, ..., 63.
 // Table is generated with Sollya as follow:
 // > display = hexadecimal;
-// > for i from 0 to 127 do { D(2^(i / 128)); };
-static constexpr double EXP_M[128] = {
-    0x1.0000000000000p0, 0x1.0163da9fb3335p0, 0x1.02c9a3e778061p0,
-    0x1.04315e86e7f85p0, 0x1.059b0d3158574p0, 0x1.0706b29ddf6dep0,
-    0x1.0874518759bc8p0, 0x1.09e3ecac6f383p0, 0x1.0b5586cf9890fp0,
-    0x1.0cc922b7247f7p0, 0x1.0e3ec32d3d1a2p0, 0x1.0fb66affed31bp0,
-    0x1.11301d0125b51p0, 0x1.12abdc06c31ccp0, 0x1.1429aaea92de0p0,
-    0x1.15a98c8a58e51p0, 0x1.172b83c7d517bp0, 0x1.18af9388c8deap0,
-    0x1.1a35beb6fcb75p0, 0x1.1bbe084045cd4p0, 0x1.1d4873168b9aap0,
-    0x1.1ed5022fcd91dp0, 0x1.2063b88628cd6p0, 0x1.21f49917ddc96p0,
-    0x1.2387a6e756238p0, 0x1.251ce4fb2a63fp0, 0x1.26b4565e27cddp0,
-    0x1.284dfe1f56381p0, 0x1.29e9df51fdee1p0, 0x1.2b87fd0dad990p0,
-    0x1.2d285a6e4030bp0, 0x1.2ecafa93e2f56p0, 0x1.306fe0a31b715p0,
-    0x1.32170fc4cd831p0, 0x1.33c08b26416ffp0, 0x1.356c55f929ff1p0,
-    0x1.371a7373aa9cbp0, 0x1.38cae6d05d866p0, 0x1.3a7db34e59ff7p0,
-    0x1.3c32dc313a8e5p0, 0x1.3dea64c123422p0, 0x1.3fa4504ac801cp0,
-    0x1.4160a21f72e2ap0, 0x1.431f5d950a897p0, 0x1.44e086061892dp0,
-    0x1.46a41ed1d0057p0, 0x1.486a2b5c13cd0p0, 0x1.4a32af0d7d3dep0,
-    0x1.4bfdad5362a27p0, 0x1.4dcb299fddd0dp0, 0x1.4f9b2769d2ca7p0,
-    0x1.516daa2cf6642p0, 0x1.5342b569d4f82p0, 0x1.551a4ca5d920fp0,
-    0x1.56f4736b527dap0, 0x1.58d12d497c7fdp0, 0x1.5ab07dd485429p0,
-    0x1.5c9268a5946b7p0, 0x1.5e76f15ad2148p0, 0x1.605e1b976dc09p0,
-    0x1.6247eb03a5585p0, 0x1.6434634ccc320p0, 0x1.6623882552225p0,
-    0x1.68155d44ca973p0, 0x1.6a09e667f3bcdp0, 0x1.6c012750bdabfp0,
-    0x1.6dfb23c651a2fp0, 0x1.6ff7df9519484p0, 0x1.71f75e8ec5f74p0,
-    0x1.73f9a48a58174p0, 0x1.75feb564267c9p0, 0x1.780694fde5d3fp0,
-    0x1.7a11473eb0187p0, 0x1.7c1ed0130c132p0, 0x1.7e2f336cf4e62p0,
-    0x1.80427543e1a12p0, 0x1.82589994cce13p0, 0x1.8471a4623c7adp0,
-    0x1.868d99b4492edp0, 0x1.88ac7d98a6699p0, 0x1.8ace5422aa0dbp0,
-    0x1.8cf3216b5448cp0, 0x1.8f1ae99157736p0, 0x1.9145b0b91ffc6p0,
-    0x1.93737b0cdc5e5p0, 0x1.95a44cbc8520fp0, 0x1.97d829fde4e50p0,
-    0x1.9a0f170ca07bap0, 0x1.9c49182a3f090p0, 0x1.9e86319e32323p0,
-    0x1.a0c667b5de565p0, 0x1.a309bec4a2d33p0, 0x1.a5503b23e255dp0,
-    0x1.a799e1330b358p0, 0x1.a9e6b5579fdbfp0, 0x1.ac36bbfd3f37ap0,
-    0x1.ae89f995ad3adp0, 0x1.b0e07298db666p0, 0x1.b33a2b84f15fbp0,
-    0x1.b59728de5593ap0, 0x1.b7f76f2fb5e47p0, 0x1.ba5b030a1064ap0,
-    0x1.bcc1e904bc1d2p0, 0x1.bf2c25bd71e09p0, 0x1.c199bdd85529cp0,
-    0x1.c40ab5fffd07ap0, 0x1.c67f12e57d14bp0, 0x1.c8f6d9406e7b5p0,
-    0x1.cb720dcef9069p0, 0x1.cdf0b555dc3fap0, 0x1.d072d4a07897cp0,
-    0x1.d2f87080d89f2p0, 0x1.d5818dcfba487p0, 0x1.d80e316c98398p0,
-    0x1.da9e603db3285p0, 0x1.dd321f301b460p0, 0x1.dfc97337b9b5fp0,
-    0x1.e264614f5a129p0, 0x1.e502ee78b3ff6p0, 0x1.e7a51fbc74c83p0,
-    0x1.ea4afa2a490dap0, 0x1.ecf482d8e67f1p0, 0x1.efa1bee615a27p0,
-    0x1.f252b376bba97p0, 0x1.f50765b6e4540p0, 0x1.f7bfdad9cbe14p0,
-    0x1.fa7c1819e90d8p0, 0x1.fd3c22b8f71f1p0,
+// > for i from 0 to 63 do { D(2^(i / 64)); };
+static constexpr double EXP_M[64] = {
+    0x1.0000000000000p0, 0x1.02c9a3e778061p0, 0x1.059b0d3158574p0,
+    0x1.0874518759bc8p0, 0x1.0b5586cf9890fp0, 0x1.0e3ec32d3d1a2p0,
+    0x1.11301d0125b51p0, 0x1.1429aaea92de0p0, 0x1.172b83c7d517bp0,
+    0x1.1a35beb6fcb75p0, 0x1.1d4873168b9aap0, 0x1.2063b88628cd6p0,
+    0x1.2387a6e756238p0, 0x1.26b4565e27cddp0, 0x1.29e9df51fdee1p0,
+    0x1.2d285a6e4030bp0, 0x1.306fe0a31b715p0, 0x1.33c08b26416ffp0,
+    0x1.371a7373aa9cbp0, 0x1.3a7db34e59ff7p0, 0x1.3dea64c123422p0,
+    0x1.4160a21f72e2ap0, 0x1.44e086061892dp0, 0x1.486a2b5c13cd0p0,
+    0x1.4bfdad5362a27p0, 0x1.4f9b2769d2ca7p0, 0x1.5342b569d4f82p0,
+    0x1.56f4736b527dap0, 0x1.5ab07dd485429p0, 0x1.5e76f15ad2148p0,
+    0x1.6247eb03a5585p0, 0x1.6623882552225p0, 0x1.6a09e667f3bcdp0,
+    0x1.6dfb23c651a2fp0, 0x1.71f75e8ec5f74p0, 0x1.75feb564267c9p0,
+    0x1.7a11473eb0187p0, 0x1.7e2f336cf4e62p0, 0x1.82589994cce13p0,
+    0x1.868d99b4492edp0, 0x1.8ace5422aa0dbp0, 0x1.8f1ae99157736p0,
+    0x1.93737b0cdc5e5p0, 0x1.97d829fde4e50p0, 0x1.9c49182a3f090p0,
+    0x1.a0c667b5de565p0, 0x1.a5503b23e255dp0, 0x1.a9e6b5579fdbfp0,
+    0x1.ae89f995ad3adp0, 0x1.b33a2b84f15fbp0, 0x1.b7f76f2fb5e47p0,
+    0x1.bcc1e904bc1d2p0, 0x1.c199bdd85529cp0, 0x1.c67f12e57d14bp0,
+    0x1.cb720dcef9069p0, 0x1.d072d4a07897cp0, 0x1.d5818dcfba487p0,
+    0x1.da9e603db3285p0, 0x1.dfc97337b9b5fp0, 0x1.e502ee78b3ff6p0,
+    0x1.ea4afa2a490dap0, 0x1.efa1bee615a27p0, 0x1.f50765b6e4540p0,
+    0x1.fa7c1819e90d8p0,
 };
 
 INLINE_FMA
 LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
-  // When x =< -150 or nan
-  if (unlikely(xbits.uintval() >= 0xc316'0000U)) {
+  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;
+  }
+
+  // // 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
+    if (!xbits.get_sign()) {
+      // x is finite
+      if (x_u < 0x7f80'0000U) {
+        int rounding = fputil::get_round();
+        if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+          return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
+
+        errno = ERANGE;
+      }
+      // x is +inf or nan
+      return x + static_cast<float>(FPBits::inf());
+    }
+    // x < -150
+    if (x_u >= 0xc316'0000U) {
       // exp(-Inf) = 0
       if (xbits.is_inf())
         return 0.0f;
       // exp(nan) = nan
       if (xbits.is_nan())
         return x;
       if (fputil::get_round() == FE_UPWARD)
         return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL));
       if (x != 0.0f)
         errno = ERANGE;
       return 0.0f;
     }
-  // x >= 128 or nan
-  if (unlikely(!xbits.get_sign() && (xbits.uintval() >= 0x4300'0000U))) {
-    if (xbits.uintval() < 0x7f80'0000U) {
-      int rounding = fputil::get_round();
-      if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
-        return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
-
-      errno = ERANGE;
-    }
-    return x + static_cast<float>(FPBits::inf());
-  }
-  // |x| < 2^-25
-  if (unlikely(xbits.get_unbiased_exponent() <= 101)) {
-    return 1.0f + x;
   }
-  // Exceptional values.
-  switch (xbits.uintval()) {
-  case 0x3b42'9d37U: // x = 0x1.853a6ep-9f
-    if (fputil::get_round() == FE_TONEAREST)
-      return 0x1.00870ap+0f;
-    break;
-  case 0xbcf3'a937U: // x = -0x1.e7526ep-6f
-    if (fputil::get_round() == FE_TONEAREST)
-      return 0x1.f58d62p-1f;
-    break;
-  }
-
   // 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,
-  //     mid * 2^7 is an integer
-  //     -2^(-8) <= lo < 2^-8.
+  //     mid * 2^6 is an integer
+  //     -2^(-7) <= lo < 2^-7.
   // In particular,
-  //   hi + mid = round(x * 2^7) * 2^(-7).
+  //   hi + mid = round(x * 2^6) * 2^(-6).
   // Then,
   //   2^(x) = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
   // 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-7
+  // exp(mid) in the lookup tables EXP_M.  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);
-  // Subtract (hi + mid) from x to get lo.
-  x -= static_cast<float>(x_hi) * 0x1.0p-7f;
-  double xd = static_cast<double>(x);
-  // Make sure that -2^(-8) <= lo < 2^-8.
-  if (x >= 0x1.0p-8f) {
-    ++x_hi;
-    xd -= 0x1.0p-7;
-  }
-  if (x < -0x1.0p-8f) {
-    --x_hi;
-    xd += 0x1.0p-7;
-  }
+  // 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));
   // 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 hi = x_hi >> 7;
-  // mid = x_hi & 0x0000'007fU;
-  double exp_mid = EXP_M[x_hi & 0x7f];
-  // Degree-6 minimax polynomial generated by Sollya with the following
+  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);
+  // Degree-4 minimax polynomial generated by Sollya with the following
   // commands:
   //   > display = hexadecimal;
-  //   > Q = fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-8, 2^-8]);
+  //   > Q = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-7, 2^-7]);
   //   > Q;
   double exp_lo =
-      fputil::polyeval(xd, 0x1p0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58ep-3,
-                       0x1.c6b08d711fe2fp-5, 0x1.3b2ab6fe3deb5p-7,
-                       0x1.5d72a05f45c04p-10, 0x1.4284d40c33326p-13);
-  fputil::FPBits<double> result(exp_mid * exp_lo);
-  result.set_unbiased_exponent(static_cast<uint16_t>(
-      static_cast<int>(result.get_unbiased_exponent()) + hi));
-  return static_cast<float>(static_cast<double>(result));
+      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/test/src/math/exp2f_test.cpp

TEST(LlvmLibcExp2fTest, Overflow) {
 
   EXPECT_FP_EQ(inf, __llvm_libc::exp2f(float(FPBits(0x43000000U))));
   EXPECT_MATH_ERRNO(ERANGE);
 
   EXPECT_FP_EQ(inf, __llvm_libc::exp2f(float(FPBits(0x43000001U))));
   EXPECT_MATH_ERRNO(ERANGE);
 }
 
-// Test with inputs which are the borders of underflow/overflow but still
-// produce valid results without setting errno.
-TEST(LlvmLibcExp2fTest, Borderline) {
-  float x;
-
+TEST(LlvmLibcExp2fTest, TrickyInputs) {
+  constexpr int N = 12;
+  constexpr uint32_t INPUTS[N] = {
+      0x3b429d37U, /*0x1.853a6ep-9f*/
+      0x3c02a9adU, /*0x1.05535ap-7f*/
+      0x3ca66e26U, /*0x1.4cdc4cp-6f*/
+      0x3d92a282U, /*0x1.254504p-4f*/
+      0x42fa0001U, /*0x1.f40002p+6f*/
+      0x42ffffffU, /*0x1.fffffep+6f*/
+      0xb8d3d026U, /*-0x1.a7a04cp-14f*/
+      0xbcf3a937U, /*-0x1.e7526ep-6f*/
+      0xc2fa0001U, /*-0x1.f40002p+6f*/
+      0xc2fc0000U, /*-0x1.f8p+6f*/
+      0xc2fc0001U, /*-0x1.f80002p+6f*/
+      0xc3150000U, /*-0x1.2ap+7f*/
+  };
+  for (int i = 0; i < N; ++i) {
     errno = 0;
-  x = float(FPBits(0x42fa0001U));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0x42ffffffU));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0xc2fa0001U));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0xc2fc0000U));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0xc2fc0001U));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0xc3150000U));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0x3b42'9d37U));
-  EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
-                                 __llvm_libc::exp2f(x), 0.5);
-  EXPECT_MATH_ERRNO(0);
-
-  x = float(FPBits(0xbcf3'a937U));
+    float x = float(FPBits(INPUTS[i]));
     EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
                                    __llvm_libc::exp2f(x), 0.5);
     EXPECT_MATH_ERRNO(0);
   }
+}
 
 TEST(LlvmLibcExp2fTest, Underflow) {
   errno = 0;
   EXPECT_FP_EQ(0.0f, __llvm_libc::exp2f(float(FPBits(0xff7fffffU))));
   EXPECT_MATH_ERRNO(ERANGE);
 
   float x = float(FPBits(0xc3158000U));
   EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,