[APFloat] Add exp function for APFloat::IEEESsingle using expf implementation from LLVM libc. (#143959)

Discourse RFC:
https://discourse.llvm.org/t/rfc-make-clang-builtin-math-functions-constexpr-with-llvm-libc-to-support-c-23-constexpr-math-functions/86450

- The implementation in LLVM libc is header-only.
- expf implementation in LLVM libc is correctly rounded for all rounding
modes.
- LLVM libc implementation will round to the floating point
environment's rounding mode.
- No cmake build dependency between LLVM and LLVM libc, only requires
LLVM libc source presents in llvm-project/libc folder.

Author: lntue

libc/shared/math/expf.h

@@ -11,6 +11,7 @@
 
 #include "shared/libc_common.h"
 #include "src/__support/math/expf.h"
+#include "src/__support/math/expf_static_rounding.h"
 
 namespace LIBC_NAMESPACE_DECL {
 namespace shared {

libc/src/__support/math/expf_static_rounding.h

@@ -0,0 +1,43 @@
+//===-- Implementation header for expf --------------------------*- 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___SUPPORT_MATH_EXPF_STATIC_ROUNDING_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_STATIC_ROUNDING_H
+
+#include "expf.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+// Remark: "#pragma STDC FENV_ACCESS" might generate unsupported warnings on
+// certain platforms.
+#pragma STDC FENV_ACCESS ON
+
+// Directional rounding version of expf.
+LIBC_INLINE static float expf(float x, int rounding_mode) {
+  int current_rounding_mode = fputil::get_round();
+  if (rounding_mode == current_rounding_mode)
+    return expf(x);
+
+  fputil::set_round(rounding_mode);
+  float result = expf(x);
+  fputil::set_round(current_rounding_mode);
+  return result;
+}
+
+#pragma STDC FENV_ACCESS DEFAULT
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_STATIC_ROUNDING_H

llvm/CMakeLists.txt

@@ -657,6 +657,11 @@ endif()
 
 set(LLVM_ENABLE_Z3_SOLVER_DEFAULT "${Z3_FOUND}")
 
+include(FindLibcCommonUtils)
+if(NOT TARGET llvm-libc-common-utilities)
+  message(FATAL_ERROR "LLVM libc is not found at ${libc_path}.")
+endif()
+
 
 if( LLVM_TARGETS_TO_BUILD STREQUAL "all" )
   set( LLVM_TARGETS_TO_BUILD ${LLVM_ALL_TARGETS} )

llvm/include/llvm/ADT/APFloat.h

@@ -1493,6 +1493,8 @@ class APFloat : public APFloatBase {
   friend APFloat frexp(const APFloat &X, int &Exp, roundingMode RM);
   friend IEEEFloat;
   friend DoubleAPFloat;
+
+  friend APFloat exp(const APFloat &X, roundingMode RM);
 };
 
 static_assert(sizeof(APFloat) == sizeof(detail::IEEEFloat),
@@ -1645,6 +1647,10 @@ inline APFloat maximumnum(const APFloat &A, const APFloat &B) {
   return A < B ? B : A;
 }
 
+/// Implement IEEE 754-2019 exp functions.
+LLVM_READONLY
+APFloat exp(const APFloat &X, RoundingMode RM = APFloat::rmNearestTiesToEven);
+
 inline raw_ostream &operator<<(raw_ostream &OS, const APFloat &V) {
   V.print(OS);
   return OS;

llvm/lib/Support/APFloat.cpp

@@ -28,6 +28,9 @@
 #include <cstring>
 #include <limits.h>
 
+// Shared headers from LLVM libc
+#include "shared/math.h"
+
 #define APFLOAT_DISPATCH_ON_SEMANTICS(METHOD_CALL)                             \
   do {                                                                         \
     if (usesLayout<IEEEFloat>(getSemantics()))                                 \
@@ -6155,6 +6158,29 @@ float APFloat::convertToFloat() const {
   return Temp.getIEEE().convertToFloat();
 }
 
+static constexpr int getFEnvRoundingMode(llvm::RoundingMode rm) {
+  switch (rm) {
+  case APFloat::rmTowardPositive:
+    return FE_UPWARD;
+  case APFloat::rmTowardNegative:
+    return FE_DOWNWARD;
+  case APFloat::rmTowardZero:
+    return FE_TOWARDZERO;
+  default:
+    // TODO: fix rmNearestTiesToAway for platform without FE_TONEARESTFROMZERO.
+    return FE_TONEAREST;
+  };
+}
+
+APFloat exp(const APFloat &X, RoundingMode rounding_mode) {
+  if (&X.getSemantics() == &APFloat::IEEEsingle()) {
+    float result = LIBC_NAMESPACE::shared::expf(
+        X.convertToFloat(), getFEnvRoundingMode(rounding_mode));
+    return APFloat(result);
+  }
+  llvm_unreachable("Unsupported floating-point semantics");
+}
+
 APFloat::Storage::~Storage() {
   if (usesLayout<IEEEFloat>(*semantics)) {
     IEEE.~IEEEFloat();

llvm/lib/Support/CMakeLists.txt

@@ -395,3 +395,9 @@ target_include_directories(LLVMSupport
   PRIVATE
   ${LLVM_THIRD_PARTY_DIR}/siphash/include
 )
+
+# Integrating LLVM libc's math functions
+target_include_directories(LLVMSupport PRIVATE "${LLVM_INCLUDE_DIR}/../../libc")
+if(NOT MSVC)
+  target_compile_options(LLVMSupport PRIVATE "-Wno-c99-extensions") # _Complex warnings.
+endif()

llvm/lib/Target/AMDGPU/AMDGPULibCalls.cpp

@@ -1531,7 +1531,7 @@ bool AMDGPULibCalls::evaluateScalarMathFunc(const FuncInfo &FInfo, double &Res0,
     return true;
 
   case AMDGPULibFunc::EI_EXP:
-    Res0 = exp(opr0);
+    Res0 = std::exp(opr0);
     return true;
 
   case AMDGPULibFunc::EI_EXP2:

llvm/unittests/ADT/APFloatTest.cpp

@@ -10182,4 +10182,78 @@ TEST(APFloatTest, FrexpQuietSNaN) {
   EXPECT_FALSE(Result.isSignaling());
 }
 
+TEST(APFloatTest, expf) {
+  std::array<llvm::RoundingMode, 4> allRoundingModes = {
+      APFloat::rmNearestTiesToEven, APFloat::rmTowardPositive,
+      APFloat::rmTowardNegative, APFloat::rmTowardZero};
+  for (auto rm : allRoundingModes) {
+    // exp(+-0) = 1 for all rounding modes.
+    EXPECT_EQ(1.0f, llvm::exp(APFloat(0.0f), rm).convertToFloat());
+    EXPECT_EQ(1.0f, llvm::exp(APFloat(-0.0f), rm).convertToFloat());
+    // exp(+Inf) = +Inf for all rounding modes.
+    EXPECT_EQ(std::numeric_limits<float>::infinity(),
+              llvm::exp(APFloat::getInf(APFloat::IEEEsingle(), false), rm)
+                  .convertToFloat());
+    // exp(-Inf) = 0 for all rounding modes.
+    EXPECT_EQ(0.0f, llvm::exp(APFloat::getInf(APFloat::IEEEsingle(), true), rm)
+                        .convertToFloat());
+    // exp(NaN) = NaN for all rounding modes.
+    EXPECT_TRUE(llvm::exp(APFloat::getNaN(APFloat::IEEEsingle()), rm).isNaN());
+  }
+  // exp(1)
+  EXPECT_EQ(
+      0x1.5bf0a8p1f,
+      llvm::exp(APFloat(1.0f), APFloat::rmNearestTiesToEven).convertToFloat());
+  EXPECT_EQ(
+      0x1.5bf0aap1f,
+      llvm::exp(APFloat(1.0f), APFloat::rmTowardPositive).convertToFloat());
+  EXPECT_EQ(
+      0x1.5bf0a8p1f,
+      llvm::exp(APFloat(1.0f), APFloat::rmTowardNegative).convertToFloat());
+  EXPECT_EQ(0x1.5bf0a8p1f,
+            llvm::exp(APFloat(1.0f), APFloat::rmTowardZero).convertToFloat());
+  // exp(float max)
+  EXPECT_EQ(std::numeric_limits<float>::infinity(),
+            llvm::exp(APFloat::getLargest(APFloat::IEEEsingle(), false),
+                      APFloat::rmNearestTiesToEven)
+                .convertToFloat());
+  EXPECT_EQ(std::numeric_limits<float>::infinity(),
+            llvm::exp(APFloat::getLargest(APFloat::IEEEsingle(), false),
+                      APFloat::rmTowardPositive)
+                .convertToFloat());
+  EXPECT_EQ(std::numeric_limits<float>::max(),
+            llvm::exp(APFloat::getLargest(APFloat::IEEEsingle(), false),
+                      APFloat::rmTowardNegative)
+                .convertToFloat());
+  EXPECT_EQ(std::numeric_limits<float>::max(),
+            llvm::exp(APFloat::getLargest(APFloat::IEEEsingle(), false),
+                      APFloat::rmTowardZero)
+                .convertToFloat());
+  // exp(min_denormal)
+  EXPECT_EQ(1.0f, llvm::exp(APFloat::getSmallest(APFloat::IEEEsingle(), false),
+                            APFloat::rmNearestTiesToEven)
+                      .convertToFloat());
+  EXPECT_EQ(0x1.000002p0f,
+            llvm::exp(APFloat::getSmallest(APFloat::IEEEsingle(), false),
+                      APFloat::rmTowardPositive)
+                .convertToFloat());
+  EXPECT_EQ(1.0f, llvm::exp(APFloat::getSmallest(APFloat::IEEEsingle(), false),
+                            APFloat::rmTowardNegative)
+                      .convertToFloat());
+  EXPECT_EQ(1.0f, llvm::exp(APFloat::getSmallest(APFloat::IEEEsingle(), false),
+                            APFloat::rmTowardZero)
+                      .convertToFloat());
+  // Default rounding mode.
+  // exp(-1)
+  EXPECT_EQ(0x1.78b564p-2f, llvm::exp(APFloat(-1.0f)).convertToFloat());
+  EXPECT_EQ(
+      0x1.78b564p-2f,
+      llvm::exp(APFloat(-1.0f), APFloat::rmTowardPositive).convertToFloat());
+  EXPECT_EQ(
+      0x1.78b562p-2f,
+      llvm::exp(APFloat(-1.0f), APFloat::rmTowardNegative).convertToFloat());
+  EXPECT_EQ(0x1.78b562p-2f,
+            llvm::exp(APFloat(-1.0f), APFloat::rmTowardZero).convertToFloat());
+}
+
 } // namespace