Skip to content

Commit 565230a

Browse files
overmightyovermighty
committed
fixup! [libc][math][c23] Add exp2m1f C23 math function
1 parent a813c1f commit 565230a

File tree

4 files changed

+42
-36
lines changed

4 files changed

+42
-36
lines changed

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

Lines changed: 18 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,7 @@ static constexpr fputil::ExceptValues<float, N_EXCEPTS_LO> EXP2M1F_EXCEPTS_LO =
4545
{0xbcf3'a937U, 0xbca7'29efU, 0U, 1U, 1U},
4646
}};
4747

48-
static constexpr size_t N_EXCEPTS_HI = 2;
48+
static constexpr size_t N_EXCEPTS_HI = 3;
4949

5050
static constexpr fputil::ExceptValues<float, N_EXCEPTS_HI> EXP2M1F_EXCEPTS_HI =
5151
{{
@@ -54,6 +54,8 @@ static constexpr fputil::ExceptValues<float, N_EXCEPTS_HI> EXP2M1F_EXCEPTS_HI =
5454
{0x3f0b'54b9U, 0x3eea'a2d9U, 1U, 0U, 0U},
5555
// x = -0x1.9f12acp-5, exp2m1f(x) = -0x1.1ab68cp-5 (RZ)
5656
{0xbd4f'8956U, 0xbd0d'5b46U, 0U, 1U, 0U},
57+
// x = -0x1.de7b9cp-5, exp2m1f(x) = -0x1.4508f4p-5 (RZ)
58+
{0xbd6f'3dceU, 0xbd22'847aU, 0U, 1U, 1U},
5759
}};
5860

