Skip to content

Commit 9d1279d

Browse files
bassiounixbassiounix
committed
[libc][math] Refactor exp10f16 implementation to header-only in src/__support/math folder.
1 parent 26646c5 commit 9d1279d

File tree

11 files changed

+357
-194
lines changed

11 files changed

+357
-194
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
#include "math/exp.h"
1515
#include "math/exp10.h"
1616
#include "math/exp10f.h"
17+
#include "math/exp10f16.h"
1718
#include "math/expf.h"
1819
#include "math/expf16.h"
1920
#include "math/frexpf.h"

‎libc/shared/math/exp10f16.h

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
//===-- Shared exp10f16 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_EXP10F_H
10+
#define LLVM_LIBC_SHARED_MATH_EXP10F_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "shared/libc_common.h"
17+
#include "src/__support/math/exp10f16.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
namespace shared {
21+
22+
using math::exp10f16;
23+
24+
} // namespace shared
25+
} // namespace LIBC_NAMESPACE_DECL
26+
27+
#endif // LIBC_TYPES_HAS_FLOAT16
28+
29+
#endif // LLVM_LIBC_SHARED_MATH_EXP10F_H

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

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -198,3 +198,37 @@ add_header_library(
198198
libc.src.__support.FPUtil.rounding_mode
199199
libc.src.__support.macros.optimization
200200
)
201+
202+
add_header_library(
203+
exp10_float16_constants
204+
HDRS
205+
exp10_float16_constants.h
206+
DEPENDS
207+
libc.src.__support.CPP.array
208+
)
209+
210+
add_header_library(
211+
exp10f16_utils
212+
HDRS
213+
exp10f16_utils.h
214+
DEPENDS
215+
.expf16_utils
216+
.exp10_float16_constants
217+
libc.src.__support.FPUtil.fp_bits
218+
)
219+
220+
add_header_library(
221+
exp10f16
222+
HDRS
223+
exp10f16.h
224+
DEPENDS
225+
.exp10f16_utils
226+
libc.src.__support.FPUtil.fp_bits
227+
src.__support.FPUtil.FEnvImpl
228+
src.__support.FPUtil.FPBits
229+
src.__support.FPUtil.cast
230+
src.__support.FPUtil.rounding_mode
231+
src.__support.FPUtil.except_value_utils
232+
src.__support.macros.optimization
233+
src.__support.macros.properties.cpu_features
234+
)

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

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
//===-- Constants for exp10f16 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_SRC___SUPPORT_MATH_EXP10_FLOAT16_CONSTANTS_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10_FLOAT16_CONSTANTS_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
#include <stdint.h>
14+
15+
#ifdef LIBC_TYPES_HAS_FLOAT16
16+
17+
#include "src/__support/CPP/array.h"
18+
19+
namespace LIBC_NAMESPACE_DECL {
20+
21+
// Generated by Sollya with the following commands:
22+
// > display = hexadecimal;
23+
// > for i from 0 to 7 do printsingle(round(2^(i * 2^-3), SG, RN));
24+
static constexpr cpp::array<uint32_t, 8> EXP2_MID_BITS = {
25+
0x3f80'0000U, 0x3f8b'95c2U, 0x3f98'37f0U, 0x3fa5'fed7U,
26+
0x3fb5'04f3U, 0x3fc5'672aU, 0x3fd7'44fdU, 0x3fea'c0c7U,
27+
};
28+
29+
// Generated by Sollya with the following commands:
30+
// > display = hexadecimal;
31+
// > round(log2(10), SG, RN);
32+
static constexpr float LOG2F_10 = 0x1.a934fp+1f;
33+
34+
// Generated by Sollya with the following commands:
35+
// > display = hexadecimal;
36+
// > round(log10(2), SG, RN);
37+
static constexpr float LOG10F_2 = 0x1.344136p-2f;
38+
39+
} // namespace LIBC_NAMESPACE_DECL
40+
41+
#endif // LIBC_TYPES_HAS_FLOAT16
42+
43+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F16_H

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

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
//===-- Implementation header for exp10f16 ----------------------*- 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_EXP10F16_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F16_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "exp10f16_utils.h"
17+
#include "src/__support/FPUtil/FEnvImpl.h"
18+
#include "src/__support/FPUtil/FPBits.h"
19+
#include "src/__support/FPUtil/cast.h"
20+
#include "src/__support/FPUtil/except_value_utils.h"
21+
#include "src/__support/FPUtil/rounding_mode.h"
22+
#include "src/__support/macros/config.h"
23+
#include "src/__support/macros/optimization.h"
24+
#include "src/__support/macros/properties/cpu_features.h"
25+
26+
namespace LIBC_NAMESPACE_DECL {
27+
28+
namespace math {
29+
30+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
31+
#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
32+
static constexpr size_t N_EXP10F16_EXCEPTS = 5;
33+
#else
34+
static constexpr size_t N_EXP10F16_EXCEPTS = 8;
35+
#endif
36+
37+
static constexpr fputil::ExceptValues<float16, N_EXP10F16_EXCEPTS>
38+
EXP10F16_EXCEPTS = {{
39+
// x = 0x1.8f4p-2, exp10f16(x) = 0x1.3ap+1 (RZ)
40+
{0x363dU, 0x40e8U, 1U, 0U, 1U},
41+
// x = 0x1.95cp-2, exp10f16(x) = 0x1.3ecp+1 (RZ)
42+
{0x3657U, 0x40fbU, 1U, 0U, 0U},
43+
// x = -0x1.018p-4, exp10f16(x) = 0x1.bbp-1 (RZ)
44+
{0xac06U, 0x3aecU, 1U, 0U, 0U},
45+
// x = -0x1.c28p+0, exp10f16(x) = 0x1.1ccp-6 (RZ)
46+
{0xbf0aU, 0x2473U, 1U, 0U, 0U},
47+
// x = -0x1.e1cp+1, exp10f16(x) = 0x1.694p-13 (RZ)
48+
{0xc387U, 0x09a5U, 1U, 0U, 0U},
49+
#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
50+
// x = 0x1.0cp+1, exp10f16(x) = 0x1.f04p+6 (RZ)
51+
{0x4030U, 0x57c1U, 1U, 0U, 1U},
52+
// x = 0x1.1b8p+1, exp10f16(x) = 0x1.47cp+7 (RZ)
53+
{0x406eU, 0x591fU, 1U, 0U, 1U},
54+
// x = 0x1.1b8p+2, exp10f16(x) = 0x1.a4p+14 (RZ)
55+
{0x446eU, 0x7690U, 1U, 0U, 1U},
56+
#endif
57+
}};
58+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
59+
60+
static constexpr float16 exp10f16(float16 x) {
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| >= 5, or x is NaN.
68+
if (LIBC_UNLIKELY(x_abs >= 0x4500U)) {
69+
// exp10(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+
// When x >= 5.
80+
if (x_bits.is_pos()) {
81+
// exp10(+inf) = +inf
82+
if (x_bits.is_inf())
83+
return FPBits::inf().get_val();
84+
85+
switch (fputil::quick_get_round()) {
86+
case FE_TONEAREST:
87+
case FE_UPWARD:
88+
fputil::set_errno_if_required(ERANGE);
89+
fputil::raise_except_if_required(FE_OVERFLOW);
90+
return FPBits::inf().get_val();
91+
default:
92+
return FPBits::max_normal().get_val();
93+
}
94+
}
95+
96+
// When x <= -8.
97+
if (x_u >= 0xc800U) {
98+
// exp10(-inf) = +0
99+
if (x_bits.is_inf())
100+
return FPBits::zero().get_val();
101+
102+
fputil::set_errno_if_required(ERANGE);
103+
fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
104+
105+
if (fputil::fenv_is_round_up())
106+
return FPBits::min_subnormal().get_val();
107+
return FPBits::zero().get_val();
108+
}
109+
}
110+
111+
// When x is 1, 2, 3, or 4. These are hard-to-round cases with exact results.
112+
if (LIBC_UNLIKELY((x_u & ~(0x3c00U | 0x4000U | 0x4200U | 0x4400U)) == 0)) {
113+
switch (x_u) {
114+
case 0x3c00U: // x = 1.0f16
115+
return fputil::cast<float16>(10.0);
116+
case 0x4000U: // x = 2.0f16
117+
return fputil::cast<float16>(100.0);
118+
case 0x4200U: // x = 3.0f16
119+
return fputil::cast<float16>(1'000.0);
120+
case 0x4400U: // x = 4.0f16
121+
return fputil::cast<float16>(10'000.0);
122+
}
123+
}
124+
125+
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
126+
if (auto r = EXP10F16_EXCEPTS.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
127+
return r.value();
128+
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
129+
130+
// 10^x = 2^((hi + mid) * log2(10)) * 10^lo
131+
auto [exp2_hi_mid, exp10_lo] = exp10_range_reduction(x);
132+
return fputil::cast<float16>(exp2_hi_mid * exp10_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_EXP10F16_H

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

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
//===-- Common utils for exp10f16 -------------------------------*- 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_EXP_FLOAT_CONSTANTS_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H
11+
12+
#include "include/llvm-libc-macros/float16-macros.h"
13+
14+
#ifdef LIBC_TYPES_HAS_FLOAT16
15+
16+
#include "exp10_float16_constants.h"
17+
#include "expf16_utils.h"
18+
#include "src/__support/FPUtil/FPBits.h"
19+
20+
namespace LIBC_NAMESPACE_DECL {
21+
22+
LIBC_INLINE static constexpr ExpRangeReduction
23+
exp10_range_reduction(float16 x) {
24+
// For -8 < x < 5, to compute 10^x, we perform the following range reduction:
25+
// find hi, mid, lo, such that:
26+
// x = (hi + mid) * log2(10) + lo, in which
27+
// hi is an integer,
28+
// mid * 2^3 is an integer,
29+
// -2^(-4) <= lo < 2^(-4).
30+
// In particular,
31+
// hi + mid = round(x * 2^3) * 2^(-3).
32+
// Then,
33+
// 10^x = 10^(hi + mid + lo) = 2^((hi + mid) * log2(10)) + 10^lo
34+
// We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
35+
// by adding hi to the exponent field of 2^mid. 10^lo is computed using a
36+
// degree-4 minimax polynomial generated by Sollya.
37+
38+
float xf = x;
39+
float kf = fputil::nearest_integer(xf * (LOG2F_10 * 0x1.0p+3f));
40+
int x_hi_mid = static_cast<int>(kf);
41+
unsigned x_hi = static_cast<unsigned>(x_hi_mid) >> 3;
42+
unsigned x_mid = static_cast<unsigned>(x_hi_mid) & 0x7;
43+
// lo = x - (hi + mid) = round(x * 2^3 * log2(10)) * log10(2) * (-2^(-3)) + x
44+
float lo = fputil::multiply_add(kf, LOG10F_2 * -0x1.0p-3f, xf);
45+
46+
uint32_t exp2_hi_mid_bits =
47+
EXP2_MID_BITS[x_mid] +
48+
static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
49+
float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
50+
// Degree-4 minimax polynomial generated by Sollya with the following
51+
// commands:
52+
// > display = hexadecimal;
53+
// > P = fpminimax((10^x - 1)/x, 3, [|SG...|], [-2^-4, 2^-4]);
54+
// > 1 + x * P;
55+
float exp10_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.26bb14p+1f, 0x1.53526p+1f,
56+
0x1.04b434p+1f, 0x1.2bcf9ep+0f);
57+
return {exp2_hi_mid, exp10_lo};
58+
}
59+
60+
} // namespace LIBC_NAMESPACE_DECL
61+
62+
#endif // LIBC_TYPES_HAS_FLOAT16
63+
64+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F16_H

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

Lines changed: 5 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1477,20 +1477,8 @@ add_entrypoint_object(
14771477
HDRS
14781478
../exp10f16.h
14791479
DEPENDS
1480-
.expxf16
1481-
libc.hdr.errno_macros
1482-
libc.hdr.fenv_macros
1483-
libc.src.__support.CPP.array
1484-
libc.src.__support.FPUtil.cast
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.nearest_integer
1490-
libc.src.__support.FPUtil.polyeval
1491-
libc.src.__support.FPUtil.rounding_mode
1492-
libc.src.__support.macros.optimization
1493-
libc.src.__support.macros.properties.cpu_features
1480+
libc.src.__support.math.exp10f16
1481+
libc.src.errno.errno
14941482
)
14951483

14961484
add_entrypoint_object(
@@ -1519,7 +1507,6 @@ add_entrypoint_object(
15191507
HDRS
15201508
../exp10m1f16.h
15211509
DEPENDS
1522-
.expxf16
15231510
libc.hdr.errno_macros
15241511
libc.hdr.fenv_macros
15251512
libc.src.__support.FPUtil.cast
@@ -1531,6 +1518,7 @@ add_entrypoint_object(
15311518
libc.src.__support.FPUtil.rounding_mode
15321519
libc.src.__support.macros.optimization
15331520
libc.src.__support.macros.properties.cpu_features
1521+
libc.src.__support.math.exp10f16_utils
15341522
)
15351523

15361524
add_entrypoint_object(
@@ -5025,10 +5013,11 @@ add_header_library(
50255013
HDRS
50265014
expxf16.h
50275015
DEPENDS
5028-
libc.src.__support.FPUtil.cast
50295016
libc.src.__support.FPUtil.fp_bits
5017+
libc.src.__support.FPUtil.cast
50305018
libc.src.__support.FPUtil.multiply_add
50315019
libc.src.__support.FPUtil.nearest_integer
50325020
libc.src.__support.macros.attributes
50335021
libc.src.__support.math.expf16_utils
5022+
libc.src.__support.math.exp10_float16_constants
50345023
)

0 commit comments

Comments
 (0)