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/// 128-bit accurate path for double-precision exp2(x).
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_ACCURATE_128_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_ACCURATE_128_H
16
17#include "src/__support/FPUtil/FPBits.h"
18#include "src/__support/FPUtil/dyadic_float.h"
19#include "src/__support/common.h"
20#include "src/__support/frac128.h"
21#include "src/__support/macros/attributes.h"
22#include "src/__support/macros/config.h"
23#include "src/__support/macros/optimization.h"
24#include "src/__support/uint128.h"
25
26namespace LIBC_NAMESPACE_DECL {
27namespace math {
28namespace double_eval {
29
30using DFloat128 = typename fputil::DyadicFloat<128>;
31
32// Table of 2^(i / 16) with i = 0..15 in 128-bit fixed-point format:
33// Generated by Sollya with:
34// > display = hexadecimal!;
35// > prec = 256;
36// > for i from 0 to 15 do {
37// a = 2^(i / 16);
38// c_int = round(a * 2^127, 256, RN);
39// hi = floor(c_int * 2^(-64));
40// lo = c_int - hi * 2^64;
41// print(" Frac128({", lo, "ULL,", hi, "ULL}), // 2^(", i, "/16)");
42// };
43LIBC_INLINE_VAR constexpr Frac128 EXP2_MID_FRAC128[16] = {
44 // lsb = 2^-127
45 Frac128({0x0000000000000000ULL, 0x8000000000000000ULL}), // 2^(0/16)
46 Frac128({0xc5c95b8c2154c1b2ULL, 0x85aac367cc487b14ULL}), // 2^(1/16)
47 Frac128({0xfbe4628758a53c90ULL, 0x8b95c1e3ea8bd6e6ULL}), // 2^(2/16)
48 Frac128({0x0fd6d8e0ae5ac9d8ULL, 0x91c3d373ab11c336ULL}), // 2^(3/16)
49 Frac128({0x46ad23182e42f6f6ULL, 0x9837f0518db8a96fULL}), // 2^(4/16)
50 Frac128({0xa0911f09ebb9fdd1ULL, 0x9ef5326091a111adULL}), // 2^(5/16)
51 Frac128({0x1cbd7f621710701bULL, 0xa5fed6a9b15138eaULL}), // 2^(6/16)
52 Frac128({0x4980a8c8f59a2ec4ULL, 0xad583eea42a14ac6ULL}), // 2^(7/16)
53 Frac128({0x597d89b3754abe9fULL, 0xb504f333f9de6484ULL}), // 2^(8/16)
54 Frac128({0xa8811fb66d0faf7aULL, 0xbd08a39f580c36beULL}), // 2^(9/16)
55 Frac128({0x3e2ad0c964dd9f37ULL, 0xc5672a115506daddULL}), // 2^(10/16)
56 Frac128({0xe235838f95f2c6edULL, 0xce248c151f8480e3ULL}), // 2^(11/16)
57 Frac128({0x39a68bb9902d3fdeULL, 0xd744fccad69d6af4ULL}), // 2^(12/16)
58 Frac128({0x065895048dd333caULL, 0xe0ccdeec2a94e111ULL}), // 2^(13/16)
59 Frac128({0xd02d75b3706e54fbULL, 0xeac0c6e7dd24392eULL}), // 2^(14/16)
60 Frac128({0x7b9d0c7aed980fc3ULL, 0xf5257d152486cc2cULL}), // 2^(15/16)
61};
62
63// Polynomial approximation for (2^x - 1) / x:
64// For x in [0, 2^-4], we approximate:
65// P(x) ~ (2^x - 1) / x
66// then:
67// 2^x = 1 + x * P(x).
68//
69// Minimax polynomial coefficients generated by Sollya with:
70// > prec = 500;
71// > P = fpminimax((2^x - 1)/x, 14, [|128...|], [2^(-150), 2^(-4)], fixed);
72// > for i from 0 to 14 do {
73// c = coeff(P, i);
74// c_int = round(c * 2^128, 256, RN);
75// hi = floor(c_int * 2^(-64));
76// lo = c_int - hi * 2^64;
77// print("Frac128({", lo, "ULL,", hi, "ULL}), // c_", i);
78// };
79// > dirtyinfnorm(2^x - (1 + x * P), [0, 2^(-4)]);
80// 0x1.278a...p-139 < 2^-138.
81//
82// Expanding P(x) in Horner form:
83// P(x) = c_0 + x * (c_1 + x * (c_2 + ...))
84// where c_k ~ (log(2))^(k+1) / (k+1)! are all positive and in (0, 1).
85LIBC_INLINE_VAR constexpr Frac128 EXP2_POLY_128[15] = {
86 // lsb = 2^-128
87 Frac128({0xc9e3b39803f2f6afULL, 0xb17217f7d1cf79abULL}), // c_0
88 Frac128({0xde2d60dd92e6c1ceULL, 0x3d7f7bff058b1d50ULL}), // c_1
89 Frac128({0x99d3b15d9958cfa0ULL, 0x0e35846b82505fc5ULL}), // c_2
90 Frac128({0x39977c16aee20457ULL, 0x0276556df749cee5ULL}), // c_3
91 Frac128({0x41c5fda1c13b8fb1ULL, 0x005761ff9e299cc4ULL}), // c_4
92 Frac128({0xb7a58759e18c05e7ULL, 0x000a184897c363c3ULL}), // c_5
93 Frac128({0x3434f0e418799c8bULL, 0x0000ffe5fe2c4586ULL}), // c_6
94 Frac128({0x241c41f55c65febbULL, 0x0000162c0223a5c8ULL}), // c_7
95 Frac128({0x77e77364ac415ce2ULL, 0x000001b5253d395eULL}), // c_8
96 Frac128({0xfe79c174054850e2ULL, 0x0000001e4cf5158bULL}), // c_9
97 Frac128({0x29d719d4d61d89c9ULL, 0x00000001e8cac72dULL}), // c_10
98 Frac128({0xd15314bfa6b56886ULL, 0x000000001c3bd6b4ULL}), // c_11
99 Frac128({0xd8de7d5b7789dd3aULL, 0x0000000001815e5dULL}), // c_12
100 Frac128({0x491912c310d5531bULL, 0x0000000000132403ULL}), // c_13
101 Frac128({0x8d1d877caeba1cfcULL, 0x000000000000c220ULL}), // c_14
102};
103
104struct ExpBase192 {
105 int base_exp;
106 UInt<192> log2_b; // lsb = 2^(base_exp - 191)
107};
108
109// Evaluate 2^(idx * 2^-4 + u) in [1, 2) as Frac128 (lsb = 2^-127) for
110// idx in [0, 15] and u in [0, 2^-4) (lsb = 2^-128).
111LIBC_INLINE Frac128 exp2_eval_frac128(unsigned idx, const Frac128 &u) {
112 // p ~ (2^u - 1) / u
113 Frac128 p = EXP2_POLY_128[14]; // lsb = 2^-128
114 LIBC_LOOP_NOUNROLL
115 for (int i = 13; i >= 0; --i)
116 p = p * u + EXP2_POLY_128[i]; // lsb = 2^-128
117
118 // Reconstruction:
119 // u * p ~ 2^u - 1 (lsb = 2^-128)
120 // m = 2^(idx * 2^-4) (lsb = 2^-127)
121 // m + m * (u * p) ~ m * 2^u (lsb = 2^-127)
122 Frac128 m = EXP2_MID_FRAC128[idx]; // lsb = 2^-127
123 return m + m * (u * p); // lsb = 2^-127
124}
125
126// Compute 2^z in 128-bit precision for a DyadicFloat<128> z with
127// 2^-55 < |z| < 2^11 (z.exponent in [-181, -117]) and optional extra low limb
128// lo (lsb = 2^(z.exponent - 64)):
129// z = k * 2^-4 + u, 0 <= u < 2^-4
130// hi = k >> 4
131// idx = k & 0x0f
132// 2^z = 2^hi * EXP2_MID_FRAC128[idx] * (1 + u * P(u)).
133LIBC_INLINE DFloat128 exp2_df128(const DFloat128 &z, uint64_t lo = 0) {
134 // z.mantissa has lsb = 2^z.exponent, so the 2^-4 bit of |z| is at bit
135 // (-4 - z.exponent) of z.mantissa, leaving
136 // shift = 128 - (-4 - z.exponent) = z.exponent + 132 bits
137 // at or above 2^-4 in z.mantissa.val[1].
138 // After range reduction, 2^-55 < |z| < 2^11, so z.exponent is in
139 // [-181, -117] and shift is in [-49, 15].
140 int shift = z.exponent + 132;
141 int k = 0;
142 UInt<128> frac;
143 if (LIBC_LIKELY(shift > 0)) {
144 // |z| >= 2^-4: split |z| = k * 2^-4 + frac * 2^-132.
145 // Here shift is in [1, 15] and shr = 64 - shift is in [49, 63], so all
146 // 64-bit shifts below are strictly in [1, 63] (no undefined behavior).
147 int shr = 64 - shift; // in [49, 63]
148 // k = floor(|z| * 2^4)
149 k = static_cast<int>(z.mantissa.val[1] >> shr); // lsb = 2^-4
150 // Shift z.mantissa left by shift and fill the low shift bits from lo:
151 frac = UInt<128>({(z.mantissa.val[0] << shift) | (lo >> shr),
152 (z.mantissa.val[1] << shift) |
153 (z.mantissa.val[0] >> shr)}); // lsb = 2^-132
154 } else {
155 // |z| < 2^-4: k = 0 and frac = |z| * 2^132, with -shift in [0, 49].
156 frac = (-shift < 128) ? (z.mantissa >> (-shift))
157 : UInt<128>(0); // lsb = 2^-132
158 }
159
160 if (z.sign.is_neg()) {
161 if (LIBC_UNLIKELY(frac.is_zero())) {
162 k = -k;
163 } else {
164 // One's-complement negation on frac (differs from two's-complement by
165 // 2^-132, which is 1/16 ULP of u):
166 // -|z| ~ (-k - 1) * 2^-4 + (2^-4 - (frac + 1) * 2^-132).
167 frac = ~frac;
168 k = ~k;
169 }
170 }
171
172 // Decompose k * 2^-4 = hi + idx * 2^-4 and scale frac from lsb = 2^-132 to
173 // u in [0, 2^-4) with lsb = 2^-128:
174 int hi = k >> 4; // lsb = 1
175 unsigned idx = static_cast<unsigned>(k & 0x0f); // lsb = 2^-4
176 Frac128 u((frac >> 4).val); // lsb = 2^-128
177 Frac128 m_final = exp2_eval_frac128(idx, u); // lsb = 2^-127
178
179 // Since m_final is in [1, 2) with lsb = 2^-127, its leading bit is already
180 // at bit 127, so r is normalized with exponent = hi - 127:
181 DFloat128 r;
182 r.sign = Sign::POS;
183 r.exponent = hi - 127;
184 r.mantissa = UInt<128>(m_final.val);
185 return r;
186}
187
188// Compute b^x (for b = 2, e, or 10) in 128-bit precision using direct base-2
189// 192-bit fixed-point range reduction. Here log2_base->log2_b holds
190// round(log2(b) * 2^(191 - base_exp)), or log2_base is nullptr when b = 2.
191LIBC_NOINLINE LIBC_INLINE double
192exp_accurate_128(double x, const ExpBase192 *log2_base = nullptr) {
193 using FPBits = fputil::FPBits<double>;
194 FPBits xbits(x);
195 int exp_x = xbits.get_exponent();
196
197 // Normalize |x|'s 53-bit mantissa to [2^63, 2^64):
198 uint64_t mant_x = (xbits.get_mantissa() | (1ULL << 52)) << 11;
199
200 DFloat128 z;
201 z.sign = xbits.sign();
202 uint64_t w0 = 0;
203
204 if (log2_base != nullptr) {
205 // 1-limb x 3-limb product:
206 // mant_x has lsb = 2^(exp_x - 63)
207 // log2_base->log2_b has lsb = 2^(base_exp - 191)
208 // p = mant_x * log2_base->log2_b has lsb = 2^(exp_x + base_exp - 254)
209 UInt<256> p = log2_base->log2_b.ful_mul(other: UInt<64>(mant_x));
210 // Extract the top 192 bits (p.val[3], p.val[2], p.val[1]) as
211 // (z.mantissa, w0):
212 z.exponent = exp_x + log2_base->base_exp - 126;
213 z.mantissa = UInt<128>({p.val[2], p.val[3]}); // lsb = 2^z.exponent
214 w0 = p.val[1]; // lsb = 2^(z.exponent - 64)
215 } else {
216 // For b = 2, log2(b) = 1, so z = x directly:
217 z.exponent = exp_x - 127;
218 z.mantissa = UInt<128>({0, mant_x}); // lsb = 2^(exp_x - 127)
219 }
220
221 return static_cast<double>(exp2_df128(z, lo: w0));
222}
223
224} // namespace double_eval
225} // namespace math
226} // namespace LIBC_NAMESPACE_DECL
227
228#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_ACCURATE_128_H
229