[libc][math] Improved performance of exp2f function.

[orex]

Summary

New exp2 function algorithm:

1. Improved performance: 19.211 vs 29.002 by core-math perf tool.
2. Improved accuracy. Only two special values left.
3. Lookup table size reduced twice.

Latest version benchmark:

GNU libc version: 2.31
GNU libc release: stable
8.009
5.610
8.176

[zimmermann6]

I confirm the new function is correctly rounded (for all rounding modes). For what concerns efficiency, here is what I get on a AMD EPYC 7282 with gcc 10.2.1 and clang 11.0.1-2:

zimmerma@biscotte:~/svn/core-math$ LIBM=/localdisk/zimmerma/llvm-project/build/projects/libc/lib/libllvmlibc.a CORE_MATH_PERF_MODE=rdtsc ./perf.sh exp2f
GNU libc version: 2.31
GNU libc release: stable
9.720
6.228
18.081


zimmerma@biscotte:~/svn/core-math$ LIBM=/localdisk/zimmerma/llvm-project/build/projects/libc/lib/libllvmlibc-bec8dff.a CORE_MATH_PERF_MODE=rdtsc ./perf.sh exp2f
GNU libc version: 2.31
GNU libc release: stable
9.714
6.215
23.455

The first run is with this patch, the second one with revision bec8dff of llvm-libc: the reciprocal throughput decreased from 23 to 18 cycles.

[orex]

Improved performance by fputil::nearest_integer function.

[lntue]

Overall the improvement looks good for me. Let's wait for @zimmermann6 to double check the numbers.

libc/src/math/generic/CMakeLists.txt

504

missing nearest_integer dependency.

libc/src/math/generic/exp2f.cpp

14

You should also add multiply_add header and dependency directly.

[orex]

Added explicit dependencies.

[orex]

OK. Thanks! let's wait.

[zimmermann6]

I confirm the reciprocal throughput decreased from 18 to 10 cycles (on the same machine as above):

zimmerma@biscotte:~/svn/core-math$ LIBM=/localdisk/zimmerma/llvm-project/build/projects/libc/lib/libllvmlibc.a CORE_MATH_PERF_MODE=rdtsc ./perf.sh exp2f
GNU libc version: 2.33
GNU libc release: release
9.728
7.085
10.040
zimmerma@biscotte:~/svn/core-math$ LIBM=/localdisk/zimmerma/llvm-project/build/projects/libc/lib/libllvmlibc.a CORE_MATH_PERF_MODE=rdtsc PERF_ARGS=--latency ./perf.sh exp2f
GNU libc version: 2.33
GNU libc release: release
37.106
29.520
48.515

Good work Kirill!

[orex]

Merged to last main, some cosmetic changes.

