Skip to content

Commit 392da0b

Browse files
bassiounixbassiounix
committed
[libc][math] Refactor exp10f implementation to header-only in src/__support/math folder.
1 parent dc6914c commit 392da0b

File tree

13 files changed

+281
-220
lines changed

13 files changed

+281
-220
lines changed

‎libc/shared/math.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313

1414
#include "math/exp.h"
1515
#include "math/exp10.h"
16+
#include "math/exp10f.h"
1617
#include "math/expf.h"
1718
#include "math/expf16.h"
1819
#include "math/frexpf.h"

‎libc/shared/math/exp10f.h

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
//===-- Shared exp10f function -----------------------------------*- C++
2+
//-*-===//
3+
//
4+
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
5+
// See https://llvm.org/LICENSE.txt for license information.
6+
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
7+
//
8+
//===----------------------------------------------------------------------===//
9+
10+
#ifndef LLVM_LIBC_SHARED_MATH_EXP10F_H
11+
#define LLVM_LIBC_SHARED_MATH_EXP10F_H
12+
13+
#include "shared/libc_common.h"
14+
#include "src/__support/math/exp10f.h"
15+
16+
namespace LIBC_NAMESPACE_DECL {
17+
namespace shared {
18+
19+
using math::exp10f;
20+
21+
} // namespace shared
22+
} // namespace LIBC_NAMESPACE_DECL
23+
24+
#endif // LLVM_LIBC_SHARED_MATH_EXP10F_H

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

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -170,3 +170,31 @@ add_header_library(
170170
libc.src.__support.integer_literals
171171
libc.src.__support.macros.optimization
172172
)
173+
174+
add_header_library(
175+
exp10f_utils
176+
HDRS
177+
exp10f_utils.h
178+
DEPENDS
179+
libc.src.__support.FPUtil.basic_operations
180+
libc.src.__support.FPUtil.fenv_impl
181+
libc.src.__support.FPUtil.multiply_add
182+
libc.src.__support.FPUtil.nearest_integer
183+
libc.src.__support.FPUtil.polyeval
184+
libc.src.__support.common
185+
libc.src.__support.math.exp_utils
186+
)
187+
188+
add_header_library(
189+
exp10f
190+
HDRS
191+
exp10f.h
192+
DEPENDS
193+
.exp10f_utils
194+
libc.src.__support.macros.config
195+
libc.src.__support.FPUtil.fenv_impl
196+
libc.src.__support.FPUtil.fp_bits
197+
libc.src.__support.FPUtil.multiply_add
198+
libc.src.__support.FPUtil.rounding_mode
199+
libc.src.__support.macros.optimization
200+
)

‎libc/src/math/generic/exp10f_impl.h renamed to ‎libc/src/__support/math/exp10f.h

Lines changed: 8 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,27 +1,26 @@
1-
//===-- Single-precision 10^x function ------------------------------------===//
1+
//===-- Implementation header for exp10f ------------------------*- C++ -*-===//
22
//
33
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
44
// See https://llvm.org/LICENSE.txt for license information.
55
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
66
//
77
//===----------------------------------------------------------------------===//
88

9-
#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H
10-
#define LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H
9+
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H
10+
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H
1111

12-
#include "explogxf.h"
12+
#include "exp10f_utils.h"
1313
#include "src/__support/FPUtil/FEnvImpl.h"
1414
#include "src/__support/FPUtil/FPBits.h"
1515
#include "src/__support/FPUtil/multiply_add.h"
1616
#include "src/__support/FPUtil/rounding_mode.h"
17-
#include "src/__support/common.h"
1817
#include "src/__support/macros/config.h"
1918
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
2019

