Skip to content

Commit 0f507ad

Browse files
bassiounixbassiounix
committed
[libc][math] Refactor expm1f16 implementation to header-only in src/__support/math folder.
1 parent 3d04ebf commit 0f507ad

File tree

9 files changed

+221
-140
lines changed

9 files changed

+221
-140
lines changed

‎libc/shared/math.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@
5656
#include "math/expf16.h"
5757
#include "math/expm1.h"
5858
#include "math/expm1f.h"
59+
#include "math/expm1f16.h"
5960
#include "math/frexpf.h"
6061
#include "math/frexpf128.h"
6162
#include "math/frexpf16.h"

‎libc/shared/math/expm1f16.h‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
//===-- Shared expm1f16 function --------------------------------*- C++ -*-===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#ifndef LLVM_LIBC_SHARED_MATH_EXPM1F16_H
10+
#define LLVM_LIBC_SHARED_MATH_EXPM1F16_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
#include "shared/libc_common.h"
14+
15+
#ifdef LIBC_TYPES_HAS_FLOAT16
16+
17+
#include "src/__support/math/expm1f16.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
namespace shared {
21+
22+
using math::expm1f16;
23+
24+
} // namespace shared
25+
} // namespace LIBC_NAMESPACE_DECL
26+
27+
#endif // LIBC_TYPES_HAS_FLOAT16
28+
29+
#endif // LLVM_LIBC_SHARED_MATH_EXPM1F16_H

‎libc/src/__support/math/CMakeLists.txt‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -907,6 +907,22 @@ add_header_library(
907907
libc.src.errno.errno
908908
)
909909

