Skip to content

Commit 560ed8c

Browse files
overmightyovermighty
authored
[libc][math][c23] Add expm1f16 C23 math function (#102387)
Part of #95250.
1 parent f364b2e commit 560ed8c

File tree

14 files changed

+469
-75
lines changed

14 files changed

+469
-75
lines changed

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -597,6 +597,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
597597
libc.src.math.exp10f16
598598
libc.src.math.exp2f16
599599
libc.src.math.expf16
600+
libc.src.math.expm1f16
600601
libc.src.math.f16add
601602
libc.src.math.f16addf
602603
libc.src.math.f16addl

‎libc/docs/math/index.rst

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -294,7 +294,7 @@ Higher Math Functions
294294
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
295295
| exp2m1 | |check| | | | | | 7.12.6.5 | F.10.3.5 |
296296
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
297-
| expm1 | |check| | |check| | | | | 7.12.6.6 | F.10.3.6 |
297+
| expm1 | |check| | |check| | | |check| | | 7.12.6.6 | F.10.3.6 |
298298
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
299299
| fma | |check| | |check| | | | | 7.12.13.1 | F.10.10.1 |
300300
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+

‎libc/spec/stdc.td

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -599,6 +599,7 @@ def StdC : StandardSpec<"stdc"> {
599599

600600
FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
601601
FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
602+
GuardedFunctionSpec<"expm1f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
602603

603604
FunctionSpec<"exp10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
604605
FunctionSpec<"exp10f", 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(exp10f16)
121121

122122
add_math_entrypoint_object(expm1)
123123
add_math_entrypoint_object(expm1f)
124+
add_math_entrypoint_object(expm1f16)
124125

125126
add_math_entrypoint_object(f16add)
126127
add_math_entrypoint_object(f16addf)

‎libc/src/math/expm1f16.h

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
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_MATH_EXPM1F16_H
10+
#define LLVM_LIBC_SRC_MATH_EXPM1F16_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 expm1f16(float16 x);
18+
19+
} // namespace LIBC_NAMESPACE_DECL
20+
21+
#endif // LLVM_LIBC_SRC_MATH_EXPM1F16_H

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

Lines changed: 25 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1359,14 +1359,12 @@ add_entrypoint_object(
13591359
HDRS
13601360
../expf16.h
13611361
DEPENDS
1362+
.expxf16
13621363
libc.hdr.errno_macros
13631364
libc.hdr.fenv_macros
1364-
libc.src.__support.CPP.array
13651365
libc.src.__support.FPUtil.except_value_utils
13661366
libc.src.__support.FPUtil.fenv_impl
13671367
libc.src.__support.FPUtil.fp_bits
1368-
libc.src.__support.FPUtil.multiply_add
1369-
libc.src.__support.FPUtil.nearest_integer
13701368
libc.src.__support.FPUtil.polyeval
13711369
libc.src.__support.FPUtil.rounding_mode
13721370
libc.src.__support.macros.optimization
@@ -1608,6 +1606,27 @@ add_entrypoint_object(
16081606
-O3
16091607
)
16101608

1609+
add_entrypoint_object(
1610+
expm1f16
1611+
SRCS
1612+
expm1f16.cpp
1613+
HDRS
1614+
../expm1f16.h
1615+
DEPENDS
1616+
.expxf16
1617+
libc.hdr.errno_macros
1618+
libc.hdr.fenv_macros
1619+
libc.src.__support.FPUtil.except_value_utils
1620+
libc.src.__support.FPUtil.fenv_impl
1621+
libc.src.__support.FPUtil.fp_bits
1622+
libc.src.__support.FPUtil.multiply_add
1623+
libc.src.__support.FPUtil.polyeval
1624+
libc.src.__support.FPUtil.rounding_mode
1625+
libc.src.__support.macros.optimization
1626+
COMPILE_OPTIONS
1627+
-O3
1628+
)
1629+
16111630
add_entrypoint_object(
16121631
powf
16131632
SRCS
@@ -5092,4 +5111,7 @@ add_header_library(
50925111
expxf16.h
50935112
DEPENDS
50945113
libc.src.__support.CPP.array
5114+
libc.src.__support.FPUtil.multiply_add
5115+
libc.src.__support.FPUtil.nearest_integer
5116+
libc.src.__support.FPUtil.polyeval
50955117
)

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

Lines changed: 4 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -7,15 +7,13 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/expf16.h"
10+
#include "expxf16.h"
1011
#include "hdr/errno_macros.h"
1112
#include "hdr/fenv_macros.h"
12-
#include "src/__support/CPP/array.h"
1313
#include "src/__support/FPUtil/FEnvImpl.h"
1414
#include "src/__support/FPUtil/FPBits.h"
1515
#include "src/__support/FPUtil/PolyEval.h"
1616
#include "src/__support/FPUtil/except_value_utils.h"
17-
#include "src/__support/FPUtil/multiply_add.h"
18-
#include "src/__support/FPUtil/nearest_integer.h"
1917
#include "src/__support/FPUtil/rounding_mode.h"
2018
#include "src/__support/common.h"
2119
#include "src/__support/macros/config.h"
@@ -41,28 +39,6 @@ static constexpr fputil::ExceptValues<float16, 3> EXPF16_EXCEPTS_HI = {{
4139
{0xa954U, 0x3bacU, 1U, 0U, 0U},
4240
}};
4341

44-
// Generated by Sollya with the following commands:
45-
// > display = hexadecimal;
46-
// > for i from -18 to 12 do print(round(exp(i), SG, RN));
47-
static constexpr cpp::array<float, 31> EXP_HI = {
48-
0x1.05a628p-26f, 0x1.639e32p-25f, 0x1.e355bcp-24f, 0x1.4875cap-22f,
49-
0x1.be6c7p-21f, 0x1.2f6054p-19f, 0x1.9c54c4p-18f, 0x1.183542p-16f,
50-
0x1.7cd79cp-15f, 0x1.02cf22p-13f, 0x1.5fc21p-12f, 0x1.de16bap-11f,
51-
0x1.44e52p-9f, 0x1.b993fep-8f, 0x1.2c155cp-6f, 0x1.97db0cp-5f,
52-
0x1.152aaap-3f, 0x1.78b564p-2f, 0x1p+0f, 0x1.5bf0a8p+1f,
53-
0x1.d8e64cp+2f, 0x1.415e5cp+4f, 0x1.b4c902p+5f, 0x1.28d38ap+7f,
54-
0x1.936dc6p+8f, 0x1.122886p+10f, 0x1.749ea8p+11f, 0x1.fa7158p+12f,
55-
0x1.5829dcp+14f, 0x1.d3c448p+15f, 0x1.3de166p+17f,
56-
};
57-
58-
// Generated by Sollya with the following commands:
59-
// > display = hexadecimal;
60-
// > for i from 0 to 7 do print(round(exp(i * 2^-3), SG, RN));
61-
static constexpr cpp::array<float, 8> EXP_MID = {
62-
0x1p+0f, 0x1.221604p+0f, 0x1.48b5e4p+0f, 0x1.747a52p+0f,
63-
0x1.a61298p+0f, 0x1.de455ep+0f, 0x1.0ef9dcp+1f, 0x1.330e58p+1f,
64-
};
65-
6642
LLVM_LIBC_FUNCTION(float16, expf16, (float16 x)) {
6743
using FPBits = fputil::FPBits<float16>;
6844
FPBits x_bits(x);
@@ -135,38 +111,9 @@ LLVM_LIBC_FUNCTION(float16, expf16, (float16 x)) {
135111
if (auto r = EXPF16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
136112
return r.value();
137113

138-
// For -18 < x < 12, to compute exp(x), we perform the following range
139-
// reduction: find hi, mid, lo, such that:
140-
// x = hi + mid + lo, in which
141-
// hi is an integer,
142-
// mid * 2^3 is an integer,
143-
// -2^(-4) <= lo < 2^(-4).
144-
// In particular,
145-
// hi + mid = round(x * 2^3) * 2^(-3).
146-
// Then,
147-
// exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).
148-
// We store exp(hi) and exp(mid) in the lookup tables EXP_HI and EXP_MID
149-
// respectively. exp(lo) is computed using a degree-3 minimax polynomial
150-
// generated by Sollya.
151-
152-
float xf = x;
153-
float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
154-
int x_hi_mid = static_cast<int>(kf);
155-
int x_hi = x_hi_mid >> 3;
156-
int x_mid = x_hi_mid & 0x7;
157-
// lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x
158-
float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);
159-
160-
float exp_hi = EXP_HI[x_hi + 18];
161-
float exp_mid = EXP_MID[x_mid];
162-
// Degree-3 minimax polynomial generated by Sollya with the following
163-
// commands:
164-
// > display = hexadecimal;
165-
// > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-4, 2^-4]);
166-
// > 1 + x * P;
167-
float exp_lo =
168-
fputil::polyeval(lo, 0x1p+0f, 0x1p+0f, 0x1.001p-1f, 0x1.555ddep-3f);
169-
return static_cast<float16>(exp_hi * exp_mid * exp_lo);
114+
// exp(x) = exp(hi + mid) * exp(lo)
115+
auto [exp_hi_mid, exp_lo] = exp_range_reduction(x);
116+
return static_cast<float16>(exp_hi_mid * exp_lo);
170117
}
171118

172119
} // namespace LIBC_NAMESPACE_DECL

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

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

0 commit comments

Comments
 (0)