2120
namespace LIBC_NAMESPACE_DECL {
22-
namespace generic {
21+
namespace math {
2322

24-
LIBC_INLINE float exp10f(float x) {
23+
static constexpr float exp10f(float x) {
2524
using FPBits = typename fputil::FPBits<float>;
2625
FPBits xbits(x);
2726

@@ -132,7 +131,7 @@ LIBC_INLINE float exp10f(float x) {
132131
return static_cast<float>(multiply_add(p, lo2 * rr.mh, c0 * rr.mh));
133132
}
134133

135-
} // namespace generic
134+
} // namespace math
136135
} // namespace LIBC_NAMESPACE_DECL
137136

138-
#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXP10F_IMPL_H
137+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP10F_H

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

Lines changed: 171 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
//===-- Common utils for exp10f ---------------------------------*- 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 "src/__support/FPUtil/FPBits.h"
13+
#include "src/__support/FPUtil/PolyEval.h"
14+
#include "src/__support/FPUtil/nearest_integer.h"
15+
#include "src/__support/macros/config.h"
16+
17+
namespace LIBC_NAMESPACE_DECL {
18+
19+
struct ExpBase {
20+
// Base = e
21+
static constexpr int MID_BITS = 5;
22+
static constexpr int MID_MASK = (1 << MID_BITS) - 1;
23+
// log2(e) * 2^5
24+
static constexpr double LOG2_B = 0x1.71547652b82fep+0 * (1 << MID_BITS);
25+
// High and low parts of -log(2) * 2^(-5)
26+
static constexpr double M_LOGB_2_HI = -0x1.62e42fefa0000p-1 / (1 << MID_BITS);
27+
static constexpr double M_LOGB_2_LO =
28+
-0x1.cf79abc9e3b3ap-40 / (1 << MID_BITS);
29+
// Look up table for bit fields of 2^(i/32) for i = 0..31, generated by Sollya
30+
// with:
31+
// > for i from 0 to 31 do printdouble(round(2^(i/32), D, RN));
32+
static constexpr int64_t EXP_2_MID[1 << MID_BITS] = {
33+
0x3ff0000000000000, 0x3ff059b0d3158574, 0x3ff0b5586cf9890f,
34+
0x3ff11301d0125b51, 0x3ff172b83c7d517b, 0x3ff1d4873168b9aa,
35+
0x3ff2387a6e756238, 0x3ff29e9df51fdee1, 0x3ff306fe0a31b715,
36+
0x3ff371a7373aa9cb, 0x3ff3dea64c123422, 0x3ff44e086061892d,
37+
0x3ff4bfdad5362a27, 0x3ff5342b569d4f82, 0x3ff5ab07dd485429,
38+
0x3ff6247eb03a5585, 0x3ff6a09e667f3bcd, 0x3ff71f75e8ec5f74,
39+
0x3ff7a11473eb0187, 0x3ff82589994cce13, 0x3ff8ace5422aa0db,
40+
0x3ff93737b0cdc5e5, 0x3ff9c49182a3f090, 0x3ffa5503b23e255d,
41+
0x3ffae89f995ad3ad, 0x3ffb7f76f2fb5e47, 0x3ffc199bdd85529c,
42+
0x3ffcb720dcef9069, 0x3ffd5818dcfba487, 0x3ffdfc97337b9b5f,
43+
0x3ffea4afa2a490da, 0x3fff50765b6e4540,
44+
};
45+
46+
// Approximating e^dx with degree-5 minimax polynomial generated by Sollya:
47+
// > Q = fpminimax(expm1(x)/x, 4, [|1, D...|], [-log(2)/64, log(2)/64]);
48+
// Then:
49+
// e^dx ~ P(dx) = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[3] * dx^5.
50+
static constexpr double COEFFS[4] = {
51+
0x1.ffffffffe5bc8p-2, 0x1.555555555cd67p-3, 0x1.5555c2a9b48b4p-5,
52+
0x1.11112a0e34bdbp-7};
53+
54+
LIBC_INLINE static double powb_lo(double dx) {
55+
using fputil::multiply_add;
56+
double dx2 = dx * dx;
57+
double c0 = 1.0 + dx;
58+
// c1 = COEFFS[0] + COEFFS[1] * dx
59+
double c1 = multiply_add(dx, ExpBase::COEFFS[1], ExpBase::COEFFS[0]);
60+
// c2 = COEFFS[2] + COEFFS[3] * dx
61+
double c2 = multiply_add(dx, ExpBase::COEFFS[3], ExpBase::COEFFS[2]);
62+
// r = c4 + c5 * dx^4
63+
// = 1 + dx + COEFFS[0] * dx^2 + ... + COEFFS[5] * dx^7
64+
return fputil::polyeval(dx2, c0, c1, c2);
65+
}
66+
};
67+
68+
struct Exp10Base : public ExpBase {
69+
// log2(10) * 2^5
70+
static constexpr double LOG2_B = 0x1.a934f0979a371p1 * (1 << MID_BITS);
71+
// High and low parts of -log10(2) * 2^(-5).
72+
// Notice that since |x * log2(10)| < 150:
73+
// |k| = |round(x * log2(10) * 2^5)| < 2^8 * 2^5 = 2^13
74+
// So when the FMA instructions are not available, in order for the product
75+
// k * M_LOGB_2_HI
76+
// to be exact, we only store the high part of log10(2) up to 38 bits
77+
// (= 53 - 15) of precision.
78+
// It is generated by Sollya with:
79+
// > round(log10(2), 44, RN);
80+
static constexpr double M_LOGB_2_HI = -0x1.34413509f8p-2 / (1 << MID_BITS);
81+
// > round(log10(2) - 0x1.34413509f8p-2, D, RN);
82+
static constexpr double M_LOGB_2_LO = 0x1.80433b83b532ap-44 / (1 << MID_BITS);
83+
84+
// Approximating 10^dx with degree-5 minimax polynomial generated by Sollya:
85+
// > Q = fpminimax((10^x - 1)/x, 4, [|D...|], [-log10(2)/2^6, log10(2)/2^6]);
86+
// Then:
87+
// 10^dx ~ P(dx) = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
88+
static constexpr double COEFFS[5] = {0x1.26bb1bbb55515p1, 0x1.53524c73bd3eap1,
89+
0x1.0470591dff149p1, 0x1.2bd7c0a9fbc4dp0,
90+
0x1.1429e74a98f43p-1};
91+
92+
static double powb_lo(double dx) {
93+
using fputil::multiply_add;
94+
double dx2 = dx * dx;
95+
// c0 = 1 + COEFFS[0] * dx
96+
double c0 = multiply_add(dx, Exp10Base::COEFFS[0], 1.0);
97+
// c1 = COEFFS[1] + COEFFS[2] * dx
98+
double c1 = multiply_add(dx, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]);
99+
// c2 = COEFFS[3] + COEFFS[4] * dx
100+
double c2 = multiply_add(dx, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]);
101+
// r = c0 + dx^2 * (c1 + c2 * dx^2)
102+
// = c0 + c1 * dx^2 + c2 * dx^4
103+
// = 1 + COEFFS[0] * dx + ... + COEFFS[4] * dx^5.
104+
return fputil::polyeval(dx2, c0, c1, c2);
105+
}
106+
};
107+
108+
constexpr int LOG_P1_BITS = 6;
109+
constexpr int LOG_P1_SIZE = 1 << LOG_P1_BITS;
110+
111+
// N[Table[Log[2, 1 + x], {x, 0/64, 63/64, 1/64}], 40]
112+
extern const double LOG_P1_LOG2[LOG_P1_SIZE];
113+
114+
// N[Table[1/(1 + x), {x, 0/64, 63/64, 1/64}], 40]
115+
extern const double LOG_P1_1_OVER[LOG_P1_SIZE];
116+
117+
// Taylor series expansion for Log[2, 1 + x] splitted to EVEN AND ODD numbers
118+
// K_LOG2_ODD starts from x^3
119+
extern const double K_LOG2_ODD[4];
120+
extern const double K_LOG2_EVEN[4];
121+
122+
// Output of range reduction for exp_b: (2^(mid + hi), lo)
123+
// where:
124+
// b^x = 2^(mid + hi) * b^lo
125+
struct exp_b_reduc_t {
126+
double mh; // 2^(mid + hi)
127+
double lo;
128+
};
129+
130+
// The function correctly calculates b^x value with at least float precision
131+
// in a limited range.
132+
// Range reduction:
133+
// b^x = 2^(hi + mid) * b^lo
134+
// where:
135+
// x = (hi + mid) * log_b(2) + lo
136+
// hi is an integer,
137+
// 0 <= mid * 2^MID_BITS < 2^MID_BITS is an integer
138+
// -2^(-MID_BITS - 1) <= lo * log2(b) <= 2^(-MID_BITS - 1)
139+
// Base class needs to provide the following constants:
140+
// - MID_BITS : number of bits after decimal points used for mid
141+
// - MID_MASK : 2^MID_BITS - 1, mask to extract mid bits
142+
// - LOG2_B : log2(b) * 2^MID_BITS for scaling
143+
// - M_LOGB_2_HI : high part of -log_b(2) * 2^(-MID_BITS)
144+
// - M_LOGB_2_LO : low part of -log_b(2) * 2^(-MID_BITS)
145+
// - EXP_2_MID : look up table for bit fields of 2^mid
146+
// Return:
147+
// { 2^(hi + mid), lo }
148+
template <class Base>
149+
LIBC_INLINE static constexpr exp_b_reduc_t exp_b_range_reduc(float x) {
150+
double xd = static_cast<double>(x);
151+
// kd = round((hi + mid) * log2(b) * 2^MID_BITS)
152+
double kd = fputil::nearest_integer(Base::LOG2_B * xd);
153+
// k = round((hi + mid) * log2(b) * 2^MID_BITS)
154+
int k = static_cast<int>(kd);
155+
// hi = floor(kd * 2^(-MID_BITS))
156+
// exp_hi = shift hi to the exponent field of double precision.
157+
uint64_t exp_hi = static_cast<uint64_t>(k >> Base::MID_BITS)
158+
<< fputil::FPBits<double>::FRACTION_LEN;
159+
// mh = 2^hi * 2^mid
160+
// mh_bits = bit field of mh
161+
uint64_t mh_bits = Base::EXP_2_MID[k & Base::MID_MASK] + exp_hi;
162+
double mh = fputil::FPBits<double>(mh_bits).get_val();
163+
// dx = lo = x - (hi + mid) * log(2)
164+
double dx = fputil::multiply_add(
165+
kd, Base::M_LOGB_2_LO, fputil::multiply_add(kd, Base::M_LOGB_2_HI, xd));
166+
return {mh, dx};
167+
}
168+
169+
} // namespace LIBC_NAMESPACE_DECL
170+
171+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_FLOAT_CONSTANTS_H

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