910+
add_header_library(
911+
expm1f16
912+
HDRS
913+
expm1f16.h
914+
DEPENDS
915+
.expxf16_utils
916+
libc.src.__support.FPUtil.cast
917+
libc.src.__support.FPUtil.except_value_utils
918+
libc.src.__support.FPUtil.fenv_impl
919+
libc.src.__support.FPUtil.fp_bits
920+
libc.src.__support.FPUtil.multiply_add
921+
libc.src.__support.FPUtil.polyeval
922+
libc.src.__support.FPUtil.rounding_mode
923+
libc.src.__support.macros.optimization
924+
)
925+
910926
add_header_library(
911927
range_reduction_double
912928
HDRS

‎libc/src/__support/math/expm1f16.h‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
//===-- Implementation header for expm1f16 ----------------------*- C++ -*-===//
2+
//
3+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4+
// See https://llvm.org/LICENSE.txt for license information.
5+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6+
//
7+
//===----------------------------------------------------------------------===//
8+
9+
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXPM1F16_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPM1F16_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "src/__support/FPUtil/FEnvImpl.h"
17+
#include "src/__support/FPUtil/FPBits.h"
18+
#include "src/__support/FPUtil/PolyEval.h"
19+
#include "src/__support/FPUtil/cast.h"
20+
#include "src/__support/FPUtil/except_value_utils.h"
21+
#include "src/__support/FPUtil/multiply_add.h"
22+
#include "src/__support/FPUtil/rounding_mode.h"
23+
#include "src/__support/common.h"
24+
#include "src/__support/macros/config.h"
25+
#include "src/__support/macros/optimization.h"
26+
#include "src/__support/math/expxf16_utils.h"
27+
28+
namespace LIBC_NAMESPACE_DECL {
29+
30+
namespace math {
31+
32+
LIBC_INLINE static constexpr float16 expm1f16(float16 x) {
33+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
34+
constexpr fputil::ExceptValues<float16, 1> EXPM1F16_EXCEPTS_LO = {{
35+
// (input, RZ output, RU offset, RD offset, RN offset)
36+
// x = 0x1.564p-5, expm1f16(x) = 0x1.5d4p-5 (RZ)
37+
{0x2959U, 0x2975U, 1U, 0U, 1U},
38+
}};
39+
40+
#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
41+
constexpr size_t N_EXPM1F16_EXCEPTS_HI = 2;
42+
#else
43+
constexpr size_t N_EXPM1F16_EXCEPTS_HI = 3;
44+
#endif
45+
46+
constexpr fputil::ExceptValues<float16, N_EXPM1F16_EXCEPTS_HI>
47+
EXPM1F16_EXCEPTS_HI = {{
48+
// (input, RZ output, RU offset, RD offset, RN offset)
49+
// x = 0x1.c34p+0, expm1f16(x) = 0x1.34cp+2 (RZ)
50+
{0x3f0dU, 0x44d3U, 1U, 0U, 1U},
51+
// x = -0x1.e28p-3, expm1f16(x) = -0x1.adcp-3 (RZ)
52+
{0xb38aU, 0xb2b7U, 0U, 1U, 1U},
53+
#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
54+
// x = 0x1.a08p-3, exp10m1f(x) = 0x1.cdcp-3 (RZ)
55+
{0x3282U, 0x3337U, 1U, 0U, 0U},
56+
#endif
57+
}};
58+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
59+
60+
using namespace math::expxf16_internal;
61+
using FPBits = fputil::FPBits<float16>;
62+
FPBits x_bits(x);
63+
64+
uint16_t x_u = x_bits.uintval();
65+
uint16_t x_abs = x_u & 0x7fffU;
66+
67+
// When |x| <= 2^(-3), or |x| >= -11 * log(2), or x is NaN.
68+
if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x47a0U)) {
69+
// expm1(NaN) = NaN
70+
if (x_bits.is_nan()) {
71+
if (x_bits.is_signaling_nan()) {
72+
fputil::raise_except_if_required(FE_INVALID);
73+
return FPBits::quiet_nan().get_val();
74+
}
75+
76+
return x;
77+
}
78+
79+
// expm1(+/-0) = +/-0
80+
if (x_abs == 0)
81+
return x;
82+
83+
// When x >= 16 * log(2).
84+
if (x_bits.is_pos() && x_abs >= 0x498cU) {
85+
// expm1(+inf) = +inf
86+
if (x_bits.is_inf())
87+
return FPBits::inf().get_val();
88+
89+
switch (fputil::quick_get_round()) {
90+
case FE_TONEAREST:
91+
case FE_UPWARD:
92+
fputil::set_errno_if_required(ERANGE);
93+
fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
94+
return FPBits::inf().get_val();
95+
default:
96+
return FPBits::max_normal().get_val();
97+
}
98+
}
99+
100+
// When x <= -11 * log(2).
101+
if (x_u >= 0xc7a0U) {
102+
// expm1(-inf) = -1
103+
if (x_bits.is_inf())
104+
return FPBits::one(Sign::NEG).get_val();
105+
106+
// When x > -0x1.0ap+3, round(expm1(x), HP, RN) = -1.
107+
if (x_u > 0xc828U)
108+
return fputil::round_result_slightly_up(
109+
FPBits::one(Sign::NEG).get_val());
110+
// When x <= -0x1.0ap+3, round(expm1(x), HP, RN) = -0x1.ffcp-1.
111+
return fputil::round_result_slightly_down(
112+
fputil::cast<float16>(-0x1.ffcp-1));
113+
}
114+
115+
// When 0 < |x| <= 2^(-3).
116+
if (x_abs <= 0x3000U && !x_bits.is_zero()) {
117+
118+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
119+
if (auto r = EXPM1F16_EXCEPTS_LO.lookup(x_u);
120+
LIBC_UNLIKELY(r.has_value()))
121+
return r.value();
122+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
123+
124+
float xf = x;
125+
// Degree-5 minimax polynomial generated by Sollya with the following
126+
// commands:
127+
// > display = hexadecimal;
128+
// > P = fpminimax(expm1(x)/x, 4, [|SG...|], [-2^-3, 2^-3]);
129+
// > x * P;
130+
return fputil::cast<float16>(
131+
xf * fputil::polyeval(xf, 0x1p+0f, 0x1.fffff8p-2f, 0x1.555556p-3f,
132+
0x1.55905ep-5f, 0x1.1124c2p-7f));
133+
}
134+
}
135+
136+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
137+
if (auto r = EXPM1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
138+
return r.value();
139+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
140+
141+
// exp(x) = exp(hi + mid) * exp(lo)
142+
auto [exp_hi_mid, exp_lo] = exp_range_reduction(x);
143+
// expm1(x) = exp(hi + mid) * exp(lo) - 1
144+
return fputil::cast<float16>(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0f));
145+
}
146+
147+
} // namespace math
148+
149+
} // namespace LIBC_NAMESPACE_DECL
150+
151+
#endif // LIBC_TYPES_HAS_FLOAT16
152+
153+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPM1F16_H

