1//===----------------------------------------------------------------------===//
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/// \file
10/// Shared utilities and tables for double-precision exponential functions.
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
16
17#include "hdr/errno_macros.h"
18#include "hdr/fenv_macros.h"
19#include "src/__support/CPP/bit.h"
20#include "src/__support/CPP/optional.h"
21#include "src/__support/FPUtil/FEnvImpl.h"
22#include "src/__support/FPUtil/FPBits.h"
23#include "src/__support/FPUtil/double_double.h"
24#include "src/__support/FPUtil/multiply_add.h"
25#include "src/__support/FPUtil/rounding_mode.h"
26#include "src/__support/common.h"
27#include "src/__support/macros/attributes.h"
28#include "src/__support/macros/config.h"
29#include "src/__support/macros/optimization.h"
30#include "src/__support/macros/properties/architectures.h"
31#include "src/__support/macros/properties/cpu_features.h"
32
33#ifdef LIBC_TARGET_CPU_HAS_SSE2
34#include <immintrin.h>
35#elif defined(LIBC_TARGET_ARCH_IS_AARCH64) && defined(__ARM_NEON)
36#include <arm_neon.h>
37#endif // LIBC_TARGET_CPU_HAS_SSE2
38
39#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
40#include "src/__support/math/exp2_accurate_128.h"
41#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
42
43namespace LIBC_NAMESPACE_DECL {
44
45// Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We
46// assume further that 1 <= mid < 2, mid + lo < 2, and |lo| << mid.
47// Notice that, if 0 < x < 2^-1022,
48// double(2^-1022 + x) - 2^-1022 = double(x).
49// So if we scale x up by 2^1022, we can use
50// double(1.0 + 2^1022 * x) - 1.0 to test how x is rounded in denormal range.
51template <bool SKIP_ZIV_TEST = false>
52LIBC_INLINE constexpr cpp::optional<double>
53ziv_test_denorm(int hi, double mid, double lo, double err) {
54 using FPBits = typename fputil::FPBits<double>;
55
56 // Scaling factor = 1/(min normal number) = 2^1022
57 int64_t exp_hi =
58 static_cast<int64_t>(hi + 1022) * (1LL << FPBits::FRACTION_LEN);
59 double mid_hi = cpp::bit_cast<double>(from: exp_hi + cpp::bit_cast<int64_t>(from: mid));
60 double lo_scaled =
61 (lo != 0.0) ? cpp::bit_cast<double>(from: exp_hi + cpp::bit_cast<int64_t>(from: lo))
62 : 0.0;
63
64 double extra_factor = 0.0;
65 uint64_t scale_down = 0x3FE0'0000'0000'0000; // 1022 in the exponent field.
66
67 // Result is denormal if (mid_hi + lo_scale < 1.0).
68 if ((1.0 - mid_hi) > lo_scaled) {
69 // Extra rounding step is needed, which adds more rounding errors.
70 err += 0x1.0p-52;
71 extra_factor = 1.0;
72 scale_down = 0x3FF0'0000'0000'0000; // 1023 in the exponent field.
73 }
74
75 // By adding 1.0, the results will have similar rounding points as denormal
76 // outputs.
77 if constexpr (SKIP_ZIV_TEST) {
78 double r = extra_factor + (mid_hi + lo_scaled);
79 return cpp::bit_cast<double>(from: cpp::bit_cast<uint64_t>(from: r) - scale_down);
80 } else {
81 double err_scaled =
82 cpp::bit_cast<double>(from: exp_hi + cpp::bit_cast<int64_t>(from: err));
83
84 double lo_u = lo_scaled + err_scaled;
85 double lo_l = lo_scaled - err_scaled;
86
87 double upper = extra_factor + (mid_hi + lo_u);
88 double lower = extra_factor + (mid_hi + lo_l);
89
90 if (LIBC_LIKELY(upper == lower)) {
91 return cpp::bit_cast<double>(from: cpp::bit_cast<uint64_t>(from: upper) - scale_down);
92 }
93
94 return cpp::nullopt;
95 }
96}
97
98namespace math {
99namespace double_eval {
100
101using fputil::DoubleDouble;
102
103// log(2) in double-double precision, generated by Sollya with:
104// > display = hexadecimal!;
105// > prec = 500;
106// > hi = round(log(2), D, RN);
107// > lo = round(log(2) - hi, D, RN);
108// > print("{", lo, ",", hi, "};");
109LIBC_INLINE_VAR constexpr DoubleDouble LOG_2 = {.lo: 0x1.abc9e3b39803fp-56,
110 .hi: 0x1.62e42fefa39efp-1};
111
112// Directly add an integer `hi` to the exponent field of a double precision `r`.
113// It is equivalent to scaling `r` by 2^hi when both `r` and the resulting bits
114// are normal.
115// The default implementation is simply bit cast `r` to int64`, shift `hi` to
116// the exponent field, then add `r` bits and shifted `hi`, before bit cast the
117// result back to double.
118// The addition is done in integer, and most likely generated as addition in the
119// general purpose registers (GPR), causing a full round trip from
120// floating-point/SIMD register -> GPR -> FP/SIMD.
121// When SSE2 or NEON are available, we can use the compiler intrinsic to keep
122// `r` in the FP/SIMD register, and perform the integer addition in the FP/SIMD
123// registers instead. This in general results in a reduction in latency since
124// `r` is normally on the critical path.
125LIBC_ALWAYS_INLINE double scale_exp2(double r, int hi) {
126 using FPBits = typename fputil::FPBits<double>;
127#ifdef LIBC_TARGET_CPU_HAS_SSE2
128 // Move `hi` to SIMD register and shift to exponent field position.
129 __m128i v_exp = _mm_slli_epi64(a: _mm_cvtsi32_si128(a: hi), count: FPBits::FRACTION_LEN);
130 // Load `r` and "bit_cast" it to integer SIMD type.
131 __m128i v_r = _mm_castpd_si128(a: _mm_set_sd(w: r));
132 // Perform integer addition in SIMD registers and bit cast the result back to
133 // double.
134 return _mm_cvtsd_f64(a: _mm_castsi128_pd(a: _mm_add_epi64(a: v_exp, b: v_r)));
135#elif defined(LIBC_TARGET_ARCH_IS_AARCH64) && defined(__ARM_NEON)
136 // Shift `hi` to exponent field position and move it to SIMD register.
137 int64x1_t v_exp =
138 vdup_n_s64(static_cast<int64_t>(hi) << FPBits::FRACTION_LEN);
139 // Bit cast `r` to integer SIMD type.
140 int64x1_t v_r = vreinterpret_s64_f64(vdup_n_f64(r));
141 // Perform integer addition in SIMD registers and bit cast the result back to
142 // double.
143 return vget_lane_f64(vreinterpret_f64_s64(vadd_s64(v_exp, v_r)), 0);
144#else
145 int64_t exp_hi = static_cast<int64_t>(hi) * (1LL << FPBits::FRACTION_LEN);
146 return cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(r));
147#endif // LIBC_TARGET_CPU_HAS_SSE2
148}
149
150// Check for exceptional cases (tiny, overflow, underflow, inf, or nan) for
151// exp, exp2, and exp10. Callers have already verified that x is in the exact
152// tiny, underflow, or overflow/inf/nan range for the specific function before
153// calling this helper.
154LIBC_NOINLINE LIBC_INLINE double set_exceptional(double x) {
155 using FPBits = typename fputil::FPBits<double>;
156 FPBits xbits(x);
157
158 // |x| <= 0x1.0p-53 (caller has already checked the exact tiny bound):
159 constexpr uint64_t TINY_BOUND = 0x3ca0'0000'0000'0000ULL;
160 if (xbits.abs().uintval() <= TINY_BOUND) {
161 // b^x ~ 1 + x
162 return 1.0 + x;
163 }
164
165 // Negative underflow or -inf/nan (caller has already excluded
166 // non-underflowing negative inputs):
167 if (xbits.is_neg()) {
168 // b^(-Inf) = 0
169 if (xbits.is_inf())
170 return 0.0;
171
172 // b^(nan) = nan
173 // Adding 0.0 to silence signaling nan.
174 if (xbits.is_nan())
175 return x + 0.0;
176
177 fputil::set_errno_if_required(ERANGE);
178 fputil::raise_underflow_except_if_required<double>();
179#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
180 if (fputil::quick_get_round() == FE_UPWARD)
181 return FPBits::min_subnormal().get_val();
182#endif
183 return 0.0;
184 }
185
186 // Positive overflow or +inf/nan (caller has already checked the exact
187 // overflow bound)
188
189 // x is finite and positive
190 if (xbits.uintval() < FPBits::inf().uintval()) {
191 fputil::set_errno_if_required(ERANGE);
192 fputil::raise_overflow_except_if_required<double>();
193#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
194 int rounding = fputil::quick_get_round();
195 if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
196 return FPBits::max_normal().get_val();
197#endif
198 return FPBits::inf().get_val();
199 }
200 // x is +inf or nan
201 // Adding `inf` to silence signaling nan.
202 return x + FPBits::inf().get_val();
203}
204
205#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) || \
206 !defined(LIBC_MATH_HAS_SMALL_TABLES)
207
208LIBC_INLINE_VAR constexpr int EXP_N_BITS = 7;
209LIBC_INLINE_VAR constexpr double EXP_SCALE = 128.0; // 2^EXP_N_BITS = 2^7
210LIBC_INLINE_VAR constexpr double EXP_SCALE_INV = 0x1.0p-7;
211LIBC_INLINE_VAR constexpr int EXP_MASK = 0x7f;
212
213// Polynomial approximation for (exp(x) - 1 - x) / x^2 on
214// [-B, B] with B = 2^-8 * log(2):
215//
216// 1) A degree-3 minimax polynomial P0 for (exp(x) - 1 - x) / x^2 on [-B, B]
217// (which minimizes the error of exp(x) - (1 + x + x^2 * P0(x)) weighted by
218// 1 / x^2, i.e., relative error with respect to x^2):
219// > display = hexadecimal!;
220// > prec = 500;
221// > B = 2^(-8) * log(2);
222// > P0 = fpminimax((exp(x) - 1 - x)/x^2, 3, [|D...|], [-B, B]);
223// > dirtyinfnorm(exp(x) - (1 + x + x^2 * P0), [-B, B]);
224// 0x1.4362...p-64 < 2^-63.66.
225//
226// Some explanation for the errors (from Chebyshev polynomials):
227// Because P0 minimizes |(exp(x) - 1 - x) / x^2 - P0(x)|, its error is governed
228// by the degree-4 Chebyshev polynomial:
229// T_4(x / B) = 8*(x / B)^4 - 8*(x / B)^2 + 1
230// approximating the leading omitted Taylor term x^4 / 720, with magnitude:
231// B^4 / (720 * 8).
232// Multiplying by x^2 gives an absolute error of:
233// x^2 * B^4 / (720 * 8) * T_4(x / B),
234// peaking at x = +-B with:
235// B^6 / (720 * 8) ~ 0x1.4309...p-64 < 2^-63.66.
236//
237// 2) If we instead minimize the absolute error of exp(x) - (1 + x + P(x))
238// directly on [-B, B] with monomials [|2, 3, 4, 5|]:
239// > P = fpminimax(exp(x) - 1 - x, [|2, 3, 4, 5|], [|D...|], [-B, B],
240// absolute);
241// > dirtyinfnorm(exp(x) - (1 + x + P), [-B, B]);
242// 0x1.8e05...p-66 < 2^-65.36.
243//
244// Some explanation for the better errors from Chebyshev polynomials:
245// The free x^2 and x^4 coefficients in P(x) approximate the omitted x^6 / 720
246// term on [-B, B] subject to vanishing at x = 0. Setting t = (x / B)^2 in
247// [0, 1], this is equivalent to approximating t^3 on [0, 1] by
248// a_2 * t^2 + a_1 * t vanishing at t = 0.
249// The degree-3 Chebyshev polynomial T_3(z) = 4*z^3 - 3*z has a root at
250// z = -sqrt(3)/2 and three alternating extrema +-1 on [-sqrt(3)/2, 1] at
251// z = -1/2, 1/2, 1.
252// Mapping t in [0, 1] affinely onto z in [-sqrt(3)/2, 1] via:
253// z = ((2 + sqrt(3)) / 2) * t - sqrt(3) / 2
254// gives an error curve T_3(z) that vanishes at t = 0 and equioscillates three
255// times on (0, 1], with leading coefficient:
256// 4 * ((2 + sqrt(3)) / 2)^3 = (2 + sqrt(3))^3 / 2.
257// Dividing by this leading coefficient gives a maximum error on [-B, B] of:
258// (2 / (2 + sqrt(3))^3) * B^6 / 720 = 2 * (2 - sqrt(3))^3 * B^6 / 720
259// ~ 0x1.8dba...p-66 < 2^-65.36,
260// reducing the maximum error compared to P0 by (2 + sqrt(3))^3 / 16 ~ 3.25x.
261LIBC_INLINE_VAR constexpr double EXP_P_COEFFS[4] = {
262 0x1.ffffffffffdbdp-2, // p_0
263 0x1.555555555543dp-3, // p_1
264 0x1.55555cf16d6eep-5, // p_2
265 0x1.111116782220ap-7, // p_3
266};
267
268// Evaluate P(u) ~ (e^u - 1 - u) / u^2 using Estrin's scheme:
269LIBC_ALWAYS_INLINE double exp_poly(double u, double u2) {
270 double c0 = fputil::multiply_add(x: u, y: EXP_P_COEFFS[1], z: EXP_P_COEFFS[0]);
271 double c1 = fputil::multiply_add(x: u, y: EXP_P_COEFFS[3], z: EXP_P_COEFFS[2]);
272 return fputil::multiply_add(x: u2, y: c1, z: c0);
273}
274
275#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS || !LIBC_MATH_HAS_SMALL_TABLES
276
277#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
278
279// Fixed-point splitting constants (ulp = 2^-22, 2^-18) and pre-split
280// 12-bit + 53-bit logarithm constants for exact non-FMA multiplication without
281// Veltkamp splitting:
282LIBC_INLINE_VAR constexpr double SPLIT_22 = 0x1.8p30; // 1.5 * 2^(-22 + 52)
283LIBC_INLINE_VAR constexpr double SPLIT_18 = 0x1.8p34; // 1.5 * 2^(-18 + 52)
284
285// Splitting log(2) ~ hi + lo with 12-bit-precision hi part.
286// > hi = round(log(2), 12, RN);
287// > lo = round(log(2) - hi, D, RN);
288LIBC_INLINE_VAR constexpr double LOG_2_HI_12 = 0x1.62ep-1;
289LIBC_INLINE_VAR constexpr double LOG_2_LO_12 = 0x1.0bfbe8e7bcd5ep-15;
290
291// Error bound for double precision fast path (N = 7, Ziv failure rate <
292// 0.05%):
293LIBC_INLINE_VAR constexpr double EXP2_ERR_D = 0x1.0p-64;
294
295// Lookup table for 2^(i * 2^-7) with i = 0..127.
296// Generated by Sollya with:
297// > for i from 0 to 127 do {
298// a = 2^(i * 2^-7);
299// b = round(a, D, RN);
300// c = round(a - b, D, RN);
301// print("{", c, ",", b, "},");
302// };
303LIBC_INLINE_VAR constexpr DoubleDouble EXP2_MID[128] = {
304 {.lo: 0, .hi: 0x1p0},
305 {.lo: 0x1.b61299ab8cdb7p-54, .hi: 0x1.0163da9fb3335p0},
306 {.lo: -0x1.19083535b085dp-56, .hi: 0x1.02c9a3e778061p0},
307 {.lo: -0x1.0a31c1977c96ep-54, .hi: 0x1.04315e86e7f85p0},
308 {.lo: 0x1.d73e2a475b465p-55, .hi: 0x1.059b0d3158574p0},
309 {.lo: -0x1.c91dfe2b13c27p-55, .hi: 0x1.0706b29ddf6dep0},
310 {.lo: 0x1.186be4bb284ffp-57, .hi: 0x1.0874518759bc8p0},
311 {.lo: 0x1.1487818316136p-54, .hi: 0x1.09e3ecac6f383p0},
312 {.lo: 0x1.8a62e4adc610bp-54, .hi: 0x1.0b5586cf9890fp0},
313 {.lo: 0x1.01edc16e24f71p-54, .hi: 0x1.0cc922b7247f7p0},
314 {.lo: 0x1.03a1727c57b53p-59, .hi: 0x1.0e3ec32d3d1a2p0},
315 {.lo: -0x1.b9bedc44ebd7bp-57, .hi: 0x1.0fb66affed31bp0},
316 {.lo: -0x1.6c51039449b3ap-54, .hi: 0x1.11301d0125b51p0},
317 {.lo: -0x1.1b514b36ca5c7p-58, .hi: 0x1.12abdc06c31ccp0},
318 {.lo: -0x1.32fbf9af1369ep-54, .hi: 0x1.1429aaea92dep0},
319 {.lo: 0x1.2406ab9eeab0ap-55, .hi: 0x1.15a98c8a58e51p0},
320 {.lo: -0x1.19041b9d78a76p-55, .hi: 0x1.172b83c7d517bp0},
321 {.lo: -0x1.11023d1970f6cp-54, .hi: 0x1.18af9388c8deap0},
322 {.lo: 0x1.e5b4c7b4968e4p-55, .hi: 0x1.1a35beb6fcb75p0},
323 {.lo: -0x1.95386352ef607p-54, .hi: 0x1.1bbe084045cd4p0},
324 {.lo: 0x1.e016e00a2643cp-54, .hi: 0x1.1d4873168b9aap0},
325 {.lo: -0x1.1df98027bb78cp-54, .hi: 0x1.1ed5022fcd91dp0},
326 {.lo: 0x1.dc775814a8495p-55, .hi: 0x1.2063b88628cd6p0},
327 {.lo: 0x1.2a97e9494a5eep-55, .hi: 0x1.21f49917ddc96p0},
328 {.lo: 0x1.9b07eb6c70573p-54, .hi: 0x1.2387a6e756238p0},
329 {.lo: 0x1.ac155bef4f4a4p-55, .hi: 0x1.251ce4fb2a63fp0},
330 {.lo: 0x1.2bd339940e9d9p-55, .hi: 0x1.26b4565e27cddp0},
331 {.lo: -0x1.a4c3a8c3f0d7ep-54, .hi: 0x1.284dfe1f56381p0},
332 {.lo: 0x1.612e8afad1255p-55, .hi: 0x1.29e9df51fdee1p0},
333 {.lo: -0x1.10adcd6381aa4p-59, .hi: 0x1.2b87fd0dad99p0},
334 {.lo: 0x1.0024754db41d5p-54, .hi: 0x1.2d285a6e4030bp0},
335 {.lo: 0x1.1ca0f45d52383p-56, .hi: 0x1.2ecafa93e2f56p0},
336 {.lo: 0x1.6f46ad23182e4p-55, .hi: 0x1.306fe0a31b715p0},
337 {.lo: 0x1.a9ce78e18047cp-55, .hi: 0x1.32170fc4cd831p0},
338 {.lo: 0x1.32721843659a6p-54, .hi: 0x1.33c08b26416ffp0},
339 {.lo: -0x1.b5cee5c4e4628p-55, .hi: 0x1.356c55f929ff1p0},
340 {.lo: -0x1.63aeabf42eae2p-54, .hi: 0x1.371a7373aa9cbp0},
341 {.lo: -0x1.e958d3c9904bdp-54, .hi: 0x1.38cae6d05d866p0},
342 {.lo: -0x1.5e436d661f5e3p-56, .hi: 0x1.3a7db34e59ff7p0},
343 {.lo: -0x1.efff8375d29c3p-54, .hi: 0x1.3c32dc313a8e5p0},
344 {.lo: 0x1.ada0911f09ebcp-55, .hi: 0x1.3dea64c123422p0},
345 {.lo: -0x1.7d023f956f9f3p-54, .hi: 0x1.3fa4504ac801cp0},
346 {.lo: -0x1.ef3691c309278p-58, .hi: 0x1.4160a21f72e2ap0},
347 {.lo: -0x1.1c7dde35f7999p-55, .hi: 0x1.431f5d950a897p0},
348 {.lo: 0x1.89b7a04ef80dp-59, .hi: 0x1.44e086061892dp0},
349 {.lo: 0x1.c944bd1648a76p-54, .hi: 0x1.46a41ed1d0057p0},
350 {.lo: 0x1.3c1a3b69062fp-56, .hi: 0x1.486a2b5c13cdp0},
351 {.lo: 0x1.9cb62f3d1be56p-54, .hi: 0x1.4a32af0d7d3dep0},
352 {.lo: 0x1.d4397afec42e2p-56, .hi: 0x1.4bfdad5362a27p0},
353 {.lo: 0x1.8ecdbbc6a7833p-54, .hi: 0x1.4dcb299fddd0dp0},
354 {.lo: -0x1.4b309d25957e3p-54, .hi: 0x1.4f9b2769d2ca7p0},
355 {.lo: -0x1.f768569bd93efp-55, .hi: 0x1.516daa2cf6642p0},
356 {.lo: -0x1.07abe1db13cadp-55, .hi: 0x1.5342b569d4f82p0},
357 {.lo: -0x1.d689cefede59bp-55, .hi: 0x1.551a4ca5d920fp0},
358 {.lo: 0x1.9bb2c011d93adp-54, .hi: 0x1.56f4736b527dap0},
359 {.lo: 0x1.295e15b9a1de8p-55, .hi: 0x1.58d12d497c7fdp0},
360 {.lo: 0x1.6324c054647adp-54, .hi: 0x1.5ab07dd485429p0},
361 {.lo: 0x1.c4b1b816986a2p-60, .hi: 0x1.5c9268a5946b7p0},
362 {.lo: 0x1.ba6f93080e65ep-54, .hi: 0x1.5e76f15ad2148p0},
363 {.lo: -0x1.3e2429b56de47p-54, .hi: 0x1.605e1b976dc09p0},
364 {.lo: -0x1.383c17e40b497p-54, .hi: 0x1.6247eb03a5585p0},
365 {.lo: -0x1.c483c759d8933p-55, .hi: 0x1.6434634ccc32p0},
366 {.lo: -0x1.bb60987591c34p-54, .hi: 0x1.6623882552225p0},
367 {.lo: 0x1.038ae44f73e65p-57, .hi: 0x1.68155d44ca973p0},
368 {.lo: -0x1.bdd3413b26456p-54, .hi: 0x1.6a09e667f3bcdp0},
369 {.lo: -0x1.2895667ff0b0dp-56, .hi: 0x1.6c012750bdabfp0},
370 {.lo: -0x1.bbe3a683c88abp-57, .hi: 0x1.6dfb23c651a2fp0},
371 {.lo: -0x1.83c0f25860ef6p-55, .hi: 0x1.6ff7df9519484p0},
372 {.lo: -0x1.16e4786887a99p-55, .hi: 0x1.71f75e8ec5f74p0},
373 {.lo: -0x1.0a8d96c65d53cp-54, .hi: 0x1.73f9a48a58174p0},
374 {.lo: -0x1.0245957316dd3p-54, .hi: 0x1.75feb564267c9p0},
375 {.lo: 0x1.866b80a02162dp-54, .hi: 0x1.780694fde5d3fp0},
376 {.lo: -0x1.41577ee04992fp-55, .hi: 0x1.7a11473eb0187p0},
377 {.lo: 0x1.f124cd1164dd6p-54, .hi: 0x1.7c1ed0130c132p0},
378 {.lo: 0x1.05d02ba15797ep-56, .hi: 0x1.7e2f336cf4e62p0},
379 {.lo: -0x1.27c86626d972bp-54, .hi: 0x1.80427543e1a12p0},
380 {.lo: -0x1.d4c1dd41532d8p-54, .hi: 0x1.82589994cce13p0},
381 {.lo: -0x1.8d684a341cdfbp-55, .hi: 0x1.8471a4623c7adp0},
382 {.lo: -0x1.fc6f89bd4f6bap-54, .hi: 0x1.868d99b4492edp0},
383 {.lo: 0x1.994c2f37cb53ap-54, .hi: 0x1.88ac7d98a6699p0},
384 {.lo: 0x1.6e9f156864b27p-54, .hi: 0x1.8ace5422aa0dbp0},
385 {.lo: -0x1.0d55e32e9e3aap-56, .hi: 0x1.8cf3216b5448cp0},
386 {.lo: 0x1.5cc13a2e3976cp-55, .hi: 0x1.8f1ae99157736p0},
387 {.lo: -0x1.dd6792e582524p-54, .hi: 0x1.9145b0b91ffc6p0},
388 {.lo: -0x1.75fc781b57ebcp-57, .hi: 0x1.93737b0cdc5e5p0},
389 {.lo: -0x1.64b7c96a5f039p-56, .hi: 0x1.95a44cbc8520fp0},
390 {.lo: -0x1.d185b7c1b85d1p-54, .hi: 0x1.97d829fde4e5p0},
391 {.lo: -0x1.173bd91cee632p-54, .hi: 0x1.9a0f170ca07bap0},
392 {.lo: 0x1.c7c46b071f2bep-56, .hi: 0x1.9c49182a3f09p0},
393 {.lo: 0x1.824ca78e64c6ep-56, .hi: 0x1.9e86319e32323p0},
394 {.lo: -0x1.359495d1cd533p-54, .hi: 0x1.a0c667b5de565p0},
395 {.lo: 0x1.6305c7ddc36abp-54, .hi: 0x1.a309bec4a2d33p0},
396 {.lo: -0x1.d2f6edb8d41e1p-54, .hi: 0x1.a5503b23e255dp0},
397 {.lo: 0x1.bcb7ecac563c7p-54, .hi: 0x1.a799e1330b358p0},
398 {.lo: 0x1.0fac90ef7fd31p-54, .hi: 0x1.a9e6b5579fdbfp0},
399 {.lo: -0x1.f9234cae76cdp-55, .hi: 0x1.ac36bbfd3f37ap0},
400 {.lo: 0x1.7a1cd345dcc81p-54, .hi: 0x1.ae89f995ad3adp0},
401 {.lo: -0x1.bdef54c80e425p-54, .hi: 0x1.b0e07298db666p0},
402 {.lo: -0x1.2805e3084d708p-57, .hi: 0x1.b33a2b84f15fbp0},
403 {.lo: -0x1.c71dfbbba6de3p-54, .hi: 0x1.b59728de5593ap0},
404 {.lo: -0x1.5584f7e54ac3bp-56, .hi: 0x1.b7f76f2fb5e47p0},
405 {.lo: -0x1.efcd30e54292ep-54, .hi: 0x1.ba5b030a1064ap0},
406 {.lo: 0x1.23dd07a2d9e84p-55, .hi: 0x1.bcc1e904bc1d2p0},
407 {.lo: -0x1.efdca3f6b9c73p-54, .hi: 0x1.bf2c25bd71e09p0},
408 {.lo: 0x1.11065895048ddp-55, .hi: 0x1.c199bdd85529cp0},
409 {.lo: 0x1.b4537e083c60ap-54, .hi: 0x1.c40ab5fffd07ap0},
410 {.lo: 0x1.2884dff483cadp-54, .hi: 0x1.c67f12e57d14bp0},
411 {.lo: 0x1.1acbc48805c44p-56, .hi: 0x1.c8f6d9406e7b5p0},
412 {.lo: 0x1.503cbd1e949dbp-56, .hi: 0x1.cb720dcef9069p0},
413 {.lo: -0x1.dd83b53829d72p-55, .hi: 0x1.cdf0b555dc3fap0},
414 {.lo: -0x1.cbc3743797a9cp-54, .hi: 0x1.d072d4a07897cp0},
415 {.lo: -0x1.d487b719d8578p-54, .hi: 0x1.d2f87080d89f2p0},
416 {.lo: 0x1.2ed02d75b3707p-55, .hi: 0x1.d5818dcfba487p0},
417 {.lo: -0x1.11ec18beddfe8p-54, .hi: 0x1.d80e316c98398p0},
418 {.lo: 0x1.c2300696db532p-54, .hi: 0x1.da9e603db3285p0},
419 {.lo: 0x1.2da5778f018c3p-54, .hi: 0x1.dd321f301b46p0},
420 {.lo: -0x1.1a5cd4f184b5cp-54, .hi: 0x1.dfc97337b9b5fp0},
421 {.lo: -0x1.7b627817a1496p-54, .hi: 0x1.e264614f5a129p0},
422 {.lo: 0x1.39e8980a9cc8fp-55, .hi: 0x1.e502ee78b3ff6p0},
423 {.lo: 0x1.2d522ca0c8de2p-54, .hi: 0x1.e7a51fbc74c83p0},
424 {.lo: -0x1.e9c23179c2893p-54, .hi: 0x1.ea4afa2a490dap0},
425 {.lo: -0x1.c93f3b411ad8cp-54, .hi: 0x1.ecf482d8e67f1p0},
426 {.lo: 0x1.dc7f486a4b6bp-54, .hi: 0x1.efa1bee615a27p0},
427 {.lo: 0x1.3a1a5bf0d8e43p-54, .hi: 0x1.f252b376bba97p0},
428 {.lo: 0x1.9d3e12dd8a18bp-54, .hi: 0x1.f50765b6e454p0},
429 {.lo: -0x1.dbb12d006350ap-54, .hi: 0x1.f7bfdad9cbe14p0},
430 {.lo: 0x1.74853f3a5931ep-55, .hi: 0x1.fa7c1819e90d8p0},
431 {.lo: 0x1.2eb74966579e7p-57, .hi: 0x1.fd3c22b8f71f1p0},
432};
433
434LIBC_NOINLINE LIBC_INLINE double
435exp_denorm(double x, int hi, double r_hi, double lo_u, double err_d,
436 const ExpBase192 *log2_base = nullptr) {
437 using FPBits = typename fputil::FPBits<double>;
438 // Exact subnormal power of 2 when x is an integer in [-1074, -1022].
439 if (log2_base == nullptr && LIBC_UNLIKELY(static_cast<double>(hi) == x))
440 return FPBits(uint64_t(1) << (hi + 1074)).get_val();
441
442 if (auto r_denorm = ziv_test_denorm(hi, mid: r_hi, lo: lo_u - err_d, err: err_d);
443 LIBC_LIKELY(r_denorm.has_value())) {
444 double res = r_denorm.value();
445 if (LIBC_UNLIKELY(FPBits(res).is_normal()))
446 return res;
447
448 fputil::set_errno_if_required(ERANGE);
449 fputil::raise_underflow_except_if_required<double>();
450 return res;
451 }
452 return exp_accurate_128(x, log2_base);
453}
454
455LIBC_ALWAYS_INLINE double exp2_eval(double x, DoubleDouble u, int hi,
456 unsigned idx, double err_d = EXP2_ERR_D,
457 const ExpBase192 *log2_base = nullptr) {
458 DoubleDouble exp2_mid = EXP2_MID[idx];
459
460 // Reconstruction of 2^(idx / 2^7) * exp(u) ~ exp2_mid * (1 + u + u^2 * P(u))
461 // where:
462 // exp2_mid = 2^(idx / 2^7) = exp2_mid.hi + exp2_mid.lo
463 // u = dx * log(b) = u.hi + u.lo
464#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
465 // Decompose exp2_mid * (1 + u + u^2 * P(u)) ~ r_hi + (r_lo + tail):
466 // r_hi + r_lo = exp2_mid.hi * (1 + u.hi)
467 // tail = exp2_mid.hi * (u.lo + u^2 * P(u))
468 // + exp2_mid.lo * (1 + u.hi)
469 // where the omitted cross term exp2_mid.lo * (u.lo + u^2 * P(u)) is bounded
470 // by 2^(-54 - 18) = 2^-72 and absorbed into err_d.
471
472 // p ~ (exp(u) - 1 - u) / u^2
473 double u2 = u.hi * u.hi;
474 double p = exp_poly(u.hi, u2);
475
476 // r_hi = round(exp2_mid.hi * (1 + u.hi))
477 double r_hi = fputil::multiply_add(exp2_mid.hi, u.hi, exp2_mid.hi);
478 // mid_lo ~ exp2_mid.lo * (1 + u.hi)
479 double mid_lo = fputil::multiply_add(exp2_mid.lo, u.hi, exp2_mid.lo);
480 // poly_lo ~ u.lo + u^2 * P(u)
481 double poly_lo = fputil::multiply_add(u2, p, u.lo);
482
483 // Fold +/-err_d into mid_lo so tail_u, tail_l, and r_lo can be evaluated in
484 // parallel:
485 // mid_u = mid_lo + err_d
486 // mid_l = mid_lo - err_d
487 double mid_u = mid_lo + err_d;
488 double mid_l = mid_lo - err_d;
489 // tail_{u,l} ~ exp2_mid.hi * (u.lo + u^2 * P(u)) + exp2_mid.lo * (1 + u.hi)
490 // +/- err_d
491 double tail_u = fputil::multiply_add(exp2_mid.hi, poly_lo, mid_u);
492 double tail_l = fputil::multiply_add(exp2_mid.hi, poly_lo, mid_l);
493 // r_lo = exp2_mid.hi * (1 + u.hi) - r_hi (exact to within 2^-105)
494 double r_lo = fputil::multiply_add(exp2_mid.hi, u.hi, exp2_mid.hi - r_hi);
495
496 // lo_{u,l} ~ r_lo + tail +/- err_d
497 double lo_u = r_lo + tail_u;
498 double lo_l = r_lo + tail_l;
499#else
500 // p ~ (exp(u) - 1 - u) / u^2
501 double u_sum = u.hi + u.lo;
502 double u2 = u_sum * u_sum;
503 double p = exp_poly(u: u_sum, u2);
504
505 // Split exp2_mid = m_hi + m_lo with lsb(m_hi) >= 2^-22 and |m_lo| <= 2^-23:
506 double m_hi = (exp2_mid.hi + SPLIT_22) - SPLIT_22; // lsb = 2^-22
507 double m_lo = (exp2_mid.hi - m_hi) + exp2_mid.lo;
508
509 // Since lsb(u.hi) >= 2^-30 and r_hi is in [0.997, 2.0), both m_hi * u.hi and
510 // r_hi = m_hi * (1 + u.hi) are exact:
511 double prod_hi = m_hi * u.hi; // exact, lsb = 2^-52
512 double r_hi = m_hi + prod_hi; // exact
513
514 // poly_lo ~ u.lo + u^2 * P(u)
515 double u2_p = u2 * p;
516 double poly_lo = u.lo + u2_p;
517
518 // r_lo ~ exp2_mid.hi * (u.lo + u^2 * P(u)) + m_lo * (1 + u.hi)
519 double l0 = exp2_mid.hi * poly_lo;
520 double l1 = m_lo * (1.0 + u.hi);
521 double r_lo = l0 + l1;
522
523 // lo_{u,l} = r_lo +/- err_d
524 double lo_u = r_lo + err_d;
525 double lo_l = r_lo - err_d;
526#endif
527 // Ziv's rounding test in the current rounding mode:
528 double upper = r_hi + lo_u;
529 double lower = r_hi + lo_l;
530
531 if (LIBC_LIKELY(upper == lower && hi > -1022))
532 return scale_exp2(r: upper, hi);
533
534 if (LIBC_UNLIKELY(hi <= -1022))
535 return exp_denorm(x, hi, r_hi, lo_u, err_d, log2_base);
536
537 return exp_accurate_128(x, log2_base);
538}
539
540#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
541
542} // namespace double_eval
543} // namespace math
544} // namespace LIBC_NAMESPACE_DECL
545
546#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
547