Skip to content

Commit c25a5e0

Browse files
bassiounixbassiounix
authored
[libc][math] Refactor expf16 implementation to header-only in src/__support/math folder. (#147428)
Part of #147386 in preparation for: https://discourse.llvm.org/t/rfc-make-clang-builtin-math-functions-constexpr-with-llvm-libc-to-support-c-23-constexpr-math-functions/86450
1 parent 82acb59 commit c25a5e0

File tree

9 files changed

+336
-193
lines changed

9 files changed

+336
-193
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,5 +12,6 @@
1212
#include "libc_common.h"
1313

1414
#include "math/expf.h"
15+
#include "math/expf16.h"
1516

1617
#endif // LLVM_LIBC_SHARED_MATH_H

‎libc/shared/math/expf16.h

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
//===-- Shared expf16 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_EXPF16_H
10+
#define LLVM_LIBC_SHARED_MATH_EXPF16_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/expf16.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
namespace shared {
21+
22+
using math::expf16;
23+
24+
} // namespace shared
25+
} // namespace LIBC_NAMESPACE_DECL
26+
27+
#endif // LIBC_TYPES_HAS_FLOAT16
28+
29+
#endif // LLVM_LIBC_SHARED_MATH_EXPF16_H

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

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,3 +22,36 @@ add_header_library(
2222
libc.src.__support.macros.config
2323
libc.src.__support.macros.optimization
2424
)
25+
26+
add_header_library(
27+
expf16_utils
28+
HDRS
29+
expf16_utils.h
30+
DEPENDS
31+
libc.src.__support.CPP.array
32+
libc.src.__support.FPUtil.nearest_integer
33+
libc.src.__support.FPUtil.polyeval
34+
libc.src.__support.macros.attributes
35+
libc.include.llvm-libc-macros.float16_macros
36+
)
37+
38+
add_header_library(
39+
expf16
40+
HDRS
41+
expf16.h
42+
DEPENDS
43+
.expf16_utils
44+
libc.hdr.errno_macros
45+
libc.hdr.fenv_macros
46+
libc.src.__support.CPP.array
47+
libc.src.__support.FPUtil.cast
48+
libc.src.__support.FPUtil.except_value_utils
49+
libc.src.__support.FPUtil.fenv_impl
50+
libc.src.__support.FPUtil.fp_bits
51+
libc.src.__support.FPUtil.multiply_add
52+
libc.src.__support.FPUtil.nearest_integer
53+
libc.src.__support.FPUtil.polyeval
54+
libc.src.__support.FPUtil.rounding_mode
55+
libc.src.__support.macros.optimization
56+
libc.include.llvm-libc-macros.float16_macros
57+
)

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

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
//===-- Implementation header for expf16 ------------------------*- 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_EXPF16_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "hdr/errno_macros.h"
17+
#include "hdr/fenv_macros.h"
18+
#include "src/__support/FPUtil/FEnvImpl.h"
19+
#include "src/__support/FPUtil/FPBits.h"
20+
#include "src/__support/FPUtil/PolyEval.h"
21+
#include "src/__support/FPUtil/cast.h"
22+
#include "src/__support/FPUtil/except_value_utils.h"
23+
#include "src/__support/FPUtil/rounding_mode.h"
24+
#include "src/__support/common.h"
25+
#include "src/__support/macros/config.h"
26+
#include "src/__support/macros/optimization.h"
27+
28+
#include "expf16_utils.h"
29+
30+
namespace LIBC_NAMESPACE_DECL {
31+
32+
namespace math {
33+
34+
static constexpr float16 expf16(float16 x) {
35+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
36+
constexpr fputil::ExceptValues<float16, 2> EXPF16_EXCEPTS_LO = {{
37+
// (input, RZ output, RU offset, RD offset, RN offset)
38+
// x = 0x1.de4p-8, expf16(x) = 0x1.01cp+0 (RZ)
39+
{0x1f79U, 0x3c07U, 1U, 0U, 0U},
40+
// x = 0x1.73cp-6, expf16(x) = 0x1.05cp+0 (RZ)
41+
{0x25cfU, 0x3c17U, 1U, 0U, 0U},
42+
}};
43+
44+
constexpr fputil::ExceptValues<float16, 3> EXPF16_EXCEPTS_HI = {{
45+
// (input, RZ output, RU offset, RD offset, RN offset)
46+
// x = 0x1.c34p+0, expf16(x) = 0x1.74cp+2 (RZ)
47+
{0x3f0dU, 0x45d3U, 1U, 0U, 1U},
48+
// x = -0x1.488p-5, expf16(x) = 0x1.ebcp-1 (RZ)
49+
{0xa922U, 0x3bafU, 1U, 0U, 0U},
50+
// x = -0x1.55p-5, expf16(x) = 0x1.ebp-1 (RZ)
51+
{0xa954U, 0x3bacU, 1U, 0U, 0U},
52+
}};
53+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
54+
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 0 < |x| <= 2^(-5), or |x| >= 12, or x is NaN.
62+
if (LIBC_UNLIKELY(x_abs <= 0x2800U || x_abs >= 0x4a00U)) {
63+
// exp(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+
// When x >= 12.
74+
if (x_bits.is_pos() && x_abs >= 0x4a00U) {
75+
// exp(+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);
84+
return FPBits::inf().get_val();
85+
default:
86+
return FPBits::max_normal().get_val();
87+
}
88+
}
89+
90+
// When x <= -18.
91+
if (x_u >= 0xcc80U) {
92+
// exp(-inf) = +0
93+
if (x_bits.is_inf())
94+
return FPBits::zero().get_val();
95+
96+
fputil::set_errno_if_required(ERANGE);
97+
fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
98+
99+
switch (fputil::quick_get_round()) {
100+
case FE_UPWARD:
101+
return FPBits::min_subnormal().get_val();
102+
default:
103+
return FPBits::zero().get_val();
104+
}
105+
}
106+
107+
// When 0 < |x| <= 2^(-5).
108+
if (x_abs <= 0x2800U && !x_bits.is_zero()) {
109+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
110+
if (auto r = EXPF16_EXCEPTS_LO.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
111+
return r.value();
112+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
113+
114+
float xf = x;
115+
// Degree-3 minimax polynomial generated by Sollya with the following
116+
// commands:
117+
// > display = hexadecimal;
118+
// > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
119+
// > 1 + x * P;
120+
return fputil::cast<float16>(
121+
fputil::polyeval(xf, 0x1p+0f, 0x1p+0f, 0x1.0004p-1f, 0x1.555778p-3f));
122+
}
123+
}
124+
125+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
126+
if (auto r = EXPF16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
127+
return r.value();
128+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
129+
130+
// exp(x) = exp(hi + mid) * exp(lo)
131+
auto [exp_hi_mid, exp_lo] = exp_range_reduction(x);
132+
return fputil::cast<float16>(exp_hi_mid * exp_lo);
133+
}
134+
135+
} // namespace math
136+
137+
} // namespace LIBC_NAMESPACE_DECL
138+
139+
#endif // LIBC_TYPES_HAS_FLOAT16
140+
141+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_H

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

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
//===-- Common utils for expf16 functions -----------------------*- 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_EXPF16_UTILS_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "src/__support/CPP/array.h"
17+
#include "src/__support/FPUtil/PolyEval.h"
18+
#include "src/__support/FPUtil/nearest_integer.h"
19+
#include "src/__support/macros/properties/types.h"
20+
21+
namespace LIBC_NAMESPACE_DECL {
22+
23+
// Generated by Sollya with the following commands:
24+
// > display = hexadecimal;
25+
// > for i from -18 to 12 do print(round(exp(i), SG, RN));
26+
static constexpr cpp::array<float, 31> EXP_HI = {
27+
0x1.05a628p-26f, 0x1.639e32p-25f, 0x1.e355bcp-24f, 0x1.4875cap-22f,
28+
0x1.be6c7p-21f, 0x1.2f6054p-19f, 0x1.9c54c4p-18f, 0x1.183542p-16f,
29+
0x1.7cd79cp-15f, 0x1.02cf22p-13f, 0x1.5fc21p-12f, 0x1.de16bap-11f,
30+
0x1.44e52p-9f, 0x1.b993fep-8f, 0x1.2c155cp-6f, 0x1.97db0cp-5f,
31+
0x1.152aaap-3f, 0x1.78b564p-2f, 0x1p+0f, 0x1.5bf0a8p+1f,
32+
0x1.d8e64cp+2f, 0x1.415e5cp+4f, 0x1.b4c902p+5f, 0x1.28d38ap+7f,
33+
0x1.936dc6p+8f, 0x1.122886p+10f, 0x1.749ea8p+11f, 0x1.fa7158p+12f,
34+
0x1.5829dcp+14f, 0x1.d3c448p+15f, 0x1.3de166p+17f,
35+
};
36+
37+
// Generated by Sollya with the following commands:
38+
// > display = hexadecimal;
39+
// > for i from 0 to 7 do print(round(exp(i * 2^-3), SG, RN));
40+
static constexpr cpp::array<float, 8> EXP_MID = {
41+
0x1p+0f, 0x1.221604p+0f, 0x1.48b5e4p+0f, 0x1.747a52p+0f,
42+
0x1.a61298p+0f, 0x1.de455ep+0f, 0x1.0ef9dcp+1f, 0x1.330e58p+1f,
43+
};
44+
45+
struct ExpRangeReduction {
46+
float exp_hi_mid;
47+
float exp_lo;
48+
};
49+
50+
static constexpr ExpRangeReduction exp_range_reduction(float16 x) {
51+
// For -18 < x < 12, to compute exp(x), we perform the following range
52+
// reduction: find hi, mid, lo, such that:
53+
// x = hi + mid + lo, in which
54+
// hi is an integer,
55+
// mid * 2^3 is an integer,
56+
// -2^(-4) <= lo < 2^(-4).
57+
// In particular,
58+
// hi + mid = round(x * 2^3) * 2^(-3).
59+
// Then,
60+
// exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).
61+
// We store exp(hi) and exp(mid) in the lookup tables EXP_HI and EXP_MID
62+
// respectively. exp(lo) is computed using a degree-3 minimax polynomial
63+
// generated by Sollya.
64+
65+
float xf = x;
66+
float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
67+
int x_hi_mid = static_cast<int>(kf);
68+
int x_hi = x_hi_mid >> 3;
69+
int x_mid = x_hi_mid & 0x7;
70+
// lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x
71+
float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);
72+
73+
float exp_hi = EXP_HI[x_hi + 18];
74+
float exp_mid = EXP_MID[x_mid];
75+
// Degree-3 minimax polynomial generated by Sollya with the following
76+
// commands:
77+
// > display = hexadecimal;
78+
// > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-4, 2^-4]);
79+
// > 1 + x * P;
80+
float exp_lo =
81+
fputil::polyeval(lo, 0x1p+0f, 0x1p+0f, 0x1.001p-1f, 0x1.555ddep-3f);
82+
return {exp_hi * exp_mid, exp_lo};
83+
}
84+
85+
} // namespace LIBC_NAMESPACE_DECL
86+
87+
#endif // LIBC_TYPES_HAS_FLOAT16
88+
89+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H

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