‎libc/src/math/generic/CMakeLists.txt‎

Lines changed: 1 addition & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1581,17 +1581,7 @@ add_entrypoint_object(
15811581
HDRS
15821582
../expm1f16.h
15831583
DEPENDS
1584-
libc.hdr.errno_macros
1585-
libc.hdr.fenv_macros
1586-
libc.src.__support.FPUtil.cast
1587-
libc.src.__support.FPUtil.except_value_utils
1588-
libc.src.__support.FPUtil.fenv_impl
1589-
libc.src.__support.FPUtil.fp_bits
1590-
libc.src.__support.FPUtil.multiply_add
1591-
libc.src.__support.FPUtil.polyeval
1592-
libc.src.__support.FPUtil.rounding_mode
1593-
libc.src.__support.macros.optimization
1594-
libc.src.__support.math.expxf16_utils
1584+
libc.src.__support.math.expm1f16
15951585
)
15961586

15971587
add_entrypoint_object(

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

Lines changed: 2 additions & 128 deletions
Original file line numberDiff line numberDiff line change
@@ -7,135 +7,9 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/expm1f16.h"
10-
#include "hdr/errno_macros.h"
11-
#include "hdr/fenv_macros.h"
12-
#include "src/__support/FPUtil/FEnvImpl.h"
13-
#include "src/__support/FPUtil/FPBits.h"
14-
#include "src/__support/FPUtil/PolyEval.h"
15-
#include "src/__support/FPUtil/cast.h"
16-
#include "src/__support/FPUtil/except_value_utils.h"
17-
#include "src/__support/FPUtil/multiply_add.h"
18-
#include "src/__support/FPUtil/rounding_mode.h"
19-
#include "src/__support/common.h"
20-
#include "src/__support/macros/config.h"
21-
#include "src/__support/macros/optimization.h"
22-
#include "src/__support/math/expxf16_utils.h"
10+
#include "src/__support/math/expm1f16.h"
2311

2412
namespace LIBC_NAMESPACE_DECL {
25-
26-
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
27-
static constexpr fputil::ExceptValues<float16, 1> EXPM1F16_EXCEPTS_LO = {{
28-
// (input, RZ output, RU offset, RD offset, RN offset)
29-
// x = 0x1.564p-5, expm1f16(x) = 0x1.5d4p-5 (RZ)
30-
{0x2959U, 0x2975U, 1U, 0U, 1U},
31-
}};
32-
33-
#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
34-
static constexpr size_t N_EXPM1F16_EXCEPTS_HI = 2;
35-
#else
36-
static constexpr size_t N_EXPM1F16_EXCEPTS_HI = 3;
37-
#endif
38-
39-
static constexpr fputil::ExceptValues<float16, N_EXPM1F16_EXCEPTS_HI>
40-
EXPM1F16_EXCEPTS_HI = {{
41-
// (input, RZ output, RU offset, RD offset, RN offset)
42-
// x = 0x1.c34p+0, expm1f16(x) = 0x1.34cp+2 (RZ)
43-
{0x3f0dU, 0x44d3U, 1U, 0U, 1U},
44-
// x = -0x1.e28p-3, expm1f16(x) = -0x1.adcp-3 (RZ)
45-
{0xb38aU, 0xb2b7U, 0U, 1U, 1U},
46-
#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
47-
// x = 0x1.a08p-3, exp10m1f(x) = 0x1.cdcp-3 (RZ)
48-
{0x3282U, 0x3337U, 1U, 0U, 0U},
49-
#endif
50-
}};
51-
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
52-
53-
LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) {
54-
using namespace math::expxf16_internal;
55-
using FPBits = fputil::FPBits<float16>;
56-
FPBits x_bits(x);
57-
58-
uint16_t x_u = x_bits.uintval();
59-
uint16_t x_abs = x_u & 0x7fffU;
60-
61-
// When |x| <= 2^(-3), or |x| >= -11 * log(2), or x is NaN.
62-
if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x47a0U)) {
63-
// expm1(NaN) = NaN
64-
if (x_bits.is_nan()) {
65-
if (x_bits.is_signaling_nan()) {
66-
fputil::raise_except_if_required(FE_INVALID);
67-
return FPBits::quiet_nan().get_val();
68-
}
69-
70-
return x;
71-
}
72-
73-
// expm1(+/-0) = +/-0
74-
if (x_abs == 0)
75-
return x;
76-
77-
// When x >= 16 * log(2).
78-
if (x_bits.is_pos() && x_abs >= 0x498cU) {
79-
// expm1(+inf) = +inf
80-
if (x_bits.is_inf())
81-
return FPBits::inf().get_val();
82-
83-
switch (fputil::quick_get_round()) {
84-
case FE_TONEAREST:
85-
case FE_UPWARD:
86-
fputil::set_errno_if_required(ERANGE);
87-
fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
88-
return FPBits::inf().get_val();
89-
default:
90-
return FPBits::max_normal().get_val();
91-
}
92-
}
93-
94-
// When x <= -11 * log(2).
95-
if (x_u >= 0xc7a0U) {
96-
// expm1(-inf) = -1
97-
if (x_bits.is_inf())
98-
return FPBits::one(Sign::NEG).get_val();
99-
100-
// When x > -0x1.0ap+3, round(expm1(x), HP, RN) = -1.
101-
if (x_u > 0xc828U)
102-
return fputil::round_result_slightly_up(
103-
FPBits::one(Sign::NEG).get_val());
104-
// When x <= -0x1.0ap+3, round(expm1(x), HP, RN) = -0x1.ffcp-1.
105-
return fputil::round_result_slightly_down(
106-
fputil::cast<float16>(-0x1.ffcp-1));
107-
}
108-
109-
// When 0 < |x| <= 2^(-3).
110-
if (x_abs <= 0x3000U && !x_bits.is_zero()) {
111-
112-
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
113-
if (auto r = EXPM1F16_EXCEPTS_LO.lookup(x_u);
114-
LIBC_UNLIKELY(r.has_value()))
115-
return r.value();
116-
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
117-
118-
float xf = x;
119-
// Degree-5 minimax polynomial generated by Sollya with the following
120-
// commands:
121-
// > display = hexadecimal;
122-
// > P = fpminimax(expm1(x)/x, 4, [|SG...|], [-2^-3, 2^-3]);
123-
// > x * P;
124-
return fputil::cast<float16>(
125-
xf * fputil::polyeval(xf, 0x1p+0f, 0x1.fffff8p-2f, 0x1.555556p-3f,
126-
0x1.55905ep-5f, 0x1.1124c2p-7f));
127-
}
128-
}
129-
130-
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
131-
if (auto r = EXPM1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
132-
return r.value();
133-
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
134-
135-
// exp(x) = exp(hi + mid) * exp(lo)
136-
auto [exp_hi_mid, exp_lo] = exp_range_reduction(x);
137-
// expm1(x) = exp(hi + mid) * exp(lo) - 1
138-
return fputil::cast<float16>(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0f));
139-
}
13+
LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) { return math::expm1f16(x); }
14014

14115
} // namespace LIBC_NAMESPACE_DECL

‎libc/test/shared/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@ add_fp_unittest(
4747
libc.src.__support.math.exp2m1f16
4848
libc.src.__support.math.expm1
4949
libc.src.__support.math.expm1f
50+
libc.src.__support.math.expm1f16
5051
libc.src.__support.math.exp10
5152
libc.src.__support.math.exp10f
5253
libc.src.__support.math.exp10f16

‎libc/test/shared/shared_math_test.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) {
3131
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp2f16(0.0f16));
3232
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::exp2m1f16(0.0f16));
3333
EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::expf16(0.0f16));
34+
EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::expm1f16(0.0f16));
3435

3536
ASSERT_FP_EQ(float16(8 << 5), LIBC_NAMESPACE::shared::ldexpf16(8.0f16, 5));
3637
ASSERT_FP_EQ(float16(-1 * (8 << 5)),

0 commit comments

Comments
 (0)