Lines changed: 6 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -358,7 +358,6 @@ add_entrypoint_object(
358358
libc.src.__support.FPUtil.fp_bits
359359
libc.src.__support.FPUtil.except_value_utils
360360
libc.src.__support.FPUtil.fma
361-
libc.src.__support.FPUtil.multiply_add
362361
libc.src.__support.FPUtil.polyeval
363362
libc.src.__support.macros.optimization
364363
)
@@ -448,7 +447,6 @@ add_entrypoint_object(
448447
libc.src.__support.FPUtil.fenv_impl
449448
libc.src.__support.FPUtil.fp_bits
450449
libc.src.__support.FPUtil.fma
451-
libc.src.__support.FPUtil.multiply_add
452450
libc.src.__support.FPUtil.polyeval
453451
libc.src.__support.FPUtil.rounding_mode
454452
libc.src.__support.macros.optimization
@@ -1461,29 +1459,15 @@ add_entrypoint_object(
14611459
libc.src.errno.errno
14621460
)
14631461

1464-
add_header_library(
1465-
exp10f_impl
1466-
HDRS
1467-
exp10f_impl.h
1468-
DEPENDS
1469-
.explogxf
1470-
libc.src.__support.FPUtil.fenv_impl
1471-
libc.src.__support.FPUtil.fp_bits
1472-
libc.src.__support.FPUtil.multiply_add
1473-
libc.src.__support.FPUtil.rounding_mode
1474-
libc.src.__support.macros.optimization
1475-
libc.src.__support.common
1476-
libc.src.errno.errno
1477-
)
1478-
14791462
add_entrypoint_object(
14801463
exp10f
14811464
SRCS
14821465
exp10f.cpp
14831466
HDRS
14841467
../exp10f.h
14851468
DEPENDS
1486-
.exp10f_impl
1469+
libc.src.__support.math.exp10f
1470+
libc.src.errno.errno
14871471
)
14881472

14891473
add_entrypoint_object(
@@ -1620,17 +1604,15 @@ add_entrypoint_object(
16201604
../powf.h
16211605
DEPENDS
16221606
.common_constants
1623-
.exp10f_impl
16241607
.exp2f_impl
16251608
.explogxf
1609+
libc.src.__support.math.exp10f
16261610
libc.src.__support.CPP.bit
1627-
libc.src.__support.CPP.optional
16281611
libc.src.__support.FPUtil.fenv_impl
16291612
libc.src.__support.FPUtil.fp_bits
16301613
libc.src.__support.FPUtil.multiply_add
16311614
libc.src.__support.FPUtil.nearest_integer
16321615
libc.src.__support.FPUtil.polyeval
1633-
libc.src.__support.FPUtil.rounding_mode
16341616
libc.src.__support.FPUtil.sqrt
16351617
libc.src.__support.FPUtil.triple_double
16361618
libc.src.__support.macros.optimization
@@ -3792,13 +3774,9 @@ add_object_library(
37923774
explogxf.cpp
37933775
DEPENDS
37943776
.common_constants
3795-
libc.src.__support.FPUtil.basic_operations
3796-
libc.src.__support.FPUtil.fenv_impl
3797-
libc.src.__support.FPUtil.multiply_add
3798-
libc.src.__support.FPUtil.nearest_integer
3799-
libc.src.__support.FPUtil.polyeval
3800-
libc.src.__support.common
38013777
libc.src.__support.math.exp_utils
3778+
libc.src.__support.math.exp10f_utils
3779+
libc.src.__support.macros.properties.cpu_features
38023780
libc.src.errno.errno
38033781
)
38043782

@@ -3981,6 +3959,7 @@ add_entrypoint_object(
39813959
DEPENDS
39823960
.explogxf
39833961
libc.src.__support.FPUtil.fp_bits
3962+
libc.src.__support.FPUtil.fenv_impl
39843963
libc.src.__support.macros.optimization
39853964
)
39863965

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/atanhf.h"
10+
#include "src/__support/FPUtil/FEnvImpl.h"
1011
#include "src/__support/FPUtil/FPBits.h"
1112
#include "src/__support/macros/config.h"
1213
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,8 @@
77
//===----------------------------------------------------------------------===//
88

99
#include "src/math/coshf.h"
10+
#include "src/__support/FPUtil/FEnvImpl.h"
1011
#include "src/__support/FPUtil/FPBits.h"
11-
#include "src/__support/FPUtil/multiply_add.h"
1212
#include "src/__support/FPUtil/rounding_mode.h"
1313
#include "src/__support/macros/config.h"
1414
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY

0 commit comments

Comments
 (0)