Skip to content
Open
Show file tree
Hide file tree
Changes from 11 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions libc/src/__support/math/expf.h
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,22 @@ LIBC_INLINE static constexpr float expf(float x) {
return static_cast<float>(exp_hi * exp_mid * exp_lo);
}

#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();
Copy link
Contributor

Choose a reason for hiding this comment

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

Use the pragma to enable fenv_access?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Done.

Copy link
Contributor

Choose a reason for hiding this comment

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

I would still be more comfortable if you just add an exp implementation with a template argument for the rounding mode, and didn't use fenv access anywhere in the stack. This has just merely shuffled the problem into the libc code directly

Copy link
Contributor Author

Choose a reason for hiding this comment

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

At some point I'll implement a version of directed rounding math functions without relying on adjusting the floating point environment.

Copy link
Contributor

Choose a reason for hiding this comment

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

I'm kind of surprised that setting the rounding mode for a foo() implementation results in a foo() result with the correct rounding mode. Are the transcendental functions in libc specifically implemented in a way that satisfies that property?

Can you point me to where libc is testing these functions in different rounding modes?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@nikic: Yes, at the moment, all LLVM libc transcendental math functions are rounded toward the current (dynamic) floating point environment rounding mode (https://libc.llvm.org/headers/math/index.html#implementation-requirements-goals), and we test them for all 4 rounding modes (FE_NEAREST, FE_UPWARD, FE_DOWNWARD, and FE_TOWARDZERO) against MPFR outputs.

Some examples in our test suites:

Copy link
Contributor

Choose a reason for hiding this comment

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

Awesome, thanks for the references!

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
Expand Down
5 changes: 5 additions & 0 deletions llvm/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -659,6 +659,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} )
Expand Down
6 changes: 6 additions & 0 deletions llvm/include/llvm/ADT/APFloat.h
Original file line number Diff line number Diff line change
Expand Up @@ -1525,6 +1525,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),
Expand Down Expand Up @@ -1658,6 +1660,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;
Expand Down
26 changes: 26 additions & 0 deletions llvm/lib/Support/APFloat.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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())) \
Expand Down Expand Up @@ -5601,6 +5604,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() == &semIEEEsingle) {
float result = LIBC_NAMESPACE::shared::expf(
X.convertToFloat(), getFEnvRoundingMode(rounding_mode));
return APFloat(result);
Comment on lines +6175 to +6177
Copy link
Contributor

Choose a reason for hiding this comment

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

Won't this provide different results depending on whether or not the platform we are running against uses DAZ or FTZ?

It seems a bit confusing to rely on host numerics in a function which goes from APFloat -> APFloat. I would recommend that APFloat just have its own copy of the implementation, that way the result of the computation won't depend on the host's ambient conditions.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If someone compiles or runs clang/llvm with FTZ/DAZ, do they really expect the results related to floating point denormals would be consistent with clang/llvm without FTZ/DAZ? Even with a perfect implementation of these functions using entirely integer arithmetics, every time a convertToDouble() or convertToFloat() are called then all bets are off.

The implementations we have for these math functions are as consistent as we could, and the unit tests added in this PRs also guard against being compiled or run with FTZ/DAZ mode. If some users really want to build or run clang/llvm with FTZ/DAZ for specific targets, we could in theory perform a clear-and-restore those flags in the directed rounding functions. But there will be a lot more places that they need to worry about in that case.

I don't think even MPFR could guarantee correct rounding when building/running under FTZ/DAZ flags.

Copy link
Contributor

Choose a reason for hiding this comment

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

It seems a bit confusing to rely on host numerics in a function which goes from APFloat -> APFloat. I would recommend that APFloat just have its own copy of the implementation, that way the result of the computation won't depend on the host's ambient conditions.

I'm somewhat concerned about this too. But I've already manually verified that the proposed tests for this method will fail if DAZ/FTZ are enabled, and given that the current state of affairs is that we just blindly call host exp or similar for constant folding, this is at least a strict improvement over the current state of affairs

Personally, I'd think it best if the implementation just worked off of APFloat values directly, which it likely will have to for ppc_fp128, x86_fp80, and fp128 types anyways.

}
assert(false && "Unsupported floating-point semantics");
}

} // namespace llvm

#undef APFLOAT_DISPATCH_ON_SEMANTICS
4 changes: 4 additions & 0 deletions llvm/lib/Support/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -384,3 +384,7 @@ 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")
target_compile_options(LLVMSupport PRIVATE "-Wno-c99-extensions") # _Complex warnings.
2 changes: 1 addition & 1 deletion llvm/lib/Target/AMDGPU/AMDGPULibCalls.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1540,7 +1540,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:
Expand Down
75 changes: 75 additions & 0 deletions llvm/unittests/ADT/APFloatTest.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "llvm/Support/FormatVariadic.h"
#include "gtest/gtest.h"
#include <cmath>
#include <limits>
#include <ostream>
#include <string>
#include <tuple>
Expand Down Expand Up @@ -8356,4 +8357,78 @@ TEST(APFloatTest, hasSignBitInMSB) {
EXPECT_FALSE(APFloat::hasSignBitInMSB(APFloat::Float8E8M0FNU()));
}

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
Loading