Skip to content

Commit a5466f0

Browse files
lntuelntue
committed
[libc] Improve the performance of expm1f.
Improve the performance of expm1f: - Rearrange the selection logic for different cases to improve the overall throughput. - Use the same degree-4 polynomial for large inputs as `expf` (https://reviews.llvm.org/D122418), reduced from a degree-7 polynomial. Performance benchmark using perf tool from CORE-MATH project (https://gitlab.inria.fr/core-math/core-math/-/tree/master): Before this patch: ``` $ ./perf.sh expm1f CORE-MATH reciprocal throughput : 15.362 System LIBC reciprocal throughput : 53.288 LIBC reciprocal throughput : 54.572 $ ./perf.sh expm1f --latency CORE-MATH latency : 57.759 System LIBC latency : 147.146 LIBC latency : 118.057 ``` After this patch: ``` $ ./perf.sh expm1f CORE-MATH reciprocal throughput : 15.359 System LIBC reciprocal throughput : 53.188 LIBC reciprocal throughput : 14.600 $ ./perf.sh expm1f --latency CORE-MATH latency : 57.774 System LIBC latency : 147.119 LIBC latency : 60.280 ``` Reviewed By: michaelrj, santoshn, zimmermann6 Differential Revision: https://reviews.llvm.org/D122538
1 parent f6740fe commit a5466f0

File tree

4 files changed

+82
-60
lines changed

4 files changed

+82
-60
lines changed

‎libc/src/math/generic/expm1f.cpp

Lines changed: 56 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -24,41 +24,59 @@ LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
2424
using FPBits = typename fputil::FPBits<float>;
2525
FPBits xbits(x);
2626

27-
// When x < log(2^-25) or nan
28-
if (unlikely(xbits.uintval() >= 0xc18a'a123U)) {
29-
// exp(-Inf) = 0
30-
if (xbits.is_inf())
31-
return -1.0f;
32-
// exp(nan) = nan
33-
if (xbits.is_nan())
34-
return x;
27+
uint32_t x_u = xbits.uintval();
28+
uint32_t x_abs = x_u & 0x7fff'ffffU;
29+
30+
// Exceptional value
31+
if (unlikely(x_u == 0x3e35'bec5U)) { // x = 0x1.6b7d8ap-3f
3532
int round_mode = fputil::get_round();
36-
if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO)
37-
return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
38-
return -1.0f;
33+
if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
34+
return 0x1.8dbe64p-3f;
35+
return 0x1.8dbe62p-3f;
3936
}
40-
// x >= 89 or nan
41-
if (unlikely(!xbits.get_sign() && (xbits.uintval() >= 0x42b2'0000))) {
42-
if (xbits.uintval() < 0x7f80'0000U) {
43-
int rounding = fputil::get_round();
44-
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
45-
return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
4637

47-
errno = ERANGE;
38+
// When |x| > 25*log(2), or nan
39+
if (unlikely(x_abs >= 0x418a'a123U)) {
40+
// x < log(2^-25)
41+
if (xbits.get_sign()) {
42+
// exp(-Inf) = 0
43+
if (xbits.is_inf())
44+
return -1.0f;
45+
// exp(nan) = nan
46+
if (xbits.is_nan())
47+
return x;
48+
int round_mode = fputil::get_round();
49+
if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO)
50+
return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
51+
return -1.0f;
52+
} else {
53+
// x >= 89 or nan
54+
if (xbits.uintval() >= 0x42b2'0000) {
55+
if (xbits.uintval() < 0x7f80'0000U) {
56+
int rounding = fputil::get_round();
57+
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
58+
return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
59+
60+
errno = ERANGE;
61+
}
62+
return x + static_cast<float>(FPBits::inf());
63+
}
4864
}
49-
return x + static_cast<float>(FPBits::inf());
5065
}
5166

52-
int unbiased_exponent = static_cast<int>(xbits.get_unbiased_exponent());
5367
// |x| < 2^-4
54-
if (unbiased_exponent < 123) {
68+
if (x_abs < 0x3d80'0000U) {
5569
// |x| < 2^-25
56-
if (unbiased_exponent < 102) {
70+
if (x_abs < 0x3300'0000U) {
5771
// x = -0.0f
5872
if (unlikely(xbits.uintval() == 0x8000'0000U))
5973
return x;
60-
// When |x| < 2^-25, the relative error:
61-
// |(e^x - 1) - x| / |x| < |x^2| / |x| = |x| < 2^-25 < epsilon(1)/2.
74+
// When |x| < 2^-25, the relative error of the approximation e^x - 1 ~ x
75+
// is:
76+
// |(e^x - 1) - x| / |e^x - 1| < |x^2| / |x|
77+
// = |x|
78+
// < 2^-25
79+
// < epsilon(1)/2.
6280
// So the correctly rounded values of expm1(x) are:
6381
// = x + eps(x) if rounding mode = FE_UPWARD,
6482
// or (rounding mode = FE_TOWARDZERO and x is negative),
@@ -67,20 +85,21 @@ LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
6785
// fma(x, x, x) ~ x + x^2 instead.
6886
return fputil::fma(x, x, x);
6987
}
88+
7089
// 2^-25 <= |x| < 2^-4
7190
double xd = static_cast<double>(x);
7291
double xsq = xd * xd;
7392
// Degree-8 minimax polynomial generated by Sollya with:
7493
// > display = hexadecimal;
75-
// > P = fpminimax(expm1(x)/x, 7, [|D...|], [-2^-4, 2^-4]);
94+
// > P = fpminimax((expm1(x) - x)/x^2, 6, [|D...|], [-2^-4, 2^-4]);
7695
double r =
77-
fputil::polyeval(xd, 0x1p-1, 0x1.55555555559abp-3, 0x1.55555555551a7p-5,
78-
0x1.111110f70f2a4p-7, 0x1.6c16c17639e82p-10,
79-
0x1.a02526febbea6p-13, 0x1.a01dc40888fcdp-16);
96+
fputil::polyeval(xd, 0x1p-1, 0x1.55555555557ddp-3, 0x1.55555555552fap-5,
97+
0x1.111110fcd58b7p-7, 0x1.6c16c1717660bp-10,
98+
0x1.a0241f0006d62p-13, 0x1.a01e3f8d3c06p-16);
8099
return static_cast<float>(fputil::fma(r, xsq, xd));
81100
}
82101

83-
// For -18 < x < 89, to compute exp(x), we perform the following range
102+
// For -18 < x < 89, to compute expm1(x), we perform the following range
84103
// reduction: find hi, mid, lo such that:
85104
// x = hi + mid + lo, in which
86105
// hi is an integer,
@@ -89,48 +108,30 @@ LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
89108
// In particular,
90109
// hi + mid = round(x * 2^7) * 2^(-7).
91110
// Then,
92-
// exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).
111+
// expm1(x) = exp(hi + mid + lo) - 1 = exp(hi) * exp(mid) * exp(lo) - 1.
93112
// We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2
94-
// respectively. exp(lo) is computed using a degree-7 minimax polynomial
113+
// respectively. exp(lo) is computed using a degree-4 minimax polynomial
95114
// generated by Sollya.
96115

97-
// Exceptional value
98-
if (xbits.uintval() == 0xbdc1'c6cbU) {
99-
// x = -0x1.838d96p-4f
100-
int round_mode = fputil::get_round();
101-
if (round_mode == FE_TONEAREST || round_mode == FE_DOWNWARD)
102-
return -0x1.71c884p-4f;
103-
return -0x1.71c882p-4f;
104-
}
105-
106116
// x_hi = hi + mid.
107-
int x_hi = static_cast<int>(x * 0x1.0p7f);
117+
int x_hi = static_cast<int>(x * 0x1.0p7f + (xbits.get_sign() ? -0.5f : 0.5f));
108118
// Subtract (hi + mid) from x to get lo.
109119
x -= static_cast<float>(x_hi) * 0x1.0p-7f;
110120
double xd = static_cast<double>(x);
111-
// Make sure that -2^(-8) <= lo < 2^-8.
112-
if (x >= 0x1.0p-8f) {
113-
++x_hi;
114-
xd -= 0x1.0p-7;
115-
}
116-
if (x < -0x1.0p-8f) {
117-
--x_hi;
118-
xd += 0x1.0p-7;
119-
}
120121
x_hi += 104 << 7;
121122
// hi = x_hi >> 7
122123
double exp_hi = EXP_M1[x_hi >> 7];
123124
// lo = x_hi & 0x0000'007fU;
124125
double exp_mid = EXP_M2[x_hi & 0x7f];
125126
double exp_hi_mid = exp_hi * exp_mid;
126-
// Degree-7 minimax polynomial generated by Sollya with the following
127+
// Degree-4 minimax polynomial generated by Sollya with the following
127128
// commands:
128129
// > display = hexadecimal;
129-
// > Q = fpminimax(expm1(x)/x, 6, [|D...|], [-2^-8, 2^-8]);
130+
// > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]);
130131
// > Q;
131-
double exp_lo = fputil::polyeval(
132-
xd, 0x1p0, 0x1p0, 0x1p-1, 0x1.5555555555555p-3, 0x1.55555555553ap-5,
133-
0x1.1111111204dfcp-7, 0x1.6c16cb2da593ap-10, 0x1.9ff1648996d2ep-13);
132+
double exp_lo =
133+
fputil::polyeval(xd, 0x1.0p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1,
134+
0x1.555566668e5e7p-3, 0x1.55555555ef243p-5);
134135
return static_cast<float>(fputil::fma(exp_hi_mid, exp_lo, -1.0));
135136
}
136137

‎libc/test/src/math/exhaustive/exhaustive_test.cpp

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include <vector>
1616

1717
#include "exhaustive_test.h"
18+
#include "src/__support/FPUtil/FPBits.h"
1819
#include "utils/UnitTest/Test.h"
1920

2021
template <typename T>
@@ -28,14 +29,26 @@ void LlvmLibcExhaustiveTest<T>::test_full_range(T start, T stop, int nthreads,
2829
thread_list.emplace_back([this, begin, end, rounding]() {
2930
std::stringstream msg;
3031
msg << "-- Testing from " << begin << " to " << end << " [0x" << std::hex
31-
<< begin << ", 0x" << end << ") ..." << std::endl;
32+
<< begin << ", 0x" << end << "), [" << std::hexfloat
33+
<< float(__llvm_libc::fputil::FPBits<float>(
34+
static_cast<uint32_t>(begin)))
35+
<< ", "
36+
<< float(
37+
__llvm_libc::fputil::FPBits<float>(static_cast<uint32_t>(end)))
38+
<< ") ..." << std::endl;
3239
std::cout << msg.str();
3340
msg.str("");
3441

3542
bool result = check(begin, end, rounding);
3643

3744
msg << "** Finished testing from " << std::dec << begin << " to " << end
38-
<< " [0x" << std::hex << begin << ", 0x" << end
45+
<< " [0x" << std::hex << begin << ", 0x" << end << "), ["
46+
<< std::hexfloat
47+
<< float(__llvm_libc::fputil::FPBits<float>(
48+
static_cast<uint32_t>(begin)))
49+
<< ", "
50+
<< float(
51+
__llvm_libc::fputil::FPBits<float>(static_cast<uint32_t>(end)))
3952
<< ") : " << (result ? "PASSED" : "FAILED") << std::endl;
4053
std::cout << msg.str();
4154
});

‎libc/test/src/math/exhaustive/expm1f_test.cpp

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@
1212
#include "utils/MPFRWrapper/MPFRUtils.h"
1313
#include "utils/UnitTest/FPMatcher.h"
1414

15+
#include <thread>
16+
1517
using FPBits = __llvm_libc::fputil::FPBits<float>;
1618

1719
namespace mpfr = __llvm_libc::testing::mpfr;
@@ -20,18 +22,19 @@ struct LlvmLibcExpfExhaustiveTest : public LlvmLibcExhaustiveTest<uint32_t> {
2022
bool check(uint32_t start, uint32_t stop,
2123
mpfr::RoundingMode rounding) override {
2224
mpfr::ForceRoundingMode r(rounding);
23-
uint32_t bits = start;
25+
uint32_t bits = stop;
2426
bool result = true;
2527
do {
2628
FPBits xbits(bits);
2729
float x = float(xbits);
2830
result &= EXPECT_MPFR_MATCH(mpfr::Operation::Expm1, x,
2931
__llvm_libc::expm1f(x), 0.5, rounding);
30-
} while (bits++ < stop);
32+
} while (bits-- > start);
33+
return result;
3134
}
3235
};
3336

34-
static constexpr int NUM_THREADS = 16;
37+
static const int NUM_THREADS = std::thread::hardware_concurrency();
3538

3639
// Range: [0, 89];
3740
static constexpr uint32_t POS_START = 0x0000'0000U;

‎libc/test/src/math/expm1f_test.cpp

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,11 @@ TEST(LlvmLibcExpm1fTest, Borderline) {
9292
ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x,
9393
__llvm_libc::expm1f(x), 0.5);
9494
EXPECT_MATH_ERRNO(0);
95+
96+
x = float(FPBits(0x3e35bec5U));
97+
ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x,
98+
__llvm_libc::expm1f(x), 0.5);
99+
EXPECT_MATH_ERRNO(0);
95100
}
96101

97102
TEST(LlvmLibcExpm1fTest, InFloatRange) {

0 commit comments

Comments
 (0)