libc/src/math/generic/CMakeLists.txt

 
 add_entrypoint_object(
   exp2f
   SRCS
     exp2f.cpp
   HDRS
     ../exp2f.h
   DEPENDS
+    .common_constants
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
 )
 
 add_entrypoint_object(
   expm1f
   SRCS
     expm1f.cpp

libc/src/math/generic/common_constants.h

 extern const double EXP_M1[195];
 
 // Lookup table for exp(m * 2^(-7)) with m = 0, ..., 127.
 // Table is generated with Sollya as follow:
 // > display = hexadecimal;
 // > for i from 0 to 127 do { D(exp(i / 128)); };
 extern const double EXP_M2[128];
 
+static constexpr int EXP_bits_p = 5;
+static constexpr int EXP_num_p = 1 << EXP_bits_p;
+
+// Value to add to make integer digits always positive.
+// Needed for correct rounding to nearest value in both positive and negative
+// range.
+static constexpr unsigned int ADD_TO_POS = (1 << 14);
+
+// Max possible exp abs value
+// 150 is abs of maximum value, which can be passed to "2^" function
+static_assert(150 * EXP_num_p < ADD_TO_POS, "Incorrect ADD_TO_POS value.");
+
+// > for i from 0 to 31 do { D(2^(i / 32)) - 1; };
+extern const double EXP_2_POW[EXP_num_p];
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H

libc/src/math/generic/common_constants.cpp

const double EXP_M2[128] = {
     0x1.299be1f3e7f1cp1, 0x1.2bf16e7d2a38cp1, 0x1.2e4baacdb6614p1,
     0x1.30aaa04e80d05p1, 0x1.330e587b62b28p1, 0x1.3576dce33feadp1,
     0x1.37e437282d4eep1, 0x1.3a5670ff972edp1, 0x1.3ccd9432682b4p1,
     0x1.3f49aa9d30590p1, 0x1.41cabe304cb34p1, 0x1.4450d8f00edd4p1,
     0x1.46dc04f4e5338p1, 0x1.496c4c6b832dap1, 0x1.4c01b9950a111p1,
     0x1.4e9c56c731f5dp1, 0x1.513c2e6c731d7p1, 0x1.53e14b042f9cap1,
     0x1.568bb722dd593p1, 0x1.593b7d72305bbp1,
 };
+// Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
+// printf("%.13a,\n", d[i]);
+const double EXP_2_POW[EXP_num_p] = {
+    -0x1.2bec333018867p-2, -0x1.1c1142e274118p-2, -0x1.0bdd71829fcf2p-2,
+    -0x1.f69d99accc7b6p-3, -0x1.d4c6af7557c93p-3, -0x1.b23213cc8e86cp-3,
+    -0x1.8edb9f5703dc0p-3, -0x1.6abf137076a8ep-3, -0x1.45d819a94b14bp-3,
+    -0x1.20224341286e4p-3, -0x1.f332113d56b1fp-4, -0x1.a46f918837cb7p-4,
+    -0x1.53f391822dbc7p-4, -0x1.01b466423250ap-4, -0x1.5b505d5b6f268p-5,
+    -0x1.5f134923757f3p-6, 0x0.0000000000000p+0,  0x1.66c34c5615d0fp-6,
+    0x1.6ab0d9f3121ecp-5,  0x1.1301d0125b50ap-4,  0x1.72b83c7d517aep-4,
+    0x1.d4873168b9aa8p-4,  0x1.1c3d373ab11c3p-3,  0x1.4f4efa8fef709p-3,
+    0x1.837f0518db8a9p-3,  0x1.b8d39b9d54e55p-3,  0x1.ef5326091a112p-3,
+    0x1.13821818624b4p-2,  0x1.2ff6b54d8a89cp-2,  0x1.4d0ad5a753e07p-2,
+    0x1.6ac1f752150a5p-2,  0x1.891fac0e95613p-2};
 
 } // namespace __llvm_libc

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 "common_constants.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/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:
-// > display = hexadecimal;
-// > 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,
-};
+constexpr float mlp = EXP_num_p;
+constexpr float mmld = -1.0 / mlp;
+
+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;
 
-  // 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());
+      return x + FPBits::inf().get_val();
     }
     // 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));
+        return FPBits(FPBits::MIN_SUBNORMAL).get_val();
       if (x != 0.0f)
         errno = ERANGE;
       return 0.0f;
     }
   }
-  // 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^6 is an integer
-  //     -2^(-7) <= lo < 2^-7.
-  // In particular,
-  //   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-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));
-  // 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);
-  // 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);
+
+  if (unlikely(x_u & exval_mask) == exval_mask) {
+    if (unlikely(x_u == exval1)) { // x = 0x1.853a6ep-9f
+      if (fputil::get_round() == FE_TONEAREST)
+        return 0x1.00870ap+0f;
+    } else if (unlikely(x_u == exval2)) { // x = -0x1.e7526ep-6f
+      if (fputil::get_round() == FE_TONEAREST)
+        return 0x1.f58d62p-1f;
+    }
+  }
+
+  float kf = fputil::nearest_integer(x * mlp);
+  double dx = fputil::multiply_add(mmld, kf, x);
+  double mult_f, ml;
+  {
+    uint32_t ps = static_cast<int>(kf) + (1 << (EXP_bits_p - 1)) +
+                  (fputil::FPBits<double>::EXPONENT_BIAS << EXP_bits_p);
+    fputil::FPBits<double> bs;
+    bs.set_unbiased_exponent(ps >> EXP_bits_p);
+    ml = 1.0 + EXP_2_POW[ps & (EXP_num_p - 1)];
+    mult_f = bs.get_val();
+  }
+
+  // N[Table[Ln[2]^n/n!,{n,1,6}],30]
+  double pe = fputil::polyeval(
+      dx, 1.0, 0x1.62e42fefa39efp-1, 0x1.ebfbdff82c58fp-3, 0x1.c6b08d704a0c0p-5,
+      0x1.3b2ab6fba4e77p-7, 0x1.5d87fe78a6731p-10, 0x1.430912f86c787p-13);
+
+  return mult_f * ml * pe;
 }
 
 } // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

SRCS
     tan_test.cpp
   DEPENDS
     libc.src.math.tan
     libc.src.__support.FPUtil.fputil
 )
 
 add_fp_unittest(
   expm1f_test
     NEED_MPFR
   SUITE
     libc_math_unittests
   SRCS
     expm1f_test.cpp
   DEPENDS
     libc.include.errno
     libc.include.math
     libc.src.math.expm1f