Skip to content

Commit da28593

Browse files
lntuelntue
authored
[libc][math] Implement double precision expm1 function correctly rounded for all rounding modes. (#67048)
Implementing expm1 function for double precision based on exp function algorithm: - Reduced x = log2(e) * (hi + mid1 + mid2) + lo, where: * hi is an integer * mid1 * 2^-6 is an integer * mid2 * 2^-12 is an integer * |lo| < 2^-13 + 2^-30 - Then exp(x) - 1 = 2^hi * 2^mid1 * 2^mid2 * exp(lo) - 1 ~ 2^hi * (2^mid1 * 2^mid2 * (1 + lo * P(lo)) - 2^(-hi) ) - We evaluate fast pass with P(lo) is a degree-3 Taylor polynomial of (e^lo - 1) / lo in double precision - If the Ziv accuracy test fails, we use degree-6 Taylor polynomial of (e^lo - 1) / lo in double double precision - If the Ziv accuracy test still fails, we re-evaluate everything in 128-bit precision.
1 parent d28a782 commit da28593

File tree

21 files changed

+781
-16
lines changed

21 files changed

+781
-16
lines changed

‎libc/config/darwin/arm/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -135,6 +135,7 @@ set(TARGET_LIBM_ENTRYPOINTS
135135
libc.src.math.exp10f
136136
libc.src.math.exp2
137137
libc.src.math.exp2f
138+
libc.src.math.expm1
138139
libc.src.math.expm1f
139140
libc.src.math.fabs
140141
libc.src.math.fabsf

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -252,6 +252,7 @@ set(TARGET_LIBM_ENTRYPOINTS
252252
libc.src.math.exp10f
253253
libc.src.math.exp2
254254
libc.src.math.exp2f
255+
libc.src.math.expm1
255256
libc.src.math.expm1f
256257
libc.src.math.fabs
257258
libc.src.math.fabsf

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -261,6 +261,7 @@ set(TARGET_LIBM_ENTRYPOINTS
261261
libc.src.math.exp10f
262262
libc.src.math.exp2
263263
libc.src.math.exp2f
264+
libc.src.math.expm1
264265
libc.src.math.expm1f
265266
libc.src.math.fabs
266267
libc.src.math.fabsf

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -265,6 +265,7 @@ set(TARGET_LIBM_ENTRYPOINTS
265265
libc.src.math.exp10f
266266
libc.src.math.exp2
267267
libc.src.math.exp2f
268+
libc.src.math.expm1
268269
libc.src.math.expm1f
269270
libc.src.math.fabs
270271
libc.src.math.fabsf

‎libc/config/windows/entrypoints.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -134,6 +134,7 @@ set(TARGET_LIBM_ENTRYPOINTS
134134
libc.src.math.exp10f
135135
libc.src.math.exp2
136136
libc.src.math.exp2f
137+
libc.src.math.expm1
137138
libc.src.math.expm1f
138139
libc.src.math.fabs
139140
libc.src.math.fabsf

‎libc/docs/math/index.rst

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -370,7 +370,7 @@ Higher Math Functions
370370
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
371371
| exp2l | | | | | | | | | | | | |
372372
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
373-
| expm1 | | | | | | | | | | | | |
373+
| expm1 | |check| | |check| | | |check| | |check| | | | |check| | | | | |
374374
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
375375
| expm1f | |check| | |check| | | |check| | |check| | | | |check| | | | | |
376376
+------------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+---------+
@@ -484,9 +484,9 @@ cos |check| large
484484
cosh |check|
485485
erf |check|
486486
exp |check| |check|
487-
exp10 |check|
487+
exp10 |check| |check|
488488
exp2 |check| |check|
489-
expm1 |check|
489+
expm1 |check| |check|
490490
fma |check| |check|
491491
hypot |check| |check|
492492
log |check| |check|

‎libc/spec/stdc.td

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -439,6 +439,8 @@ def StdC : StandardSpec<"stdc"> {
439439

440440
FunctionSpec<"exp2", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
441441
FunctionSpec<"exp2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
442+
443+
FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
442444
FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
443445

444446
FunctionSpec<"remainderf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,

‎libc/src/__support/FPUtil/FPBits.h

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ template <typename T> struct FPBits {
116116

117117
FPBits() : bits(0) {}
118118

119-
LIBC_INLINE T get_val() const { return cpp::bit_cast<T>(bits); }
119+
LIBC_INLINE constexpr T get_val() const { return cpp::bit_cast<T>(bits); }
120120

121121
LIBC_INLINE void set_val(T value) { bits = cpp::bit_cast<UIntType>(value); }
122122

@@ -185,6 +185,10 @@ template <typename T> struct FPBits {
185185
return bits;
186186
}
187187

188+
LIBC_INLINE static constexpr FPBits<T> min_normal() {
189+
return FPBits<T>(MIN_NORMAL);
190+
}
191+
188192
LIBC_INLINE static constexpr T build_nan(UIntType v) {
189193
FPBits<T> bits = inf();
190194
bits.set_mantissa(v);

‎libc/src/__support/FPUtil/except_value_utils.h

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -100,15 +100,17 @@ template <typename T, size_t N> struct ExceptValues {
100100
};
101101

102102
// Helper functions to set results for exceptional cases.
103-
LIBC_INLINE float round_result_slightly_down(float value_rn) {
104-
volatile float tmp = value_rn;
105-
tmp = tmp - 0x1.0p-100f;
103+
template <typename T> LIBC_INLINE T round_result_slightly_down(T value_rn) {
104+
volatile T tmp = value_rn;
105+
constexpr T MIN_NORMAL = FPBits<T>::min_normal().get_val();
106+
tmp = tmp - MIN_NORMAL;
106107
return tmp;
107108
}
108109

109-
LIBC_INLINE float round_result_slightly_up(float value_rn) {
110-
volatile float tmp = value_rn;
111-
tmp = tmp + 0x1.0p-100f;
110+
template <typename T> LIBC_INLINE T round_result_slightly_up(T value_rn) {
111+
volatile T tmp = value_rn;
112+
const T MIN_NORMAL = FPBits<T>::min_normal().get_val();
113+
tmp = tmp + MIN_NORMAL;
112114
return tmp;
113115
}
114116

‎libc/src/math/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@ add_math_entrypoint_object(exp2f)
8888
add_math_entrypoint_object(exp10)
8989
add_math_entrypoint_object(exp10f)
9090

91+
add_math_entrypoint_object(expm1)
9192
add_math_entrypoint_object(expm1f)
9293

9394
add_math_entrypoint_object(fabs)

‎libc/src/math/expm1.h

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
//===-- Implementation header for expm1 -------------------------*- 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_EXPM1_H
10+
#define LLVM_LIBC_SRC_MATH_EXPM1_H
11+
12+
namespace LIBC_NAMESPACE {
13+
14+
double expm1(double x);
15+
16+
} // namespace LIBC_NAMESPACE
17+
18+
#endif // LLVM_LIBC_SRC_MATH_EXPM1_H

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

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -697,6 +697,33 @@ add_entrypoint_object(
697697
-O3
698698
)
699699

700+
add_entrypoint_object(
701+
expm1
702+
SRCS
703+
expm1.cpp
704+
HDRS
705+
../expm1.h
706+
DEPENDS
707+
.common_constants
708+
.explogxf
709+
libc.src.__support.CPP.bit
710+
libc.src.__support.CPP.optional
711+
libc.src.__support.FPUtil.dyadic_float
712+
libc.src.__support.FPUtil.fenv_impl
713+
libc.src.__support.FPUtil.fp_bits
714+
libc.src.__support.FPUtil.multiply_add
715+
libc.src.__support.FPUtil.nearest_integer
716+
libc.src.__support.FPUtil.polyeval
717+
libc.src.__support.FPUtil.rounding_mode
718+
libc.src.__support.FPUtil.triple_double
719+
libc.src.__support.macros.optimization
720+
libc.include.errno
721+
libc.src.errno.errno
722+
libc.include.math
723+
COMPILE_OPTIONS
724+
-O3
725+
)
726+
700727
add_entrypoint_object(
701728
expm1f
702729
SRCS

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -104,7 +104,7 @@ Float128 poly_approx_f128(const Float128 &dx) {
104104
{false, -132, MType({0xaaaaaaaaaaaaaaab, 0xaaaaaaaaaaaaaaaa})}, // 1/24
105105
{false, -134, MType({0x8888888888888889, 0x8888888888888888})}, // 1/120
106106
{false, -137, MType({0x60b60b60b60b60b6, 0xb60b60b60b60b60b})}, // 1/720
107-
{false, -140, MType({0x00b00b00b00b00b0, 0xb00b00b00b00b00b})}, // 1/5040
107+
{false, -140, MType({0x00d00d00d00d00d0, 0xd00d00d00d00d00d})}, // 1/5040
108108
};
109109

110110
Float128 p = fputil::polyeval(dx, COEFFS_128[0], COEFFS_128[1], COEFFS_128[2],

0 commit comments

Comments
 (0)