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 pow(x, y).
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_128_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_128_H
16
17#include "src/__support/CPP/bit.h"
18#include "src/__support/FPUtil/FEnvImpl.h"
19#include "src/__support/FPUtil/FPBits.h"
20#include "src/__support/FPUtil/dyadic_float.h"
21#include "src/__support/FPUtil/multiply_add.h"
22#include "src/__support/FPUtil/nearest_integer.h"
23#include "src/__support/FPUtil/rounding_mode.h"
24#include "src/__support/common.h"
25#include "src/__support/frac128.h"
26#include "src/__support/frac64.h"
27#include "src/__support/integer_literals.h"
28#include "src/__support/macros/config.h"
29#include "src/__support/macros/optimization.h"
30#include "src/__support/math/common_constants.h"
31#include "src/__support/math/pow_accurate_256.h"
32#include "src/__support/math/pow_utils.h"
33#include "src/__support/uint128.h"
34
35namespace LIBC_NAMESPACE_DECL {
36namespace math {
37namespace pow_internal {
38
39using DFloat128 = typename fputil::DyadicFloat<128>;
40using MantissaType = typename DFloat128::MantissaType;
41using LIBC_NAMESPACE::operator""_u128;
42
43// exp2_f128 shares EXP2_MID_FRAC256 from pow_accurate_256.h via
44// Frac256::to_frac128().
45
46// Polynomial approximation for log2(1 + dx):
47// For dx in [-2^-8, 2^-7], we approximate:
48// P(x) ~ log2(1 + x) / (2 * x)
49// then:
50// log2(1 + dx) = 2 * dx * P(dx).
51//
52// Minimax polynomial coefficients generated by Sollya with:
53// > prec = 256;
54// > P = fpminimax(log2(1 + x) / (2 * x), 15, [|127...|], [-2^(-8), 2^(-7)],
55// fixed);
56// > for i from 0 to 15 do {
57// c = coeff(P, i);
58// c_int = round(abs(c) * 2^127, 256, RN);
59// lo = c_int mod 2^64;
60// hi = (c_int - lo) / 2^64;
61// print("Frac128({", lo, "ULL,", hi, "ULL}), // a_", i);
62// };
63// > dirtyinfnorm(log2(1 + x) - 2 * x * P, [-2^(-8), 2^(-7)]);
64// 0x1.c6524...p-139 < 2^-138.
65//
66// Expanding P(x) in Horner form:
67// P(x) = a_0 - a_1 * x + a_2 * x^2 - a_3 * x^3 + ...
68// We store the absolute values of the coefficients |a_k| as Frac128.
69//
70// With y = |dx|,
71// - if dx >= 0:
72// P(y) = |a_0| - y * (|a_1| - y * (|a_2| - ...))
73// - if dx < 0:
74// P(-y) = |a_0| + y * (|a_1| + y * (|a_2| + ...))
75LIBC_INLINE_VAR constexpr Frac128 LOG2_POLY_128[16] = {
76 Frac128({0xdf43ff68348e9f44ULL, 0x5c551d94ae0bf85dULL}), // a_0
77 Frac128({0xefa1ffb41a474fa2ULL, 0x2e2a8eca5705fc2eULL}), // a_1
78 Frac128({0x9fc15522bc34caf6ULL, 0x1ec709dc3a03fd74ULL}), // a_2
79 Frac128({0x77d0ffda0d2838eaULL, 0x171547652b82fe17ULL}), // a_3
80 Frac128({0x2ca733033277ada0ULL, 0x12776c50ef9bfe79ULL}), // a_4
81 Frac128({0x4fe0aa8003f351fbULL, 0x0f6384ee1d01febaULL}), // a_5
82 Frac128({0xb2490fe2b526029dULL, 0x0d30bb153d6f6c9fULL}), // a_6
83 Frac128({0xbc03c43080705d70ULL, 0x0b8aa3b295c17f0bULL}), // a_7
84 Frac128({0x8db7017adcd9fd4dULL, 0x0a42589ebe015466ULL}), // a_8
85 Frac128({0x7a1820b5f15e8892ULL, 0x093bb62877cdff26ULL}), // a_9
86 Frac128({0xfb43ea327ddb209cULL, 0x0864d424ca0b0f9dULL}), // a_10
87 Frac128({0xca32556a02086730ULL, 0x07b1c2770e8b0fabULL}), // a_11
88 Frac128({0xb05fb6d59979dddeULL, 0x071a3d579502ad57ULL}), // a_12
89 Frac128({0x4abbb99c78464ef1ULL, 0x06985d880590b72bULL}), // a_13
90 Frac128({0x86477a0bb2fb6379ULL, 0x06282cce144a837eULL}), // a_14
91 Frac128({0x25363223e4195930ULL, 0x05c5acde28e4f0bfULL}), // a_15
92};
93
94// Accurate log2(x) in 128-bit precision reusing range reduction from fast pass:
95// x = 2^x_e * m_x, with 1 <= m_x < 2
96// r = RD[idx_x]
97// dx = r * m_x - 1, with -2^-8 <= dx < 2^-7
98// log2(x) = x_e + (-log2(r)) + log2(1 + dx).
99// The value -log2(r) is provided by LOG2_RD_FRAC256[idx_x].
100// log2(1 + dx) is evaluated with a degree-15 minimax polynomial using a 2-tier
101// Horner scheme (Frac64 and Frac128).
102LIBC_INLINE DFloat128 log2_f128(int x_e, unsigned idx_x, double dx) {
103 using FPBits = fputil::FPBits<double>;
104
105 if (dx == 0.0) {
106 Frac128 log2_m = LOG2_RD_FRAC256[idx_x].to_frac128();
107 DFloat128 m_df(Sign::POS, -127, log2_m);
108 return fputil::quick_add(a: DFloat128(static_cast<double>(x_e)), b: m_df);
109 }
110
111 double abs_dx = (dx < 0.0) ? -dx : dx;
112 FPBits bits(abs_dx);
113 // Append hidden bit.
114 uint64_t mant = bits.get_mantissa() | (1ULL << 52);
115 int shift = 127 + bits.get_exponent() - 52;
116 Frac128 y =
117 (shift >= 0) ? Frac128((UInt<128>(mant) << shift).val) : Frac128(0);
118
119 // Evaluate log2(1 + dx) using Horner scheme in 2 stages:
120 // - Degree 8-15: 64-bit precision evaluation.
121 // - Degree 0-7: 128-bit precision evaluation.
122 //
123 // Since y = |dx| <= 2^-7, truncation errors at degree k is bounded by:
124 // y^k <= 2^(-7k).
125 //
126 // Step 1: 64-bit evaluation, with truncation error at degree 8 is bounded by:
127 // 2^-63 * y^8 <= 2^-119.
128 Frac64 y64 = y.to_frac64();
129 Frac64 p64 = LOG2_POLY_128[15].to_frac64();
130 if (dx >= 0.0) {
131 for (int k = 14; k >= 8; --k)
132 p64 = LOG2_POLY_128[k].to_frac64() - ((p64 * y64) << 1);
133 } else {
134 for (int k = 14; k >= 8; --k)
135 p64 = LOG2_POLY_128[k].to_frac64() + ((p64 * y64) << 1);
136 }
137
138 // Step 2: 128-bit evaluation.
139 Frac128 p({0, p64.val[0]});
140 if (dx >= 0.0) {
141 for (int k = 7; k >= 0; --k)
142 p = LOG2_POLY_128[k] - ((p * y) << 1);
143 } else {
144 for (int k = 7; k >= 0; --k)
145 p = LOG2_POLY_128[k] + ((p * y) << 1);
146 }
147
148 // log2(1 + dx) = 2 * y * P(y):
149 Frac128 log2_1p = (y * p) << 2;
150 Frac128 log2_rd = LOG2_RD_FRAC256[idx_x].to_frac128();
151 Frac128 log2_m = (dx >= 0.0) ? (log2_rd + log2_1p) : (log2_rd - log2_1p);
152
153 DFloat128 m_df(Sign::POS, -127, log2_m);
154 return fputil::quick_add(a: DFloat128(static_cast<double>(x_e)), b: m_df);
155}
156
157// Polynomial approximation for (2^x - 1) / x:
158// For x in [-2^-8, 2^-7], we approximate:
159// P(x) ~ (2^x - 1) / x
160// then:
161// 2^x = 1 + x * P(x).
162//
163// Minimax polynomial coefficients generated by Sollya with:
164// > prec = 500;
165// > P = fpminimax((2^x - 1)/x, 11, [|127...|], [-2^(-8), 2^(-7)], fixed);
166// > for i from 0 to 11 do {
167// c = coeff(P, i);
168// c_int = round(abs(c) * 2^127, 256, RN);
169// lo = c_int mod 2^64;
170// hi = (c_int - lo) / 2^64;
171// print("Frac128({", lo, "ULL,", hi, "ULL}), // c_", i);
172// };
173// > dirtyinfnorm(2^x - (1 + x * P), [-2^(-8), 2^(-7)]);
174// 0x1.4c648...p-139 < 2^-138.
175//
176// Expanding P(x) in Horner form:
177// P(x) = c_0 + c_1 * x + c_2 * x^2 + ...
178// where c_k ~ (log(2))^(k+1) / (k+1)! are all positive.
179//
180// With u = |x|,
181// - if x >= 0:
182// P(u) = c_0 + u * (c_1 + u * (c_2 + ...))
183// - if x < 0:
184// P(-u) = c_0 - u * (c_1 - u * (c_2 - ...))
185LIBC_INLINE_VAR constexpr Frac128 EXP2_POLY_128[12] = {
186 Frac128({0xe4f1d9cc01f97b58ULL, 0x58b90bfbe8e7bcd5ULL}), // c_0
187 Frac128({0x6f16b06ec9735fcbULL, 0x1ebfbdff82c58ea8ULL}), // c_1
188 Frac128({0xcce9d8aeccad816eULL, 0x071ac235c1282fe2ULL}), // c_2
189 Frac128({0x9ccbbe0b53ecfa40ULL, 0x013b2ab6fba4e772ULL}), // c_3
190 Frac128({0x20e2fed5a256d4e6ULL, 0x002bb0ffcf14ce62ULL}), // c_4
191 Frac128({0xdbd2c2a4af40b4b6ULL, 0x00050c244be1b1e1ULL}), // c_5
192 Frac128({0x1a199575c57ddfaaULL, 0x00007ff2ff1622c3ULL}), // c_6
193 Frac128({0x11fd7e22defbeaf5ULL, 0x00000b160111d2e4ULL}), // c_7
194 Frac128({0x7681395eb574dbfbULL, 0x000000da929e9cafULL}), // c_8
195 Frac128({0x1814285a9da99308ULL, 0x0000000f267a8ac6ULL}), // c_9
196 Frac128({0xafddba4f901557d4ULL, 0x00000000f46563d0ULL}), // c_10
197 Frac128({0x7ae715abc468c5f0ULL, 0x000000000e1de423ULL}), // c_11
198};
199
200// Compute 2^z in 128-bit precision:
201// Range reduction:
202// k = round(z * 64)
203// hi = k >> 6
204// idx = k & 0x3f
205// lo = z - k * 2^-6, with |lo| <= 2^-7.
206// Then:
207// 2^z = 2^hi * EXP2_MID_FRAC256[idx] * 2^lo
208// = 2^hi * EXP2_MID_FRAC256[idx] * (1 + lo * P(lo)).
209LIBC_INLINE DFloat128 exp2_f128(const DFloat128 &z) {
210 double z_d = static_cast<double>(z);
211 double z_scaled = z_d * 64.0;
212 double kd = fputil::nearest_integer(x: z_scaled);
213 int k = static_cast<int>(kd);
214
215 int hi = k >> 6;
216 unsigned idx = static_cast<unsigned>(k & 0x3f);
217
218 DFloat128 kd_f128(kd * 0x1.0p-6);
219 DFloat128 lo = fputil::quick_add(a: z, b: -kd_f128);
220
221 Frac128 m = EXP2_MID_FRAC256[idx].to_frac128();
222 if (LIBC_UNLIKELY(lo.mantissa.is_zero()))
223 return DFloat128(Sign::POS, hi - 127, m);
224
225 bool lo_is_neg = (lo.sign == Sign::NEG);
226 int shift = -127 - lo.exponent;
227 Frac128 u = (shift < 128) ? Frac128((lo.mantissa >> shift).val) : Frac128(0);
228
229 // Evaluate 2^lo - 1 = lo * P(lo) using Horner scheme in 2 stages:
230 // - Degree 7-11: 64-bit precision evaluation.
231 // - Degree 0-6: 128-bit precision evaluation.
232 //
233 // Since u = |lo| <= 2^-7, truncation errors at degree i is bounded by:
234 // u^i <= 2^(-7i).
235 //
236 // Step 1: 64-bit evaluation, with truncation error at degree 7 is bounded by:
237 // 2^-63 * u^7 <= 2^-112.
238 Frac64 u64 = u.to_frac64();
239 Frac64 p64 = EXP2_POLY_128[11].to_frac64();
240 if (lo_is_neg) {
241 for (int i = 10; i >= 7; --i)
242 p64 = EXP2_POLY_128[i].to_frac64() - ((p64 * u64) << 1);
243 } else {
244 for (int i = 10; i >= 7; --i)
245 p64 = EXP2_POLY_128[i].to_frac64() + ((p64 * u64) << 1);
246 }
247
248 // Step 2: 128-bit evaluation.
249 Frac128 p({0, p64.val[0]});
250 if (lo_is_neg) {
251 for (int i = 6; i >= 0; --i)
252 p = EXP2_POLY_128[i] - ((p * u) << 1);
253 } else {
254 for (int i = 6; i >= 0; --i)
255 p = EXP2_POLY_128[i] + ((p * u) << 1);
256 }
257
258 // Reconstruction:
259 // t = u * P(u) ~ 2^u - 1 as Frac128
260 // mt = m * t ~ m * (2^u - 1) as Frac128
261 // m_final = m +- mt ~ m * (1 +- t) ~ m * 2^lo
262 Frac128 t = (u * p) << 1;
263 Frac128 mt = (m * t) << 1;
264 Frac128 m_final = lo_is_neg ? (m - mt) : (m + mt);
265
266 if ((m_final.val[1] & (1ULL << 63)) == 0) {
267 m_final = m_final << 1;
268 --hi;
269 }
270
271 return DFloat128(Sign::POS, hi - 127, m_final);
272}
273
274// Accurate pow(x, y) reusing range reduction parameters from the fast pass.
275LIBC_INLINE double pow_accurate(double x, double y, bool is_neg, int x_e,
276 unsigned idx_x, double dx) {
277 DFloat128 log2_x = log2_f128(x_e, idx_x, dx);
278 DFloat128 y_f128(y);
279 DFloat128 z = fputil::quick_mul(a: y_f128, b: log2_x);
280
281 // 2^1025 > max normal double.
282 double z_d = static_cast<double>(z);
283 if (LIBC_UNLIKELY(z_d >= 1025.0))
284 return set_overflow(is_neg);
285
286 // 2^-1076 < min subnormal double.
287 if (LIBC_UNLIKELY(z_d <= -1076.0))
288 return set_underflow(is_neg);
289
290 // For 0 < |z| <= 2^-55, x^y is between 1 - 2^-54 and 1 + 2^-53.
291 if (LIBC_UNLIKELY(!z.mantissa.is_zero() && z.exponent + 127 <= -55)) {
292 volatile double one = 1.0;
293 volatile double eps = (z.sign == Sign::NEG) ? -0x1.0p-100 : 0x1.0p-100;
294 double res = one + eps;
295 return is_neg ? -res : res;
296 }
297
298 DFloat128 r = exp2_f128(z);
299 if (is_neg)
300 r.sign = Sign::NEG;
301
302 int unbiased_exp = r.exponent + 127;
303 if (LIBC_UNLIKELY(unbiased_exp >= 1024))
304 return set_overflow(is_neg);
305
306 // Check if r is close to a 54-bit rounding boundary (either an exact 53-bit
307 // float or a midpoint).
308 // For normal numbers, the 54 leading bits occupy bits [127:74] of r.mantissa,
309 // leaving shift = 128 - 54 = 74 fractional bits below the boundary.
310 // For subnormal numbers (unbiased_exp < -1022), the boundary shifts right.
311 int shift = 74;
312 if (LIBC_UNLIKELY(unbiased_exp < -1022))
313 shift = 74 - (-1022 - unbiased_exp);
314
315 bool is_boundary_candidate = false;
316 // Lauter & Lefevre (2009) showed that if x^y is not an exact 54-bit number,
317 // the distance to the nearest 54-bit boundary is at least:
318 // |x^y - o_54(x^y)| / x^y >= 2^-114.
319 // In r.mantissa, this minimum distance is 2^-114 * 2^127 = 2^13.
320 //
321 // Since eps_bound >= 2^13, when shift < 14, the maximum possible distance to
322 // the nearest boundary, 2^(shift - 1) <= 2^12, is strictly less than
323 // eps_bound, making every value a candidate. Checking shift >= 14 also
324 // prevents undefined negative shifts for denormals.
325 if (shift >= 14) {
326 // Distance from r.mantissa to the nearest multiple of 2^shift:
327 MantissaType mask = (MantissaType(1) << shift) - 1;
328 MantissaType rem = r.mantissa & mask;
329 MantissaType half = MantissaType(1) << (shift - 1);
330 MantissaType dist =
331 (rem <= half) ? rem : ((MantissaType(1) << shift) - rem);
332
333 // Evaluation error in r.mantissa is bounded by:
334 // log(2) * |y| * AbsErr(log2_x) * 2^127 + AbsErr(exp2_z)
335 // <= log(2) * 32 * |y| + 2^13
336 // < 32 * |y| + 2^13.
337 // For |y| >= 2^59, 32 * |y| >= 2^64, so we cap eps_bound at 2^64 to avoid
338 // 64-bit integer overflow.
339 double abs_y = (y < 0.0) ? -y : y;
340 MantissaType eps_bound = MantissaType(1) << 13;
341 if (LIBC_LIKELY(abs_y < 0x1.0p59))
342 eps_bound += (MantissaType(static_cast<uint64_t>(abs_y)) << 5);
343 else
344 eps_bound = (MantissaType(1) << 64);
345 is_boundary_candidate = (dist < eps_bound);
346 } else {
347 is_boundary_candidate = true;
348 }
349
350 if (LIBC_UNLIKELY(is_boundary_candidate)) {
351 uint64_t exact_m = 0;
352 int exact_exp = 0;
353 if (is_exact_rounding_boundary(x, y, exact_m, exact_exp)) {
354 int l = 64 - cpp::countl_zero(value: exact_m);
355 DFloat128 exact_f128(r.sign, exact_exp + l - 128,
356 MantissaType(exact_m) << (128 - l));
357 exact_f128.normalize();
358
359 return static_cast<double>(exact_f128);
360 }
361
362 return pow_accurate_256(y, is_neg, x_e, idx_x, dx);
363 }
364
365 double res = static_cast<double>(r);
366 if (LIBC_UNLIKELY(fputil::FPBits<double>(res).is_inf()))
367 return set_overflow(is_neg);
368 return res;
369}
370
371} // namespace pow_internal
372} // namespace math
373} // namespace LIBC_NAMESPACE_DECL
374
375#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_ACCURATE_128_H
376