Skip to content

Commit ce65d4e

Browse files
overmightyovermighty
authored
[libc][math][c23] Add exp2m1f16 C23 math function (#105690)
Part of #95250.
1 parent 8f9cce0 commit ce65d4e

File tree

14 files changed

+410
-42
lines changed

14 files changed

+410
-42
lines changed

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -611,6 +611,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
611611
libc.src.math.copysignf16
612612
libc.src.math.exp10f16
613613
libc.src.math.exp2f16
614+
libc.src.math.exp2m1f16
614615
libc.src.math.expf16
615616
libc.src.math.expm1f16
616617
libc.src.math.f16add

‎libc/docs/math/index.rst

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -296,7 +296,7 @@ Higher Math Functions
296296
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
297297
| exp2 | |check| | |check| | | |check| | | 7.12.6.4 | F.10.3.4 |
298298
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
299-
| exp2m1 | |check| | | | | | 7.12.6.5 | F.10.3.5 |
299+
| exp2m1 | |check| | | | |check| | | 7.12.6.5 | F.10.3.5 |
300300
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
301301
| expm1 | |check| | |check| | | |check| | | 7.12.6.6 | F.10.3.6 |
302302
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

‎libc/newhdrgen/yaml/math.yaml

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,13 @@ functions:
284284
return_type: float
285285
arguments:
286286
- type: float
287+
- name: exp2m1f16
288+
standards:
289+
- stdc
290+
return_type: _Float16
291+
arguments:
292+
- type: _Float16
293+
guard: LIBC_TYPES_HAS_FLOAT16
287294
- name: expf
288295
standards:
289296
- stdc

‎libc/spec/stdc.td

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

684684
FunctionSpec<"exp2m1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
685+
GuardedFunctionSpec<"exp2m1f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
685686

686687
FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
687688
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
@@ -121,6 +121,7 @@ add_math_entrypoint_object(exp2f)
121121
add_math_entrypoint_object(exp2f16)
122122

123123
add_math_entrypoint_object(exp2m1f)
124+
add_math_entrypoint_object(exp2m1f16)
124125

125126
add_math_entrypoint_object(exp10)
126127
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: 26 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1500,14 +1500,10 @@ add_entrypoint_object(
15001500
.expxf16
15011501
libc.hdr.errno_macros
15021502
libc.hdr.fenv_macros
1503-
libc.src.__support.CPP.array
15041503
libc.src.__support.FPUtil.cast
15051504
libc.src.__support.FPUtil.except_value_utils
15061505
libc.src.__support.FPUtil.fenv_impl
15071506
libc.src.__support.FPUtil.fp_bits
1508-
libc.src.__support.FPUtil.multiply_add
1509-
libc.src.__support.FPUtil.nearest_integer
1510-
libc.src.__support.FPUtil.polyeval
15111507
libc.src.__support.FPUtil.rounding_mode
15121508
libc.src.__support.macros.optimization
15131509
COMPILE_OPTIONS
@@ -1535,6 +1531,30 @@ add_entrypoint_object(
15351531
-O3
15361532
)
15371533

1534+
add_entrypoint_object(
1535+
exp2m1f16
1536+
SRCS
1537+
exp2m1f16.cpp
1538+
HDRS
1539+
../exp2m1f16.h
1540+
DEPENDS
1541+
.expxf16
1542+
libc.hdr.errno_macros
1543+
libc.hdr.fenv_macros
1544+
libc.src.__support.common
1545+
libc.src.__support.FPUtil.cast
1546+
libc.src.__support.FPUtil.except_value_utils
1547+
libc.src.__support.FPUtil.fenv_impl
1548+
libc.src.__support.FPUtil.fp_bits
1549+
libc.src.__support.FPUtil.multiply_add
1550+
libc.src.__support.FPUtil.polyeval
1551+
libc.src.__support.FPUtil.rounding_mode
1552+
libc.src.__support.macros.optimization
1553+
libc.src.__support.macros.properties.cpu_features
1554+
COMPILE_OPTIONS
1555+
-O3
1556+
)
1557+
15381558
add_entrypoint_object(
15391559
exp10
15401560
SRCS
@@ -5235,7 +5255,9 @@ add_header_library(
52355255
expxf16.h
52365256
DEPENDS
52375257
libc.src.__support.CPP.array
5258+
libc.src.__support.FPUtil.fp_bits
52385259
libc.src.__support.FPUtil.multiply_add
52395260
libc.src.__support.FPUtil.nearest_integer
52405261
libc.src.__support.FPUtil.polyeval
5262+
libc.src.__support.macros.attributes
52415263
)

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

Lines changed: 2 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -10,14 +10,10 @@
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/cast.h"
1816
#include "src/__support/FPUtil/except_value_utils.h"
19-
#include "src/__support/FPUtil/multiply_add.h"
20-
#include "src/__support/FPUtil/nearest_integer.h"
2117
#include "src/__support/FPUtil/rounding_mode.h"
2218
#include "src/__support/common.h"
2319
#include "src/__support/macros/config.h"
@@ -89,39 +85,8 @@ LLVM_LIBC_FUNCTION(float16, exp2f16, (float16 x)) {
8985
if (auto r = EXP2F16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
9086
return r.value();
9187

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

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

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

0 commit comments

Comments
 (0)