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/// Double-precision implementation of exp2.
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
16
17#include "src/__support/FPUtil/FPBits.h"
18#include "src/__support/FPUtil/double_double.h"
19#include "src/__support/FPUtil/multiply_add.h"
20#include "src/__support/FPUtil/nearest_integer.h"
21#include "src/__support/common.h"
22#include "src/__support/macros/attributes.h"
23#include "src/__support/macros/config.h"
24#include "src/__support/macros/optimization.h"
25#include "src/__support/macros/properties/cpu_features.h"
26#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
27#include "src/__support/math/exp_double_utils.h"
28#else
29#include "src/__support/math/exp_utils.h"
30#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
31
32namespace LIBC_NAMESPACE_DECL {
33namespace math {
34namespace double_eval {
35
36LIBC_ALWAYS_INLINE double exp2(double x) {
37 using FPBits = typename fputil::FPBits<double>;
38 FPBits xbits(x);
39
40 uint64_t x_u = xbits.uintval();
41 uint64_t x_abs = xbits.abs().uintval();
42
43 // Bound constants (with biased exponent and mantissa bits):
44 // - Tiny bound: nextafter(2^-54, +inf) = 0x1.0000000000001p-54
45 // (since |x| <= 2^-54 is handled by 1.0 + x in set_exceptional)
46 constexpr uint64_t TINY_BOUND = 0x3c90'0000'0000'0001ULL;
47 // - Overflow bound: 1024.0
48 constexpr uint64_t OVERFLOW_BOUND = 0x4090'0000'0000'0000ULL;
49 // -1024.0
50 constexpr uint64_t NEG_OVERFLOW_BOUND = 0xc090'0000'0000'0000ULL;
51 // - Underflow bound: -1075.0
52 constexpr uint64_t UNDERFLOW_BOUND = 0xc090'cc00'0000'0000ULL;
53
54 // Normal range: 2^-54 <= |x| < 1024.
55 if (LIBC_UNLIKELY(x_abs - TINY_BOUND >= OVERFLOW_BOUND - TINY_BOUND)) {
56 // Exclude non-underflowing negative inputs in (-1075, -1024].
57 if (LIBC_LIKELY(x_u < NEG_OVERFLOW_BOUND || x_u >= UNDERFLOW_BOUND))
58 return set_exceptional(x);
59 }
60
61 // Now -1075 < x < -2^-54 or 2^-54 < x < 1024.
62
63 // Range reduction:
64 // Let x = k * 2^-N + dx
65 // in which:
66 // k = round(x * 2^N)
67 // hi = k >> N
68 // idx = k & (2^N - 1)
69 // then:
70 // 2^x = 2^hi * 2^(idx / 2^N) * 2^dx.
71
72 double kd = fputil::nearest_integer(x: x * EXP_SCALE);
73 int k = static_cast<int>(kd);
74
75 int hi = k >> EXP_N_BITS;
76 unsigned idx = static_cast<unsigned>(k & EXP_MASK);
77
78 double dx = fputil::multiply_add(x: kd, y: -EXP_SCALE_INV, z: x); // exact
79
80#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
81 double u = dx * LOG_2.hi;
82 return exp2_eval(x, u, hi, idx, /*is_exp2=*/true);
83#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
84#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
85 DoubleDouble u = fputil::exact_mult(dx, LOG_2.hi);
86 u.lo = fputil::multiply_add(dx, LOG_2.lo, u.lo);
87#else // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
88 // Without FMA, we use fixed-point splitting to avoid Veltkamp's splitting:
89 // Let:
90 // m = 2^(idx / 2^7) = m_hi + m_lo
91 // u = dx * log(2) = u_hi + u_lo
92 //
93 // In exp2_eval, we approximate 2^(idx / 2^7) * 2^dx = m * exp(u) by:
94 // m * (1 + u + u^2 * P(u)) ~ r_hi + r_lo
95 // where:
96 // r_hi = m_hi + m_hi * u_hi
97 // r_lo = m * (u_lo + u^2 * P(u)) + m_lo * (1 + u_hi)
98 //
99 // Since r_hi is in [0.997, 2.0), for r_hi to be exact in double precision,
100 // we need:
101 // lsb(m_hi * u_hi) >= 2^-52.
102 // We choose:
103 // lsb(m_hi) >= 2^-22 => |m_lo| <= 2^-23
104 // lsb(u_hi) >= 2^-30
105 // so that lsb(m_hi * u_hi) >= 2^-52, and the rounding errors from
106 // m_lo * (1 + u_hi) are bounded by 2^(-23 - 52) = 2^-75.
107 //
108 // To get lsb(u_hi) >= 2^-30 from u = dx * log(2), since |dx| <= 2^-8 and
109 // log(2) < 1, we split dx at 2^-18 and log(2) at 2^-12:
110 // dx = dx_hi + dx_lo (|dx_lo| <= 2^-19)
111 // log(2) = LOG_2_HI_12 + LOG_2_LO_12 (|LOG_2_LO_12| < 2^-13)
112 // Then:
113 // u_hi = dx_hi * LOG_2_HI_12
114 // is an exact multiple of 2^-30, and the two remainder products in:
115 // u_lo = dx_lo * LOG_2_HI_12 + dx * LOG_2_LO_12
116 // are both bounded by 2^-19, keeping the rounding errors of u_lo below 2^-71.
117 double dx_hi = (dx + SPLIT_18) - SPLIT_18; // lsb = 2^-18
118 double dx_lo = dx - dx_hi;
119 double u_hi = dx_hi * LOG_2_HI_12; // exact, lsb = 2^-30
120 double t0 = dx_lo * LOG_2_HI_12;
121 double t1 = dx * LOG_2_LO_12;
122 double u_lo = t0 + t1;
123 DoubleDouble u = {.lo: u_lo, .hi: u_hi};
124#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
125 return exp2_eval(x, u, hi, idx);
126#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
127}
128
129} // namespace double_eval
130} // namespace math
131} // namespace LIBC_NAMESPACE_DECL
132
133#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
134