Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[libc][math][c23] Add exp2m1f C23 math function #86996

Merged
merged 6 commits into from Apr 4, 2024

Conversation

overmighty
Copy link
Contributor

@overmighty overmighty commented Mar 28, 2024

Fixes #86502.

cc @lntue

}

// x <= -25
if (x <= -25.0f) {
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In exp2f you compare x_u. I assume it's a small optimization based on integer comparisons usually being cheaper than floating-point ones. For example, on x86-64, cmp seems cheaper than {,v}ucomiss:

I can change this if you want. Also, the comment is pretty useless if I don't change it.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can benchmark later, but there are several reasons:

  • the first one is that inside CPUs, floating point instructions and integer instructions use different ports, and switching values between them does have some hidden cost, like the above we use integer comparisons to omit inf / nan.
  • second is that with multiple tight conditions, keeping them as integers will help compilers to optimize better, since the full semantic of floating point is more complicated, such as it's easy for us to see that NaNs cannot get here, but not so easy for the compilers.

And of course everything is just guessing. Once you finish the implementations and accuracy tests, you can try both with performance tests to see if it has any affects.

libc/src/math/generic/exp2m1f.cpp Outdated Show resolved Hide resolved
libc/test/src/math/exhaustive/CMakeLists.txt Show resolved Hide resolved
libc/test/src/math/smoke/exp2m1f_test.cpp Outdated Show resolved Hide resolved
libc/src/math/generic/exp2m1f.cpp Show resolved Hide resolved
@lntue lntue self-requested a review March 28, 2024 20:39
@overmighty
Copy link
Contributor Author

overmighty commented Mar 29, 2024

mpfr_exp2m1 isn't available in the CI environment.

It seems to have been added in MPFR 4.2.0: https://gitlab.inria.fr/mpfr/mpfr/-/commit/78f1866613ca67b801ce4c17f103e2bc1e576fc1.

I don't know what OS the CI runs but MPFR >= 4.2.0 isn't available on any Ubuntu LTS release: https://packages.ubuntu.com/search?suite=default&section=all&arch=any&keywords=libmpfr-dev&searchon=names.

Should I add #if guards using MPFR_VERSION_MAJOR and MPFR_VERSION_MINOR for now?

@lntue
Copy link
Contributor

lntue commented Mar 29, 2024

mpfr_exp2m1 isn't available in the CI environment.

It seems to have been added in MPFR 4.2.0: https://gitlab.inria.fr/mpfr/mpfr/-/commit/78f1866613ca67b801ce4c17f103e2bc1e576fc1.

I don't know what OS the CI runs but MPFR >= 4.2.0 isn't available on any Ubuntu LTS release: https://packages.ubuntu.com/search?suite=default&section=all&arch=any&keywords=libmpfr-dev&searchon=names.

Should I add #if guards using MPFR_VERSION_MAJOR and MPFR_VERSION_MINOR for now?

You can implement mpfr_exp2m1 manually using mpfr_exp2 and maybe using precision * 3 for the intermediate computations. And then add a TODO comment to switch to use mpfr_exp2m1 once MPFR 4.2.0 is available in all the bots.

@llvmbot
Copy link
Collaborator

llvmbot commented Apr 1, 2024

@llvm/pr-subscribers-libc

Author: OverMighty (overmighty)

Changes

Fixes #86502.

TODO:

  • Add unit tests with MPFR. (Once the smoke tests passed I directly jumped onto the exhaustive tests.)

  • Handle exceptional values found by exhaustive tests.
    <details>
    <summary>Exhaustive tests output</summary>
    <br>

    ~/projects/llvm-project libc-math-exp2m1f* ➜ time ./build/projects/libc/test/src/math/exhaustive/libc.test.src.math.exhaustive.exp2m1f_test.__unit__
    [ RUN      ] LlvmLibcExp2m1fExhaustiveTest.PostiveRange
    -- Testing for FE_TONEAREST in range [0x0, 0x7f800000) --
    Test failed for 1 inputs in range: 764411904 to 765460480 [0x2d900000, 0x2da00000), [0x1.2p-36, 0x1.4p-36)
    Test failed for 1 inputs in range: 1056964608 to 1058013184 [0x3f000000, 0x3f100000), [0x1p-1, 0x1.2p-1)
    100% is in process
    Test FAILED
    /home/overmighty/projects/llvm-project/libc/test/src/math/exhaustive/exhaustive_test.h:148: FAILURE
          Expected: failed.load()
          Which is: 2
    To be equal to: uint64_t(0)
          Which is: 0
    -- Testing for FE_UPWARD in range [0x0, 0x7f800000) --
    Test failed for 1 inputs in range: 907018240 to 908066816 [0x36100000, 0x36200000), [0x1.2p-19, 0x1.4p-19)
    Test failed for 1 inputs in range: 903872512 to 904921088 [0x35e00000, 0x35f00000), [0x1.cp-20, 0x1.ep-20)
    Test failed for 1 inputs in range: 947912704 to 948961280 [0x38800000, 0x38900000), [0x1p-14, 0x1.2p-14)
    100% is in process
    Test FAILED
    /home/overmighty/projects/llvm-project/libc/test/src/math/exhaustive/exhaustive_test.h:148: FAILURE
          Expected: failed.load()
          Which is: 3
    To be equal to: uint64_t(0)
          Which is: 0
    -- Testing for FE_DOWNWARD in range [0x0, 0x7f800000) --
    Test failed for 1 inputs in range: 903872512 to 904921088 [0x35e00000, 0x35f00000), [0x1.cp-20, 0x1.ep-20)
    Test failed for 1 inputs in range: 947912704 to 948961280 [0x38800000, 0x38900000), [0x1p-14, 0x1.2p-14)
    100% is in process
    Test FAILED
    /home/overmighty/projects/llvm-project/libc/test/src/math/exhaustive/exhaustive_test.h:148: FAILURE
          Expected: failed.load()
          Which is: 2
    To be equal to: uint64_t(0)
          Which is: 0
    -- Testing for FE_TOWARDZERO in range [0x0, 0x7f800000) --
    Test failed for 1 inputs in range: 903872512 to 904921088 [0x35e00000, 0x35f00000), [0x1.cp-20, 0x1.ep-20)
    Test failed for 1 inputs in range: 947912704 to 948961280 [0x38800000, 0x38900000), [0x1p-14, 0x1.2p-14)
    100% is in process
    Test FAILED
    /home/overmighty/projects/llvm-project/libc/test/src/math/exhaustive/exhaustive_test.h:148: FAILURE
          Expected: failed.load()
          Which is: 2
    To be equal to: uint64_t(0)
          Which is: 0
    [  FAILED  ] LlvmLibcExp2m1fExhaustiveTest.PostiveRange
    [ RUN      ] LlvmLibcExp2m1fExhaustiveTest.NegativeRange
    -- Testing for FE_TONEAREST in range [0x80000000, 0xff800000) --
    Test failed for 1 inputs in range: 2898264064 to 2899312640 [0xacc00000, 0xacd00000), [-0x1.8p-38, -0x1.ap-38)
    100% is in process
    Test FAILED
    /home/overmighty/projects/llvm-project/libc/test/src/math/exhaustive/exhaustive_test.h:148: FAILURE
          Expected: failed.load()
          Which is: 1
    To be equal to: uint64_t(0)
          Which is: 0
    -- Testing for FE_UPWARD in range [0x80000000, 0xff800000) --
    Test failed for 1 inputs in range: 2946498560 to 2947547136 [0xafa00000, 0xafb00000), [-0x1.4p-32, -0x1.6p-32)
    Test failed for 1 inputs in range: 3134193664 to 3135242240 [0xbad00000, 0xbae00000), [-0x1.ap-10, -0x1.cp-10)
    Test failed for 1 inputs in range: 3169845248 to 3170893824 [0xbcf00000, 0xbd000000), [-0x1.ep-6, -0x1p-5)
    Test failed for 1 inputs in range: 3175088128 to 3176136704 [0xbd400000, 0xbd500000), [-0x1.8p-5, -0x1.ap-5)
    100% is in process
    Test FAILED
    /home/overmighty/projects/llvm-project/libc/test/src/math/exhaustive/exhaustive_test.h:148: FAILURE
          Expected: failed.load()
          Which is: 4
    To be equal to: uint64_t(0)
          Which is: 0
    -- Testing for FE_DOWNWARD in range [0x80000000, 0xff800000) --
    100% is in process
    Test PASSED
    -- Testing for FE_TOWARDZERO in range [0x80000000, 0xff800000) --
    100% is in process
    Test PASSED
    [  FAILED  ] LlvmLibcExp2m1fExhaustiveTest.NegativeRange
    Ran 2 tests.  PASS: 0  FAIL: 2
      21232.00s user 2.42s system 1992% cpu 17:45.45 total
    ~/projects/llvm-project libc-math-exp2m1f* ➜
    

    </details>

cc @lntue


Patch is 22.57 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/86996.diff

15 Files Affected:

  • (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
  • (modified) libc/docs/math/index.rst (+2)
  • (modified) libc/spec/stdc.td (+2)
  • (modified) libc/src/math/CMakeLists.txt (+2)
  • (added) libc/src/math/exp2m1f.h (+18)
  • (modified) libc/src/math/generic/CMakeLists.txt (+21)
  • (added) libc/src/math/generic/exp2m1f.cpp (+182)
  • (modified) libc/test/src/math/CMakeLists.txt (+16)
  • (modified) libc/test/src/math/exhaustive/CMakeLists.txt (+15)
  • (added) libc/test/src/math/exhaustive/exp2m1f_test.cpp (+33)
  • (added) libc/test/src/math/exp2m1f_test.cpp (+67)
  • (modified) libc/test/src/math/smoke/CMakeLists.txt (+11)
  • (added) libc/test/src/math/smoke/exp2m1f_test.cpp (+68)
  • (modified) libc/utils/MPFRWrapper/MPFRUtils.cpp (+39)
  • (modified) libc/utils/MPFRWrapper/MPFRUtils.h (+1)
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5b428e51aee620..b1e0592278e18c 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -369,6 +369,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.exp10f
     libc.src.math.exp2
     libc.src.math.exp2f
+    libc.src.math.exp2m1f
     libc.src.math.expm1
     libc.src.math.expm1f
     libc.src.math.fabs
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
index 080b6a4427f511..02d3dbfd7ba3bf 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -480,6 +480,8 @@ Higher Math Functions
 +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
 | exp2l      |         |         |         |         |         |         |         |         |         |         |         |         |
 +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
+| exp2m1f    | |check| |         |         |         |         |         |         |         |         |         |         |         |
++------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
 | expm1      | |check| | |check| | |check| | |check| | |check| |         |         | |check| | |check| | |check| |         |         |
 +------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
 | expm1f     | |check| | |check| | |check| | |check| | |check| |         |         | |check| | |check| | |check| |         |         |
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index ac6e1d1801ba55..bd72ff9095e2ea 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -535,6 +535,8 @@ def StdC : StandardSpec<"stdc"> {
           FunctionSpec<"exp2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"exp2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
+          FunctionSpec<"exp2m1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+
           FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index b67ee3a8e0dfbd..c89792b8ac7be6 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -88,6 +88,8 @@ add_math_entrypoint_object(expf)
 add_math_entrypoint_object(exp2)
 add_math_entrypoint_object(exp2f)
 
+add_math_entrypoint_object(exp2m1f)
+
 add_math_entrypoint_object(exp10)
 add_math_entrypoint_object(exp10f)
 
diff --git a/libc/src/math/exp2m1f.h b/libc/src/math/exp2m1f.h
new file mode 100644
index 00000000000000..0eaf6b00e958c5
--- /dev/null
+++ b/libc/src/math/exp2m1f.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for exp2m1f -----------------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_EXP2M1F_H
+#define LLVM_LIBC_SRC_MATH_EXP2M1F_H
+
+namespace LIBC_NAMESPACE {
+
+float exp2m1f(float x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_EXP2M1F_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 4d9b91150d0200..c1a4e5ea6e3b9d 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -837,6 +837,27 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  exp2m1f
+  SRCS
+    exp2m1f.cpp
+  HDRS
+    ../exp2m1f.h
+  DEPENDS
+    .explogxf
+    libc.src.errno.errno
+    libc.src.__support.common
+    libc.src.__support.FPUtil.fenv_impl
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
+    libc.src.__support.FPUtil.rounding_mode
+    libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.cpu_features
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   exp10
   SRCS
diff --git a/libc/src/math/generic/exp2m1f.cpp b/libc/src/math/generic/exp2m1f.cpp
new file mode 100644
index 00000000000000..3790171ed531eb
--- /dev/null
+++ b/libc/src/math/generic/exp2m1f.cpp
@@ -0,0 +1,182 @@
+//===-- Implementation of exp2m1f 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/exp2m1f.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/rounding_mode.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/errno/libc_errno.h"
+
+#include "explogxf.h"
+
+namespace LIBC_NAMESPACE {
+
+static constexpr size_t N_EXCEPTS_LO = 8;
+
+static constexpr fputil::ExceptValues<float, N_EXCEPTS_LO> EXP2M1F_EXCEPTS_LO =
+    {{
+        // (input, RZ output, RU offset, RD offset, RN offset)
+        // x = 0x1.36dc8ep-36, exp2m1f(x) = 0x1.aef212p-37 (RZ)
+        {0x2d9b'6e47U, 0x2d57'7909U, 1U, 0U, 0U},
+        // x = 0x1.224936p-19, exp2m1f(x) = 0x1.926c0ep-20 (RZ)
+        {0x3611'249bU, 0x35c9'3607U, 1U, 0U, 1U},
+        // x = 0x1.d16d2p-20, exp2m1f(x) = 0x1.429becp-20 (RZ)
+        {0x35e8'b690U, 0x35a1'4df6U, 1U, 0U, 1U},
+        // x = 0x1.17949ep-14, exp2m1f(x) = 0x1.8397p-15 (RZ)
+        {0x388b'ca4fU, 0x3841'cb80U, 1U, 0U, 1U},
+        // x = -0x1.9c3e1ep-38, exp2m1f(x) = -0x1.1dbeacp-38 (RZ)
+        {0xacce'1f0fU, 0xac8e'df56U, 0U, 1U, 0U},
+        // x = -0x1.4d89b4p-32, exp2m1f(x) = -0x1.ce61b6p-33 (RZ)
+        {0xafa6'c4daU, 0xaf67'30dbU, 0U, 1U, 1U},
+        // x = -0x1.a6eac4p-10, exp2m1f(x) = -0x1.24fadap-10 (RZ)
+        {0xbad3'7562U, 0xba92'7d6dU, 0U, 1U, 1U},
+        // x = -0x1.e7526ep-6, exp2m1f(x) = -0x1.4e53dep-6 (RZ)
+        {0xbcf3'a937U, 0xbca7'29efU, 0U, 1U, 1U},
+    }};
+
+static constexpr size_t N_EXCEPTS_HI = 2;
+
+static constexpr fputil::ExceptValues<float, N_EXCEPTS_HI> EXP2M1F_EXCEPTS_HI =
+    {{
+        // (input, RZ output, RU offset, RD offset, RN offset)
+        // x = 0x1.16a972p-1, exp2m1f(x) = 0x1.d545b2p-2 (RZ)
+        {0x3f0b'54b9U, 0x3eea'a2d9U, 1U, 0U, 0U},
+        // x = -0x1.9f12acp-5, exp2m1f(x) = -0x1.1ab68cp-5 (RZ)
+        {0xbd4f'8956U, 0xbd0d'5b46U, 0U, 1U, 0U},
+    }};
+
+LLVM_LIBC_FUNCTION(float, exp2m1f, (float x)) {
+  using FPBits = 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) {
+      if (auto r = EXP2M1F_EXCEPTS_LO.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+        return r.value();
+
+      // Minimax polynomial generated by Sollya with:
+      // > display = hexadecimal;
+      // > 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 = 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);
+      return static_cast<float>(p * xd);
+    }
+
+    // x >= 128, or x is nan
+    if (xbits.is_pos()) {
+      if (xbits.is_finite()) {
+        int rounding = fputil::quick_get_round();
+        if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+          return FPBits::max_normal().get_val();
+
+        fputil::set_errno_if_required(ERANGE);
+        fputil::raise_except_if_required(FE_OVERFLOW);
+      }
+
+      // x >= 128 and 2^x - 1 rounds to +inf, or x is +inf or nan
+      return x + FPBits::inf().get_val();
+    }
+  }
+
+  if (LIBC_UNLIKELY(x <= -25.0f)) {
+    // 2^(-inf) - 1 = -1
+    if (xbits.is_inf())
+      return -1.0f;
+    // 2^nan - 1 = nan
+    if (xbits.is_nan())
+      return x;
+
+    int rounding = fputil::quick_get_round();
+    if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO)
+      return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
+
+    fputil::set_errno_if_required(ERANGE);
+    fputil::raise_except_if_required(FE_UNDERFLOW);
+    return -1.0f;
+  }
+
+  if (auto r = EXP2M1F_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+    return r.value();
+
+  // For -25 < 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-4 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;
+  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);
+
+  // 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::FPBits<double>::FRACTION_LEN);
+  // 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;
+  double mh = fputil::FPBits<double>(static_cast<uint64_t>(mh_bits)).get_val();
+
+  // Degree-4 polynomial approximating (2^x - 1)/x generated by Sollya with:
+  // > display = hexadecimal;
+  // > fpminimax((2^x - 1)/x, 4, [|D...|], [-2^-6, 2^-6]);
+  constexpr double COEFFS[5] = {0x1.62e42fefa39efp-1, 0x1.ebfbdff8131c4p-3,
+                                0x1.c6b08d7061695p-5, 0x1.3b2b1bee74b2ap-7,
+                                0x1.5d88091198529p-10};
+  double dx_sq = dx * dx;
+  double c1 = fputil::multiply_add(dx, COEFFS[0], 1.0);
+  double c2 = fputil::multiply_add(dx, COEFFS[2], COEFFS[1]);
+  double c3 = fputil::multiply_add(dx, COEFFS[4], COEFFS[3]);
+  double p = fputil::multiply_add(dx_sq, c3, c2);
+  // 2^x = 2^(hi + mid + lo)
+  //     = 2^(hi + mid) * 2^lo
+  //     ~ mh * (1 + lo * P(lo))
+  //     = mh + (mh*lo) * P(lo)
+  double exp2_lo = fputil::multiply_add(p, dx_sq, c1);
+  return static_cast<float>(fputil::multiply_add(exp2_lo, mh, -1.0));
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 026bcd12928bdf..d27c3218416732 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -637,6 +637,22 @@ add_fp_unittest(
    libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  exp2m1f_test
+  NEED_MPFR
+  SUITE
+    libc-math-unittests
+  SRCS
+    exp2m1f_test.cpp
+  DEPENDS
+    libc.include.llvm-libc-macros.math_macros
+    libc.include.llvm-libc-macros.stdint_macros
+    libc.src.errno.errno
+    libc.src.math.exp2m1f
+    libc.src.__support.CPP.array
+    libc.src.__support.FPUtil.fp_bits
+)
+
 add_fp_unittest(
   exp10f_test
   NEED_MPFR
diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index df32dd4f943f30..6b2f3dddcadd24 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -142,6 +142,21 @@ add_fp_unittest(
     -lpthread
 )
 
+add_fp_unittest(
+  exp2m1f_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    exp2m1f_test.cpp
+  DEPENDS
+    .exhaustive_test
+    libc.src.math.exp2m1f
+  LINK_LIBRARIES
+    -lpthread
+)
+
 add_fp_unittest(
   exp10f_test
   NO_RUN_POSTBUILD
diff --git a/libc/test/src/math/exhaustive/exp2m1f_test.cpp b/libc/test/src/math/exhaustive/exp2m1f_test.cpp
new file mode 100644
index 00000000000000..2111024cb5c03b
--- /dev/null
+++ b/libc/test/src/math/exhaustive/exp2m1f_test.cpp
@@ -0,0 +1,33 @@
+//===-- Exhaustive test for exp2m1f ---------------------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/math/exp2m1f.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LlvmLibcExp2m1fExhaustiveTest =
+    LlvmLibcUnaryOpExhaustiveMathTest<float, mpfr::Operation::Exp2m1,
+                                      LIBC_NAMESPACE::exp2m1f>;
+
+// Range: [0, Inf];
+static constexpr uint32_t POS_START = 0x0000'0000U;
+static constexpr uint32_t POS_STOP = 0x7f80'0000U;
+
+TEST_F(LlvmLibcExp2m1fExhaustiveTest, PostiveRange) {
+  test_full_range_all_roundings(POS_START, POS_STOP);
+}
+
+// Range: [-Inf, 0];
+static constexpr uint32_t NEG_START = 0x8000'0000U;
+static constexpr uint32_t NEG_STOP = 0xff80'0000U;
+
+TEST_F(LlvmLibcExp2m1fExhaustiveTest, NegativeRange) {
+  test_full_range_all_roundings(NEG_START, NEG_STOP);
+}
diff --git a/libc/test/src/math/exp2m1f_test.cpp b/libc/test/src/math/exp2m1f_test.cpp
new file mode 100644
index 00000000000000..b35d94cb0f788d
--- /dev/null
+++ b/libc/test/src/math/exp2m1f_test.cpp
@@ -0,0 +1,67 @@
+//===-- Unittests for exp2m1f ---------------------------------------------===//
+//
+// 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 "include/llvm-libc-macros/math-macros.h"
+#include "src/__support/CPP/array.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/exp2m1f.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <stdint.h>
+
+using LlvmLibcExp2m1fTest = LIBC_NAMESPACE::testing::FPTest<float>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+TEST_F(LlvmLibcExp2m1fTest, TrickyInputs) {
+  constexpr LIBC_NAMESPACE::cpp::array<float, 10> INPUTS = {
+      // EXP2M1F_EXCEPTS_LO
+      0x1.36dc8ep-36,
+      0x1.224936p-19,
+      0x1.d16d2p-20,
+      0x1.17949ep-14,
+      -0x1.9c3e1ep-38,
+      -0x1.4d89b4p-32,
+      -0x1.a6eac4p-10,
+      -0x1.e7526ep-6,
+      // EXP2M1F_EXCEPTS_HI
+      0x1.16a972p-1,
+      -0x1.9f12acp-5,
+  };
+
+  for (float x : INPUTS) {
+    LIBC_NAMESPACE::libc_errno = 0;
+    EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2m1, x,
+                                   LIBC_NAMESPACE::exp2m1f(x), 0.5);
+    EXPECT_MATH_ERRNO(0);
+  }
+}
+
+TEST_F(LlvmLibcExp2m1fTest, InFloatRange) {
+  constexpr uint32_t COUNT = 100'000;
+  constexpr uint32_t STEP = UINT32_MAX / COUNT;
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = FPBits(v).get_val();
+    if (isnan(x) || isinf(x))
+      continue;
+    LIBC_NAMESPACE::libc_errno = 0;
+    float result = LIBC_NAMESPACE::exp2m1f(x);
+
+    // If the computation resulted in an error or did not produce valid result
+    // in the single-precision floating point range, then ignore comparing with
+    // MPFR result as MPFR can still produce valid results because of its
+    // wider precision.
+    if (isnan(result) || isinf(result) || LIBC_NAMESPACE::libc_errno != 0)
+      continue;
+    ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2m1, x,
+                                   LIBC_NAMESPACE::exp2m1f(x), 0.5);
+  }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 3b756127fe21ea..c65d6339d6c671 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -769,6 +769,17 @@ add_fp_unittest(
    libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  exp2m1f_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    exp2m1f_test.cpp
+  DEPENDS
+    libc.src.errno.errno
+    libc.src.math.exp2m1f
+)
+
 add_fp_unittest(
   exp10f_test
   SUITE
diff --git a/libc/test/src/math/smoke/exp2m1f_test.cpp b/libc/test/src/math/smoke/exp2m1f_test.cpp
new file mode 100644
index 00000000000000..c144856b30b769
--- /dev/null
+++ b/libc/test/src/math/smoke/exp2m1f_test.cpp
@@ -0,0 +1,68 @@
+//===-- Unittests for exp2m1f ---------------------------------------------===//
+//
+// 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/errno/libc_errno.h"
+#include "src/math/exp2m1f.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcExp2m1fTest = LIBC_NAMESPACE::testing::FPTest<float>;
+using LIBC_NAMESPACE::fputil::testing::ForceRoundingMode;
+using LIBC_NAMESPACE::fputil::testing::RoundingMode;
+
+TEST_F(LlvmLibcExp2m1fTest, SpecialNumbers) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::exp2m1f(aNaN));
+  EXPECT_MATH_ERRNO(0);
+  EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::exp2m1f(inf));
+  EXPECT_MATH_ERRNO(0);
+  EXPECT_FP_EQ_ALL_ROUNDING(-1.0f, LIBC_NAMESPACE::exp2m1f(neg_inf));
+  EXPECT_MATH_ERRNO(0);
+  EXPECT_FP_EQ_ALL_ROUNDING(0.0f, LIBC_NAMESPACE::exp2m1f(0.0f));
+  EXPECT_MATH_ERRNO(0);
+  EXPECT_FP_EQ_ALL_ROUNDING(-0.0f, LIBC_NAMESPACE::exp2m1f(-0.0f));
+  EXPECT_MATH_ERRNO(0);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(1.0f, LIBC_NAMESPACE::exp2m1f(1.0f));
+  EXPECT_FP_EQ_ALL_ROUNDING(-0.5f, LIBC_NAMESPACE::exp2m1f(-1.0f));
+  EXPECT_FP_EQ_ALL_ROUNDING(3.0f, LIBC_NAMESPACE::exp2m1f(2.0f));
+  EXPECT_FP_EQ_ALL_ROUNDING(-0.75f, LIBC_NAMESPACE::exp2m1f(-2.0f));
+}
+
+TEST_F(LlvmLibcExp2m1fTest, Overflow) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2m1f(0x1.fffffep+127),
+                              FE_OVERFLOW);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2m1f(128.0f),
+                              FE_OVERFLOW);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2m1f(0x1.000002p+7),
+                              FE_OVERFLOW);
+  EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST_F(LlvmLibcExp2m1fTest, Underflow) {
+  LIBC_NAMESPACE::libc_errno = 0;
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(-1.0f, LIBC_NAMESPACE::exp2m1f(-0x1.fffffep+127),
+                              FE_UNDERFLOW);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(-1.0f, LIBC_NAMESPACE::exp2m1f(-25.0f),
+                              FE_UNDERFLOW);
+  EXPECT_MATH_ERRNO(ERANGE);
+
+  EXPECT_FP_EQ_WITH_EXCEPTION(-1.0f, LIBC_NAMESPACE::exp2m1f(-0x1.900002p4),
+                              FE_UNDERFLOW);
+  EXPECT_MA...
[truncated]

@overmighty overmighty marked this pull request as ready for review April 1, 2024 14:38
{{
// (input, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.36dc8ep-36, exp2m1f(x) = 0x1.aef212p-37 (RZ)
{0x2d9b'6e47U, 0x2d57'7909U, 1U, 0U, 0U},
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure when exactly to use the U suffix and when not to use it.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

U suffix means the literal is unsigned. This is opposed to unsigned x = 7; where 7 is signed with an implicit cast inserted in front during AST parsing. Dumping the ast alludes to this, and makes it more understandable.

https://godbolt.org/z/7n5MT4xYM

Notice the declaration of foo has an ImplictCast node inserted in the AST; the declaration of bar does not.

Copy link
Contributor Author

@overmighty overmighty Apr 1, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes but existing code sometimes doesn't use the U suffix in literals that initialize variables of unsigned integer types and in literals that are compared against variables of unsigned integer types, and sometimes uses it. I'm not sure if there's a rule for this or if it just doesn't really matter if it's inconsistent (as long as it doesn't change the program's behavior).

libc/test/src/math/smoke/exp2m1f_test.cpp Outdated Show resolved Hide resolved
libc/utils/MPFRWrapper/MPFRUtils.cpp Outdated Show resolved Hide resolved
Comment on lines 173 to 178
double p = fputil::multiply_add(dx_sq, c3, c2);
// 2^x = 2^(hi + mid + lo)
// = 2^(hi + mid) * 2^lo
// ~ mh * (1 + lo * P(lo))
// = mh + (mh*lo) * P(lo)
double exp2_lo = fputil::multiply_add(p, dx_sq, c1);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can combine these into double p = fputil::polyeval(dx_sq, c1, c2, c3);
Also update the comment to reflect the computation steps are done here: mh * exp2_lo - 1.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't we also rename dx and dx_sq to lo and lo_sq instead of stating that dx = lo in the comment above the declaration of dx?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1. Also update dx appearance in some of the comments.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There aren't any comments in this file that reference dx.

#include "test/UnitTest/Test.h"
#include "utils/MPFRWrapper/MPFRUtils.h"

#include <stdint.h>
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought of including include/llvm-libc-macros/stdint-macros.h instead for UINT32_MAX, but there are other includes that end up including stdint.h, which redefines UINT32_MAX. I could include stdint-macros.h after all the other headers but I'm not sure if it makes sense.

Also, some files include *.h headers (e.g., stdint.h) and some include the c* alternatives (e.g., cstdint). Is there a rule for this?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's ok, we are going to re-organize and codify how to include standard headers soon: #87017

@lntue
Copy link
Contributor

lntue commented Apr 3, 2024

@overmighty : Can you sync and resolve the merge conflict?

@overmighty
Copy link
Contributor Author

overmighty commented Apr 3, 2024

Sure. Sorry I didn't notice it earlier.

@lntue
Copy link
Contributor

lntue commented Apr 3, 2024

@zimmermann6

@zimmermann6
Copy link

I confirm all exhaustive tests are ok. On a Intel(R) Xeon(R) Silver 4214 I get the following timings:

GNU libc version: 2.37
GNU libc release: stable
[####################] 100 %
Ntrial = 20 ; Min = 11.306 + 0.522 clc/call; Median-Min = 0.384 clc/call; Max = 12.291 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 11.572 + 0.585 clc/call; Median-Min = 0.536 clc/call; Max = 12.625 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 11.416 + 0.552 clc/call; Median-Min = 0.504 clc/call; Max = 12.269 clc/call;
GNU libc version: 2.37
GNU libc release: stable
[####################] 100 %
Ntrial = 20 ; Min = 37.811 + 0.627 clc/call; Median-Min = 0.580 clc/call; Max = 39.107 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 37.707 + 0.430 clc/call; Median-Min = 0.346 clc/call; Max = 38.526 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 37.525 + 0.654 clc/call; Median-Min = 0.372 clc/call; Max = 39.523 clc/call;

The timings are very close to core-math (the "glibc" timings are those from core-math, since exp2m1f is not yet available in GNU libc).

@lntue
Copy link
Contributor

lntue commented Apr 4, 2024

I confirm all exhaustive tests are ok. On a Intel(R) Xeon(R) Silver 4214 I get the following timings:

GNU libc version: 2.37
GNU libc release: stable
[####################] 100 %
Ntrial = 20 ; Min = 11.306 + 0.522 clc/call; Median-Min = 0.384 clc/call; Max = 12.291 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 11.572 + 0.585 clc/call; Median-Min = 0.536 clc/call; Max = 12.625 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 11.416 + 0.552 clc/call; Median-Min = 0.504 clc/call; Max = 12.269 clc/call;
GNU libc version: 2.37
GNU libc release: stable
[####################] 100 %
Ntrial = 20 ; Min = 37.811 + 0.627 clc/call; Median-Min = 0.580 clc/call; Max = 39.107 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 37.707 + 0.430 clc/call; Median-Min = 0.346 clc/call; Max = 38.526 clc/call;
[####################] 100 %
Ntrial = 20 ; Min = 37.525 + 0.654 clc/call; Median-Min = 0.372 clc/call; Max = 39.523 clc/call;

The timings are very close to core-math (the "glibc" timings are those from core-math, since exp2m1f is not yet available in GNU libc).

Thanks Paul!

@lntue lntue merged commit a8c5975 into llvm:main Apr 4, 2024
5 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[libc][math][c23] Implement C23 math function exp2m1f
5 participants