Skip to content

Commit 22d56a5

Browse files
overmightyovermighty
committed
[libc][math][c23] Add exp2m1f16 C23 math function
Part of #95250.
1 parent 547917a commit 22d56a5

File tree

13 files changed

+404
-42
lines changed

13 files changed

+404
-42
lines changed

‎libc/config/linux/x86_64/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -596,6 +596,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
596596
libc.src.math.copysignf16
597597
libc.src.math.exp10f16
598598
libc.src.math.exp2f16
599+
libc.src.math.exp2m1f16
599600
libc.src.math.expf16
600601
libc.src.math.expm1f16
601602
libc.src.math.f16add

‎libc/docs/math/index.rst

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -292,7 +292,7 @@ Higher Math Functions
292292
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
293293
| exp2 | |check| | |check| | | |check| | | 7.12.6.4 | F.10.3.4 |
294294
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
295-
| exp2m1 | |check| | | | | | 7.12.6.5 | F.10.3.5 |
295+
| exp2m1 | |check| | | | |check| | | 7.12.6.5 | F.10.3.5 |
296296
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
297297
| expm1 | |check| | |check| | | |check| | | 7.12.6.6 | F.10.3.6 |
298298
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

‎libc/spec/stdc.td

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -596,6 +596,7 @@ def StdC : StandardSpec<"stdc"> {
596596
GuardedFunctionSpec<"exp2f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
597597

598598
FunctionSpec<"exp2m1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
599+
GuardedFunctionSpec<"exp2m1f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
599600

600601
FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
601602
FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,

‎libc/src/math/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,7 @@ add_math_entrypoint_object(exp2f)
115115
add_math_entrypoint_object(exp2f16)
116116

117117
add_math_entrypoint_object(exp2m1f)
118+
add_math_entrypoint_object(exp2m1f16)
118119

119120
add_math_entrypoint_object(exp10)
120121
add_math_entrypoint_object(exp10f)

‎libc/src/math/exp2m1f16.h

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
//===-- Implementation header for exp2m1f16 ---------------------*- 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_MATH_EXP2M1F16_H
10+
#define LLVM_LIBC_SRC_MATH_EXP2M1F16_H
11+
12+
#include "src/__support/macros/config.h"
13+
#include "src/__support/macros/properties/types.h"
14+
15+
namespace LIBC_NAMESPACE_DECL {
16+
17+
float16 exp2m1f16(float16 x);
18+
19+
} // namespace LIBC_NAMESPACE_DECL
20+
21+
#endif // LLVM_LIBC_SRC_MATH_EXP2M1F16_H

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

Lines changed: 22 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1441,13 +1441,9 @@ add_entrypoint_object(
14411441
.expxf16
14421442
libc.hdr.errno_macros
14431443
libc.hdr.fenv_macros
1444-
libc.src.__support.CPP.array
14451444
libc.src.__support.FPUtil.except_value_utils
14461445
libc.src.__support.FPUtil.fenv_impl
14471446
libc.src.__support.FPUtil.fp_bits
1448-
libc.src.__support.FPUtil.multiply_add
1449-
libc.src.__support.FPUtil.nearest_integer
1450-
libc.src.__support.FPUtil.polyeval
14511447
libc.src.__support.FPUtil.rounding_mode
14521448
libc.src.__support.macros.optimization
14531449
COMPILE_OPTIONS
@@ -1475,6 +1471,28 @@ add_entrypoint_object(
14751471
-O3
14761472
)
14771473

1474+
add_entrypoint_object(
1475+
exp2m1f16
1476+
SRCS
1477+
exp2m1f16.cpp
1478+
HDRS
1479+
../exp2m1f16.h
1480+
DEPENDS
1481+
.expxf16
1482+
libc.hdr.errno_macros
1483+
libc.hdr.fenv_macros
1484+
libc.src.__support.common
1485+
libc.src.__support.FPUtil.except_value_utils
1486+
libc.src.__support.FPUtil.fenv_impl
1487+
libc.src.__support.FPUtil.fp_bits
1488+
libc.src.__support.FPUtil.multiply_add
1489+
libc.src.__support.FPUtil.polyeval
1490+
libc.src.__support.FPUtil.rounding_mode
1491+
libc.src.__support.macros.optimization
1492+
COMPILE_OPTIONS
1493+
-O3
1494+
)
1495+
14781496
add_entrypoint_object(
14791497
exp10
14801498
SRCS

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

Lines changed: 2 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -10,13 +10,9 @@
1010
#include "expxf16.h"
1111
#include "hdr/errno_macros.h"
1212
#include "hdr/fenv_macros.h"
13-
#include "src/__support/CPP/array.h"
1413
#include "src/__support/FPUtil/FEnvImpl.h"
1514
#include "src/__support/FPUtil/FPBits.h"
16-
#include "src/__support/FPUtil/PolyEval.h"
1715
#include "src/__support/FPUtil/except_value_utils.h"
18-
#include "src/__support/FPUtil/multiply_add.h"
19-
#include "src/__support/FPUtil/nearest_integer.h"
2016
#include "src/__support/FPUtil/rounding_mode.h"
2117
#include "src/__support/common.h"
2218
#include "src/__support/macros/config.h"
@@ -88,39 +84,8 @@ LLVM_LIBC_FUNCTION(float16, exp2f16, (float16 x)) {
8884
if (auto r = EXP2F16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
8985
return r.value();
9086

91-
// For -25 < x < 16, to compute 2^x, we perform the following range reduction:
92-
// find hi, mid, lo, such that:
93-
// x = hi + mid + lo, in which
94-
// hi is an integer,
95-
// mid * 2^3 is an integer,
96-
// -2^(-4) <= lo < 2^(-4).
97-
// In particular,
98-
// hi + mid = round(x * 2^3) * 2^(-3).
99-
// Then,
100-
// 2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
101-
// We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
102-
// by adding hi to the exponent field of 2^mid. 2^lo is computed using a
103-
// degree-3 minimax polynomial generated by Sollya.
104-
105-
float xf = x;
106-
float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
107-
int x_hi_mid = static_cast<int>(kf);
108-
int x_hi = x_hi_mid >> 3;
109-
int x_mid = x_hi_mid & 0x7;
110-
// lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x
111-
float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);
112-
113-
uint32_t exp2_hi_mid_bits =
114-
EXP2_MID_BITS[x_mid] +
115-
static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
116-
float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
117-
// Degree-3 minimax polynomial generated by Sollya with the following
118-
// commands:
119-
// > display = hexadecimal;
120-
// > P = fpminimax((2^x - 1)/x, 2, [|SG...|], [-2^-4, 2^-4]);
121-
// > 1 + x * P;
122-
float exp2_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.62e43p-1f, 0x1.ec0aa6p-3f,
123-
0x1.c6b4a6p-5f);
87+
// exp2(x) = exp2(hi + mid) * exp2(lo)
88+
auto [exp2_hi_mid, exp2_lo] = exp2_range_reduction(x);
12489
return static_cast<float16>(exp2_hi_mid * exp2_lo);
12590
}
12691

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

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
//===-- Half-precision 2^x - 1 function -----------------------------------===//
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+
#include "src/math/exp2m1f16.h"
10+
#include "expxf16.h"
11+
#include "hdr/errno_macros.h"
12+
#include "hdr/fenv_macros.h"
13+
#include "src/__support/FPUtil/FEnvImpl.h"
14+
#include "src/__support/FPUtil/FPBits.h"
15+
#include "src/__support/FPUtil/PolyEval.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+
23+
namespace LIBC_NAMESPACE_DECL {
24+
25+
static constexpr fputil::ExceptValues<float16, 6> EXP2M1F16_EXCEPTS_LO = {{
26+
// (input, RZ output, RU offset, RD offset, RN offset)
27+
// x = 0x1.cf4p-13, exp2m1f16(x) = 0x1.41p-13 (RZ)
28+
{0x0b3dU, 0x0904U, 1U, 0U, 1U},
29+
// x = 0x1.4fcp-12, exp2m1f16(x) = 0x1.d14p-13 (RZ)
30+
{0x0d3fU, 0x0b45U, 1U, 0U, 1U},
31+
// x = 0x1.63p-11, exp2m1f16(x) = 0x1.ec4p-12 (RZ)
32+
{0x118cU, 0x0fb1U, 1U, 0U, 0U},
33+
// x = 0x1.6fp-7, exp2m1f16(x) = 0x1.fe8p-8 (RZ)
34+
{0x21bcU, 0x1ffaU, 1U, 0U, 1U},
35+
// x = -0x1.c6p-10, exp2m1f16(x) = -0x1.3a8p-10 (RZ)
36+
{0x9718U, 0x94eaU, 0U, 1U, 0U},
37+
// x = -0x1.cfcp-10, exp2m1f16(x) = -0x1.414p-10 (RZ)
38+
{0x973fU, 0x9505U, 0U, 1U, 0U},
39+
}};
40+
41+
#ifdef LIBC_TARGET_CPU_HAS_FMA
42+
static constexpr size_t N_EXP2M1F16_EXCEPTS_HI = 6;
43+
#else
44+
static constexpr size_t N_EXP2M1F16_EXCEPTS_HI = 7;
45+
#endif
46+
47+
static constexpr fputil::ExceptValues<float16, N_EXP2M1F16_EXCEPTS_HI>
48+
EXP2M1F16_EXCEPTS_HI = {{
49+
// (input, RZ output, RU offset, RD offset, RN offset)
50+
// x = 0x1.e58p-3, exp2m1f16(x) = 0x1.6dcp-3 (RZ)
51+
{0x3396U, 0x31b7U, 1U, 0U, 0U},
52+
#ifndef LIBC_TARGET_CPU_HAS_FMA
53+
// x = 0x1.2e8p-2, exp2m1f16(x) = 0x1.d14p-3 (RZ)
54+
{0x34baU, 0x3345U, 1U, 0U, 0U},
55+
#endif
56+
// x = 0x1.ad8p-2, exp2m1f16(x) = 0x1.598p-2 (RZ)
57+
{0x36b6U, 0x3566U, 1U, 0U, 0U},
58+
#ifdef LIBC_TARGET_CPU_HAS_FMA
59+
// x = 0x1.edcp-2, exp2m1f16(x) = 0x1.964p-2 (RZ)
60+
{0x37b7U, 0x3659U, 1U, 0U, 1U},
61+
#endif
62+
// x = -0x1.804p-3, exp2m1f16(x) = -0x1.f34p-4 (RZ)
63+
{0xb201U, 0xafcdU, 0U, 1U, 1U},
64+
// x = -0x1.f3p-3, exp2m1f16(x) = -0x1.3e4p-3 (RZ)
65+
{0xb3ccU, 0xb0f9U, 0U, 1U, 0U},
66+
// x = -0x1.294p-1, exp2m1f16(x) = -0x1.53p-2 (RZ)
67+
{0xb8a5U, 0xb54cU, 0U, 1U, 1U},
68+
#ifndef LIBC_TARGET_CPU_HAS_FMA
69+
// x = -0x1.a34p-1, exp2m1f16(x) = -0x1.bb4p-2 (RZ)
70+
{0xba8dU, 0xb6edU, 0U, 1U, 1U},
71+
#endif
72+
}};
73+
74+
LLVM_LIBC_FUNCTION(float16, exp2m1f16, (float16 x)) {
75+
using FPBits = fputil::FPBits<float16>;
76+
FPBits x_bits(x);
77+
78+
uint16_t x_u = x_bits.uintval();
79+
uint16_t x_abs = x_u & 0x7fffU;
80+
81+
// When |x| <= 2^(-3), or |x| >= 11, or x is NaN.
82+
if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x4980U)) {
83+
// exp2m1(NaN) = NaN
84+
if (x_bits.is_nan()) {
85+
if (x_bits.is_signaling_nan()) {
86+
fputil::raise_except_if_required(FE_INVALID);
87+
return FPBits::quiet_nan().get_val();
88+
}
89+
90+
return x;
91+
}
92+
93+
// When x >= 16.
94+
if (x_u >= 0x4c00 && x_bits.is_pos()) {
95+
// exp2m1(+inf) = +inf
96+
if (x_bits.is_inf())
97+
return FPBits::inf().get_val();
98+
99+
switch (fputil::quick_get_round()) {
100+
case FE_TONEAREST:
101+
case FE_UPWARD:
102+
fputil::set_errno_if_required(ERANGE);
103+
fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
104+
return FPBits::inf().get_val();
105+
default:
106+
return FPBits::max_normal().get_val();
107+
}
108+
}
109+
110+
// When x < -11.
111+
if (x_u > 0xc980U) {
112+
// exp2m1(-inf) = -1
113+
if (x_bits.is_inf())
114+
return FPBits::one(Sign::NEG).get_val();
115+
116+
// When -12 < x < -11, round(2^x - 1, HP, RN) = -0x1.ffcp-1.
117+
if (x_u < 0xca00U) {
118+
return fputil::round_result_slightly_down(
119+
static_cast<float16>(-0x1.ffcp-1));
120+
}
121+
122+
// When x <= -12, round(2^x - 1, HP, RN) = -1.
123+
switch (fputil::quick_get_round()) {
124+
case FE_TONEAREST:
125+
case FE_DOWNWARD:
126+
return FPBits::one(Sign::NEG).get_val();
127+
default:
128+
return static_cast<float16>(-0x1.ffcp-1);
129+
}
130+
}
131+
132+
// When |x| <= 2^(-3).
133+
if (x_abs <= 0x3000U) {
134+
if (auto r = EXP2M1F16_EXCEPTS_LO.lookup(x_u);
135+
LIBC_UNLIKELY(r.has_value()))
136+
return r.value();
137+
138+
float xf = x;
139+
// Degree-5 minimax polynomial generated by Sollya with the following
140+
// commands:
141+
// > display = hexadecimal;
142+
// > P = fpminimax((2^x - 1)/x, 4, [|SG...|], [-2^-3, 2^-3]);
143+
// > x * P;
144+
return static_cast<float16>(
145+
xf * fputil::polyeval(xf, 0x1.62e43p-1f, 0x1.ebfbdep-3f,
146+
0x1.c6af88p-5f, 0x1.3b45d6p-7f,
147+
0x1.641e7cp-10f));
148+
}
149+
}
150+
151+
if (auto r = EXP2M1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
152+
return r.value();
153+
154+
// exp2(x) = exp2(hi + mid) * exp2(lo)
155+
auto [exp2_hi_mid, exp2_lo] = exp2_range_reduction(x);
156+
// exp2m1(x) = exp2(hi + mid) * exp2(lo) - 1
157+
return static_cast<float16>(
158+
fputil::multiply_add(exp2_hi_mid, exp2_lo, -1.0f));
159+
}
160+
161+
} // namespace LIBC_NAMESPACE_DECL

‎libc/src/math/generic/expxf16.h

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H
1111

1212
#include "src/__support/CPP/array.h"
13+
#include "src/__support/FPUtil/FPBits.h"
1314
#include "src/__support/FPUtil/PolyEval.h"
1415
#include "src/__support/FPUtil/multiply_add.h"
1516
#include "src/__support/FPUtil/nearest_integer.h"
@@ -89,6 +90,43 @@ constexpr cpp::array<uint32_t, 8> EXP2_MID_BITS = {
8990
0x3fb5'04f3U, 0x3fc5'672aU, 0x3fd7'44fdU, 0x3fea'c0c7U,
9091
};
9192

93+
LIBC_INLINE ExpRangeReduction exp2_range_reduction(float16 x) {
94+
// For -25 < x < 16, to compute 2^x, we perform the following range reduction:
95+
// find hi, mid, lo, such that:
96+
// x = hi + mid + lo, in which
97+
// hi is an integer,
98+
// mid * 2^3 is an integer,
99+
// -2^(-4) <= lo < 2^(-4).
100+
// In particular,
101+
// hi + mid = round(x * 2^3) * 2^(-3).
102+
// Then,
103+
// 2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
104+
// We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
105+
// by adding hi to the exponent field of 2^mid. 2^lo is computed using a
106+
// degree-3 minimax polynomial generated by Sollya.
107+
108+
float xf = x;
109+
float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
110+
int x_hi_mid = static_cast<int>(kf);
111+
int x_hi = x_hi_mid >> 3;
112+
int x_mid = x_hi_mid & 0x7;
113+
// lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x
114+
float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);
115+
116+
uint32_t exp2_hi_mid_bits =
117+
EXP2_MID_BITS[x_mid] +
118+
static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
119+
float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
120+
// Degree-3 minimax polynomial generated by Sollya with the following
121+
// commands:
122+
// > display = hexadecimal;
123+
// > P = fpminimax((2^x - 1)/x, 2, [|SG...|], [-2^-4, 2^-4]);
124+
// > 1 + x * P;
125+
float exp2_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.62e43p-1f, 0x1.ec0aa6p-3f,
126+
0x1.c6b4a6p-5f);
127+
return {exp2_hi_mid, exp2_lo};
128+
}
129+
92130
} // namespace LIBC_NAMESPACE_DECL
93131

94132
#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H

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

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1003,6 +1003,17 @@ add_fp_unittest(
10031003
libc.src.__support.FPUtil.fp_bits
10041004
)
10051005

1006+
add_fp_unittest(
1007+
exp2m1f16_test
1008+
NEED_MPFR
1009+
SUITE
1010+
libc-math-unittests
1011+
SRCS
1012+
exp2m1f16_test.cpp
1013+
DEPENDS
1014+
libc.src.math.exp2m1f16
1015+
)
1016+
10061017
add_fp_unittest(
10071018
exp10_test
10081019
NEED_MPFR

0 commit comments

Comments
 (0)