Lines changed: 3 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1332,19 +1332,8 @@ add_entrypoint_object(
13321332
HDRS
13331333
../expf16.h
13341334
DEPENDS
1335-
.expxf16
1336-
libc.hdr.errno_macros
1337-
libc.hdr.fenv_macros
1338-
libc.src.__support.CPP.array
1339-
libc.src.__support.FPUtil.cast
1340-
libc.src.__support.FPUtil.except_value_utils
1341-
libc.src.__support.FPUtil.fenv_impl
1342-
libc.src.__support.FPUtil.fp_bits
1343-
libc.src.__support.FPUtil.multiply_add
1344-
libc.src.__support.FPUtil.nearest_integer
1345-
libc.src.__support.FPUtil.polyeval
1346-
libc.src.__support.FPUtil.rounding_mode
1347-
libc.src.__support.macros.optimization
1335+
libc.src.__support.math.expf16
1336+
libc.src.errno.errno
13481337
)
13491338

13501339
add_entrypoint_object(
@@ -5075,11 +5064,10 @@ add_header_library(
50755064
HDRS
50765065
expxf16.h
50775066
DEPENDS
5078-
libc.src.__support.CPP.array
50795067
libc.src.__support.FPUtil.cast
50805068
libc.src.__support.FPUtil.fp_bits
50815069
libc.src.__support.FPUtil.multiply_add
50825070
libc.src.__support.FPUtil.nearest_integer
5083-
libc.src.__support.FPUtil.polyeval
50845071
libc.src.__support.macros.attributes
5072+
libc.src.__support.math.expf16_utils
50855073
)

0 commit comments

Comments
 (0)