5961
LLVM_LIBC_FUNCTION(float, exp2m1f, (float x)) {
@@ -66,7 +68,7 @@ LLVM_LIBC_FUNCTION(float, exp2m1f, (float x)) {
6668
// When |x| >= 128, or x is nan, or |x| <= 2^-5
6769
if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U || x_abs <= 0x3d00'0000U)) {
6870
// |x| <= 2^-5
69-
if (x_abs <= 0x3d00'0000) {
71+
if (x_abs <= 0x3d00'0000U) {
7072
if (auto r = EXP2M1F_EXCEPTS_LO.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
7173
return r.value();
7274

@@ -147,17 +149,17 @@ LLVM_LIBC_FUNCTION(float, exp2m1f, (float x)) {
147149
kf = static_cast<float>(k);
148150
#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT
149151
150-
// dx = lo = x - (hi + mid) = x - kf * 2^(-5)
151-
double dx = fputil::multiply_add(-0x1.0p-5f, kf, x);
152+
// lo = x - (hi + mid) = x - kf * 2^(-5)
153+
double lo = fputil::multiply_add(-0x1.0p-5f, kf, x);
152154
153155
// hi = floor(kf * 2^(-4))
154-
// exp_hi = shift hi to the exponent field of double precision.
155-
int64_t exp_hi =
156+
// exp2_hi = shift hi to the exponent field of double precision.
157+
int64_t exp2_hi =
156158
static_cast<int64_t>(static_cast<uint64_t>(k >> ExpBase::MID_BITS)
157159
<< fputil::FPBits<double>::FRACTION_LEN);
158160
// mh = 2^hi * 2^mid
159161
// mh_bits = bit field of mh
160-
int64_t mh_bits = ExpBase::EXP_2_MID[k & ExpBase::MID_MASK] + exp_hi;
162+
int64_t mh_bits = ExpBase::EXP_2_MID[k & ExpBase::MID_MASK] + exp2_hi;
161163
double mh = fputil::FPBits<double>(static_cast<uint64_t>(mh_bits)).get_val();
162164
163165
// Degree-4 polynomial approximating (2^x - 1)/x generated by Sollya with:
@@ -166,16 +168,15 @@ LLVM_LIBC_FUNCTION(float, exp2m1f, (float x)) {
166168
constexpr double COEFFS[5] = {0x1.62e42fefa39efp-1, 0x1.ebfbdff8131c4p-3,
167169
0x1.c6b08d7061695p-5, 0x1.3b2b1bee74b2ap-7,
168170
0x1.5d88091198529p-10};
169-
double dx_sq = dx * dx;
170-
double c1 = fputil::multiply_add(dx, COEFFS[0], 1.0);
171-
double c2 = fputil::multiply_add(dx, COEFFS[2], COEFFS[1]);
172-
double c3 = fputil::multiply_add(dx, COEFFS[4], COEFFS[3]);
173-
double p = fputil::multiply_add(dx_sq, c3, c2);
174-
// 2^x = 2^(hi + mid + lo)
175-
// = 2^(hi + mid) * 2^lo
176-
// ~ mh * (1 + lo * P(lo))
177-
// = mh + (mh*lo) * P(lo)
178-
double exp2_lo = fputil::multiply_add(p, dx_sq, c1);
171+
double lo_sq = lo * lo;
172+
double c1 = fputil::multiply_add(lo, COEFFS[0], 1.0);
173+
double c2 = fputil::multiply_add(lo, COEFFS[2], COEFFS[1]);
174+
double c3 = fputil::multiply_add(lo, COEFFS[4], COEFFS[3]);
175+
double exp2_lo = fputil::polyeval(lo_sq, c1, c2, c3);
176+
// 2^x - 1 = 2^(hi + mid + lo) - 1
177+
// = 2^(hi + mid) * 2^lo - 1
178+
// ~ mh * (1 + lo * P(lo)) - 1
179+
// = mh * exp2_lo - 1
179180
return static_cast<float>(fputil::multiply_add(exp2_lo, mh, -1.0));
180181
}
181182

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

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,6 @@ TEST_F(LlvmLibcExp2m1fTest, TrickyInputs) {
4141
LIBC_NAMESPACE::libc_errno = 0;
4242
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2m1, x,
4343
LIBC_NAMESPACE::exp2m1f(x), 0.5);
44-
EXPECT_MATH_ERRNO(0);
4544
}
4645
}
4746

‎libc/test/src/math/smoke/exp2m1f_test.cpp

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -19,15 +19,10 @@ TEST_F(LlvmLibcExp2m1fTest, SpecialNumbers) {
1919
LIBC_NAMESPACE::libc_errno = 0;
2020

2121
EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::exp2m1f(aNaN));
22-
EXPECT_MATH_ERRNO(0);
2322
EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::exp2m1f(inf));
24-
EXPECT_MATH_ERRNO(0);
2523
EXPECT_FP_EQ_ALL_ROUNDING(-1.0f, LIBC_NAMESPACE::exp2m1f(neg_inf));
26-
EXPECT_MATH_ERRNO(0);
2724
EXPECT_FP_EQ_ALL_ROUNDING(0.0f, LIBC_NAMESPACE::exp2m1f(0.0f));
28-
EXPECT_MATH_ERRNO(0);
2925
EXPECT_FP_EQ_ALL_ROUNDING(-0.0f, LIBC_NAMESPACE::exp2m1f(-0.0f));
30-
EXPECT_MATH_ERRNO(0);
3126

3227
EXPECT_FP_EQ_ALL_ROUNDING(1.0f, LIBC_NAMESPACE::exp2m1f(1.0f));
3328
EXPECT_FP_EQ_ALL_ROUNDING(-0.5f, LIBC_NAMESPACE::exp2m1f(-1.0f));

‎libc/utils/MPFRWrapper/MPFRUtils.cpp

Lines changed: 24 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -89,7 +89,8 @@ class MPFRNumber {
8989
// precision.
9090
template <typename XType,
9191
cpp::enable_if_t<cpp::is_same_v<float, XType>, int> = 0>
92-
explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE,
92+
explicit MPFRNumber(XType x,
93+
unsigned int precision = ExtraPrecision<XType>::VALUE,
9394
RoundingMode rounding = RoundingMode::Nearest)
9495
: mpfr_precision(precision),
9596
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
@@ -99,7 +100,8 @@ class MPFRNumber {
99100

100101
template <typename XType,
101102
cpp::enable_if_t<cpp::is_same_v<double, XType>, int> = 0>
102-
explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE,
103+
explicit MPFRNumber(XType x,
104+
unsigned int precision = ExtraPrecision<XType>::VALUE,
103105
RoundingMode rounding = RoundingMode::Nearest)
104106
: mpfr_precision(precision),
105107
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
@@ -109,7 +111,8 @@ class MPFRNumber {
109111

110112
template <typename XType,
111113
cpp::enable_if_t<cpp::is_same_v<long double, XType>, int> = 0>
112-
explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE,
114+
explicit MPFRNumber(XType x,
115+
unsigned int precision = ExtraPrecision<XType>::VALUE,
113116
RoundingMode rounding = RoundingMode::Nearest)
114117
: mpfr_precision(precision),
115118
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
@@ -119,7 +122,8 @@ class MPFRNumber {
119122

120123
template <typename XType,
121124
cpp::enable_if_t<cpp::is_integral_v<XType>, int> = 0>
122-
explicit MPFRNumber(XType x, int precision = ExtraPrecision<float>::VALUE,
125+
explicit MPFRNumber(XType x,
126+
unsigned int precision = ExtraPrecision<float>::VALUE,
123127
RoundingMode rounding = RoundingMode::Nearest)
124128
: mpfr_precision(precision),
125129
mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
@@ -244,17 +248,24 @@ class MPFRNumber {
244248
mpfr_exp2m1(result.value, value, mpfr_rounding);
245249
return result;
246250
#else
247-
MPFRNumber result(*this, mpfr_precision * 8);
248-
mpfr_exp2(result.value, value, mpfr_rounding);
249-
250-
if (mpfr_underflow_p()) {
251-
mpfr_clear_underflow();
252-
253-
if (mpfr_rounding == MPFR_RNDZ)
254-
return underflow_rndz_fallback;
251+
unsigned int prec = mpfr_precision * 3;
252+
MPFRNumber result(*this, prec);
253+
254+
float f = mpfr_get_flt(abs().value, mpfr_rounding);
255+
if (f > 0.5f && f < 0x1.0p30f) {
256+
mpfr_exp2(result.value, value, mpfr_rounding);
257+
mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding);
258+
return result;
255259
}
256260

257-
mpfr_sub_ui(result.value, result.value, 1, mpfr_rounding);
261+
MPFRNumber ln2(2.0f, prec);
262+
// log(2)
263+
mpfr_log(ln2.value, ln2.value, mpfr_rounding);
264+
// x * log(2)
265+
mpfr_mul(result.value, value, ln2.value, mpfr_rounding);
266+
// e^(x * log(2)) - 1
267+
int ex = mpfr_expm1(result.value, result.value, mpfr_rounding);
268+
mpfr_subnormalize(result.value, ex, mpfr_rounding);
258269
return result;
259270
#endif
260271
}

0 commit comments

Comments
 (0)