1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
//===-- Single-precision e^x - 1 function ---------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#include "src/math/expm1f.h"
#include "common_constants.h" // Lookup tables EXP_M1 and EXP_M2.
#include "src/__support/FPUtil/BasicOperations.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FMA.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
uint32_t x_u = xbits.uintval();
uint32_t x_abs = x_u & 0x7fff'ffffU;
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Exceptional value
if (LIBC_UNLIKELY(x_u == 0x3e35'bec5U)) { // x = 0x1.6b7d8ap-3f
int round_mode = fputil::quick_get_round();
if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
return 0x1.8dbe64p-3f;
return 0x1.8dbe62p-3f;
}
#if !defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
if (LIBC_UNLIKELY(x_u == 0xbdc1'c6cbU)) { // x = -0x1.838d96p-4f
int round_mode = fputil::quick_get_round();
if (round_mode == FE_TONEAREST || round_mode == FE_DOWNWARD)
return -0x1.71c884p-4f;
return -0x1.71c882p-4f;
}
#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// When |x| > 25*log(2), or nan
if (LIBC_UNLIKELY(x_abs >= 0x418a'a123U)) {
// x < log(2^-25)
if (xbits.is_neg()) {
// exp(-Inf) = 0
if (xbits.is_inf())
return -1.0f;
// exp(nan) = nan
if (xbits.is_nan())
return x;
int round_mode = fputil::quick_get_round();
if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO)
return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
return -1.0f;
} else {
// x >= 89 or nan
if (xbits.uintval() >= 0x42b2'0000) {
if (xbits.uintval() < 0x7f80'0000U) {
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
}
return x + FPBits::inf().get_val();
}
}
}
// |x| < 2^-4
if (x_abs < 0x3d80'0000U) {
// |x| < 2^-25
if (x_abs < 0x3300'0000U) {
// x = -0.0f
if (LIBC_UNLIKELY(xbits.uintval() == 0x8000'0000U))
return x;
// When |x| < 2^-25, the relative error of the approximation e^x - 1 ~ x
// is:
// |(e^x - 1) - x| / |e^x - 1| < |x^2| / |x|
// = |x|
// < 2^-25
// < epsilon(1)/2.
// So the correctly rounded values of expm1(x) are:
// = x + eps(x) if rounding mode = FE_UPWARD,
// or (rounding mode = FE_TOWARDZERO and x is
// negative),
// = x otherwise.
// To simplify the rounding decision and make it more efficient, we use
// fma(x, x, x) ~ x + x^2 instead.
// Note: to use the formula x + x^2 to decide the correct rounding, we
// do need fma(x, x, x) to prevent underflow caused by x*x when |x| <
// 2^-76. For targets without FMA instructions, we simply use double for
// intermediate results as it is more efficient than using an emulated
// version of FMA.
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
return fputil::multiply_add(x, x, x);
#else
double xd = x;
return static_cast<float>(fputil::multiply_add(xd, xd, xd));
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
}
constexpr double COEFFS[] = {0x1p-1,
0x1.55555555557ddp-3,
0x1.55555555552fap-5,
0x1.111110fcd58b7p-7,
0x1.6c16c1717660bp-10,
0x1.a0241f0006d62p-13,
0x1.a01e3f8d3c06p-16};
// 2^-25 <= |x| < 2^-4
double xd = static_cast<double>(x);
double xsq = xd * xd;
// Degree-8 minimax polynomial generated by Sollya with:
// > display = hexadecimal;
// > P = fpminimax((expm1(x) - x)/x^2, 6, [|D...|], [-2^-4, 2^-4]);
double c0 = fputil::multiply_add(xd, COEFFS[1], COEFFS[0]);
double c1 = fputil::multiply_add(xd, COEFFS[3], COEFFS[2]);
double c2 = fputil::multiply_add(xd, COEFFS[5], COEFFS[4]);
double r = fputil::polyeval(xsq, c0, c1, c2, COEFFS[6]);
return static_cast<float>(fputil::multiply_add(r, xsq, xd));
}
// For -18 < x < 89, to compute expm1(x), we perform the following range
// reduction: find hi, mid, lo such that:
// x = hi + mid + lo, in which
// hi is an integer,
// mid * 2^7 is an integer
// -2^(-8) <= lo < 2^-8.
// In particular,
// hi + mid = round(x * 2^7) * 2^(-7).
// Then,
// expm1(x) = exp(hi + mid + lo) - 1 = exp(hi) * exp(mid) * exp(lo) - 1.
// We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2
// respectively. exp(lo) is computed using a degree-4 minimax polynomial
// generated by Sollya.
// x_hi = hi + mid.
float kf = fputil::nearest_integer(x * 0x1.0p7f);
int x_hi = static_cast<int>(kf);
// Subtract (hi + mid) from x to get lo.
double xd = static_cast<double>(fputil::multiply_add(kf, -0x1.0p-7f, x));
x_hi += 104 << 7;
// hi = x_hi >> 7
double exp_hi = EXP_M1[x_hi >> 7];
// lo = x_hi & 0x0000'007fU;
double exp_mid = EXP_M2[x_hi & 0x7f];
double exp_hi_mid = exp_hi * exp_mid;
// Degree-4 minimax polynomial generated by Sollya with the following
// commands:
// > display = hexadecimal;
// > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]);
// > Q;
double exp_lo =
fputil::polyeval(xd, 0x1.0p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1,
0x1.555566668e5e7p-3, 0x1.55555555ef243p-5);
return static_cast<float>(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0));
}
} // namespace LIBC_NAMESPACE_DECL
|
