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/// Implementation header for double-precision pow(x, y).
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_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/config.h"
23#include "src/__support/macros/optimization.h"
24#include "src/__support/macros/properties/cpu_features.h"
25#include "src/__support/math/common_constants.h"
26#include "src/__support/math/exp_constants.h"
27#include "src/__support/math/pow_utils.h"
28
29#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
30#include "src/__support/math/pow_fast.h"
31#else
32
33#include "src/__support/math/pow_accurate_128.h"
34
35namespace LIBC_NAMESPACE_DECL {
36
37namespace math {
38
39// Overview of the main part of pow(x, y) = x^y computations.
40//
41// Let x = 2^(e_x) * m_x > 0. Then:
42// x^y = 2^( y * log2(x) )
43// = 2^( y * ( e_x + log2(m_x) ) )
44// = 2^( e_h + e_l )
45// = 2^(e_h) * 2^(e_l)
46// where:
47// e_h = round(y * log2(x)),
48// e_l = {y * log2(x)} = y * log2(x) - e_h.
49//
50// In particular, e_h is an integer, and |e_l| <= 0.5.
51//
52// For the final result to be fit in double precision, the exponent field can be
53// bounded by:
54// -1075 <= e_h <= 1024,
55// and anything outside that can be passed to quick overflow/underflow logic.
56//
57// Since |e_l| <= 0.5:
58// 0.5 < 2^(e_l) < 2,
59// the relative error of x^y = 2^(e_h) * 2^(e_l) is about
60// ~ absolute error of e_l
61// ~ absolute error of y * log2(x)
62// ~ relative error(log2(x)) * |y|
63// ~ relative error(log2(x)) * |e_h|
64// < relative error(log2(x)) * 2^11.
65//
66// Roughly speaking, to compute x^y with relative error < 2^(-n), we will need
67// to compute log2(x) with relative error < 2^(-n - 11), approximately.
68//
69// To compute log2(x), we will perform range reduction for log2(m_x):
70// dx = r * m_x - 1.
71// Then m_x = (1 + dx) / r, and
72// log2(m_x) = log2( (1 + dx) / r )
73// = log2(1 + dx) - log2(r),
74// where -log2(r) is obtained from look up tables.
75// The computation of dx = r * m_x - 1 is exact, and we choose the size of the
76// look up table for r's such that:
77// -2^-8 <= dx < 2^-7.
78// Combining them together, we have that:
79// log2(x) = e_x - log2(r) + log2(1 + dx)
80// = e_x - log2(r) + log2(e) * (dx - dx^2/2 + dx^3/3 - ...)
81// So in the worst case, where e_x = log2(r) = 0, the relative error of log2(x)
82// computation will be:
83// ~ absolute_error(log2(1 + dx)) / |dx|.
84// So if we compute log2(1 + dx) accurately up to dx^n term, and use a
85// polynomial approximation for the dx^(n + 1) and higher terms:
86// log2(1 + dx) ~ log2(e) * (dx - dx^2/2 + ... + (-1)^(n + 1) dx^n / n) +
87// dx^(n + 1) * P(dx)
88// such that the absolute error of P(dx) is smaller than double precision ulp,
89// then the overall relative error of our log2(1 + dx) approximation is:
90// ~ ulp(dx^(n + 1)) / |dx| ~ 2^(-52) * |dx^n|
91//
92// Let's consider 3 cases:
93// - For a fast path of a correctly rounded implementation, we will want the
94// relative error of x^y bounded above by:
95// ~ 2^(-precision - 10) = 2^(-53 - 10) = 2^(-63)
96// so that the Ziv accuracy test passes with probability > 1 - 2^-10 ~ 99.99%.
97// And from the above argument, we need to compute log2(1 + dx) with relative
98// error ~ 2^(-63 - 11) = 2^(-74).
99// So we will need to compute log2(1 + dx) accurately up to dx^n such that:
100// 2^(-52) * |dx^n| ~ 2^(-74).
101// Or equivalently:
102// |dx^n| < 2^(-22).
103// With our bounds from range reduction |dx| < 2^(-7), n = 3 will be quite
104// close to our desired precision needed.
105//
106// In summary, for a fast path of correctly rounded implementation, we will
107// compute log2(e) * (dx - dx^2/2 + dx^3/3) accurately, and approximate higher
108// terms with dx^4 * P(dx), where:
109// |(log2(1 + dx) - log2(e) * (dx - dx^2/2 + dx^3/3))/dx^4 - P(dx)| < 2^-53.
110//
111// - For a strictly < 1 ULP error fast version across the full exponent range
112// (as implemented in pow_fast), we will want the relative error of x^y
113// bounded above by:
114// ~ 2^(-precision) = 2^(-53).
115// So the relative error needed to approximate log2(1 + dx) is:
116// ~ 2^(-53 - 11) = 2^(-64),
117// and we will need to compute log2(1 + dx) accurately up to dx^n such that:
118// 2^(-52) * |dx^n| ~ 2^(-64),
119// or equivalently:
120// |dx^n| < 2^-12.
121// With |dx| < 2^(-7), n = 2 satisfies |dx^2| < 2^(-14) < 2^(-12).
122// Hence, we compute log2(e) * (dx - dx^2/2) accurately in DoubleDouble, and
123// approximate higher terms with dx^3 * P(dx), where:
124// |(log2(1 + dx) - log2(e) * (dx - dx^2/2))/dx^3 - P(dx)| < 2^-53.
125// This guarantees error < 1 ULP everywhere (empirically bounded by 0.75 ULP).
126//
127// - If we were to choose n = 1 (computing only log2(e) * dx accurately, and
128// approximating higher terms with dx^2 * P(dx)), the relative error of
129// log2(1 + dx) is only bounded by:
130// ~ 2^(-52) * |dx| < 2^(-52) * 2^(-7) = 2^(-59).
131// The resulting relative error of x^y is about:
132// ~ 2^(-59) * |e_h|.
133// While this achieves < 1 ULP for typical inputs with |e_h| <= 64, near
134// extreme exponent boundaries (|e_h| ~ 2^10 and |dx| ~ 2^-7) the error can
135// drift up to ~7.8 ULPs (bounded theoretically by ~16 ULPs).
136
137LIBC_INLINE double pow(double x, double y) {
138 using namespace pow_internal;
139 using FPBits = fputil::FPBits<double>;
140
141 FPBits xbits(x), ybits(y);
142 uint64_t x_u = xbits.uintval();
143 uint64_t y_u = ybits.uintval();
144 uint64_t y_a = ybits.abs().uintval();
145
146 if (LIBC_UNLIKELY((x_u & 0x0003'FFFF'FFFF'FFFF) == 0 ||
147 (y_u & 0x000F'FFFF'FFFF'FFFF) == 0)) {
148 if (auto r = check_special_inputs(x, y); LIBC_UNLIKELY(r.has_value()))
149 return r.value();
150 }
151
152 double e_x = static_cast<double>(xbits.get_exponent());
153 uint64_t x_mant = xbits.get_mantissa();
154 bool is_neg = false;
155 double sign_d = 1.0;
156
157 if (LIBC_UNLIKELY(y_a <= Y_LOWER_BOUND || y_a >= Y_UPPER_BOUND ||
158 x_u >= FPBits::inf().uintval() ||
159 x_u < FPBits::min_normal().uintval())) {
160 if (auto r = check_exceptional_cases(x, y, e_x, x_mant, is_neg, sign_d);
161 LIBC_UNLIKELY(r.has_value()))
162 return r.value();
163 }
164
165 // x^y = 2^( y * log2(x) )
166 // = 2^( y * ( e_x + log2(m_x) ) )
167 // First we compute log2(x) = e_x + log2(m_x)
168
169 // Extract exponent field of x.
170
171 // Use the highest 7 fractional bits of m_x as the index for look up tables.
172 unsigned idx_x = static_cast<unsigned>(x_mant >> (FPBits::FRACTION_LEN - 7));
173 // Add the hidden bit to the mantissa.
174 // 1 <= m_x < 2
175 FPBits m_x = FPBits(x_mant | 0x3ff0'0000'0000'0000);
176
177 // Reduced argument for log2(m_x):
178 // dx = r * m_x - 1.
179 // The computation is exact, and -2^-8 <= dx < 2^-7.
180 // Then m_x = (1 + dx) / r, and
181 // log2(m_x) = log2( (1 + dx) / r )
182 // = log2(1 + dx) - log2(r).
183
184 // As analyzed in the overview comment above, we evaluate the cubic part
185 // (dx * C1 + dx^2 * C2 + dx^3 * C3) accurately using DoubleDouble (n = 3)
186 // and approximate higher terms with dx^4 * P(dx) to ensure relative error
187 // < 2^-63 for this fast path.
188
189 // Degree-5 polynomial approximation for:
190 // P(dx) ~ (log2(1 + dx) - (dx - dx^2/2 + dx^3/3)/log(2)) / dx^4
191 // Generated by Sollya with:
192 // > P = fpminimax((log2(1 + x) - (x - x^2/2 + x^3/3)/log(2))/x^4, 5,
193 // [|D...|], [-2^-8, 2^-7]);
194 // > dirtyinfnorm((log2(1 + x) - (x - x^2/2 + x^3/3)/log(2))/x - x^3*P,
195 // [-2^-8, 2^-7]);
196 // 0x1.a643c...p-74
197 // > dirtyinfnorm((log2(1 + x) - (x - x^2/2 + x^3/3)/log(2))/x^4 - P,
198 // [-2^-8, 2^-7]);
199 // 0x1.b81c5...p-53
200 constexpr double COEFFS[] = {-0x1.71547652b82fdp-2, 0x1.2776c50ef8f9bp-2,
201 -0x1.ec709dc4f0fedp-3, 0x1.a617677f716dep-3,
202 -0x1.715423d54d5p-3, 0x1.44e6355fc4d03p-3};
203
204 // Constants for C_k = (-1)^(k-1) / (k * log(2)) in DoubleDouble:
205 // C1 = 1 / log(2)
206 constexpr DoubleDouble C1 = {.lo: 0x1.777d0ffda0d24p-56, .hi: 0x1.71547652b82fep0};
207 // C2 = -1 / (2 * log(2))
208 constexpr DoubleDouble C2 = {.lo: -0x1.777d0ffda0d24p-57, .hi: -0x1.71547652b82fep-1};
209 // C3 = 1 / (3 * log(2))
210 constexpr DoubleDouble C3 = {.lo: 0x1.b749fc15522bcp-50, .hi: 0x1.ec709dc3a03e2p-2};
211
212 // Perform exact range reduction.
213#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
214 double dx = fputil::multiply_add(RD[idx_x], m_x.get_val(), -1.0); // Exact
215#else
216 double c = FPBits(m_x.uintval() & 0x3fff'e000'0000'0000).get_val();
217 double dx =
218 fputil::multiply_add(x: RD[idx_x], y: m_x.get_val() - c, z: CD[idx_x]); // Exact
219#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
220
221 // Evaluate the cubic part (dx * C1 + dx^2 * C2 + dx^3 * C3) and polynomial
222 // tail using a parallel Double-Double Estrin scheme.
223#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
224 // Error-free transformation for r = C1.hi + dx * C2.hi:
225 double r_hi = fputil::multiply_add(dx, C2.hi, C1.hi);
226 double r_lo = fputil::multiply_add(dx, C2.hi, C1.hi - r_hi); // Exact error
227#else
228 DoubleDouble dx_c2 = fputil::exact_mult(a: dx, b: C2.hi);
229 DoubleDouble r_sum = fputil::exact_add(a: C1.hi, b: dx_c2.hi);
230 double r_hi = r_sum.hi;
231 double r_lo = r_sum.lo + dx_c2.lo;
232#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
233
234 // Low parts polynomial evaluated in parallel:
235 // C1.lo + dx * (C2.lo + dx * C3.lo) + r_lo
236 double lo_tail = fputil::multiply_add(x: dx, y: C3.lo, z: C2.lo);
237 double lo_poly = fputil::multiply_add(x: dx, y: lo_tail, z: C1.lo + r_lo);
238
239 // Evaluate polynomial tail P(dx) using Estrin's scheme:
240 double dx2 = dx * dx;
241 double c0 = fputil::multiply_add(x: dx, y: COEFFS[1], z: COEFFS[0]);
242 double c1 = fputil::multiply_add(x: dx, y: COEFFS[3], z: COEFFS[2]);
243 double c2 = fputil::multiply_add(x: dx, y: COEFFS[5], z: COEFFS[4]);
244
245 double dx4 = dx2 * dx2;
246 double d0 = fputil::multiply_add(x: dx2, y: c1, z: c0);
247 double p = fputil::multiply_add(x: dx4, y: c2, z: d0);
248
249 // High part of cubic term: C3.hi + dx * P(dx)
250 double q = fputil::multiply_add(x: dx, y: p, z: C3.hi);
251
252 // Combine dx^2 * q into lo_poly:
253 lo_poly = fputil::multiply_add(x: dx2, y: q, z: lo_poly);
254
255 // Multiply by dx to get log2(1 + dx) in DoubleDouble:
256 DoubleDouble log2_1p = fputil::exact_mult(a: dx, b: r_hi);
257 log2_1p.lo = fputil::multiply_add(x: dx, y: lo_poly, z: log2_1p.lo);
258
259 // Combine with e_x - log2(r):
260 DoubleDouble log2_x_hi =
261 fputil::exact_add(a: e_x + LOG2_R_DD[idx_x].hi, b: log2_1p.hi);
262 double log2_x_lo = log2_1p.lo + LOG2_R_DD[idx_x].lo;
263 DoubleDouble log2_x = fputil::exact_add(a: log2_x_hi.hi, b: log2_x_lo);
264 log2_x.lo += log2_x_hi.lo;
265
266 // To compute 2^(y * log2(x)), we break the exponent into 3 parts:
267 // y * log2(x) = hi + mid + lo, where
268 // hi is an integer
269 // mid * 2^6 is an integer
270 // |lo| <= 2^-7
271 // Then:
272 // x^y = 2^(y * log2(x)) = 2^hi * 2^mid * 2^lo,
273 // In which 2^mid is obtained from a look-up table of size 2^6 = 64 elements,
274 // and 2^lo ~ 1 + lo * P(lo).
275 // Thus, we have:
276 // hi + mid = 2^-6 * round( 2^6 * y * log2(x) )
277 // If we restrict the output such that |hi| < 512, (hi + mid) uses (9 + 6)
278 // bits, hence, if we use double precision to perform
279 // round( 2^6 * y * log2(x))
280 // the lo part is bounded by 2^-7 + 2^(-(52 - 15)) = 2^-7 + 2^-37
281
282 // In the following computations:
283 // y6 = 2^6 * y
284 // hm = 2^6 * (hi + mid) = round(2^6 * y * log2(x)) ~ round(y6 * s)
285 // lo6 = 2^6 * lo = 2^6 * (y - (hi + mid)) = y6 * log2(x) - hm.
286 constexpr double SCALE = 0x1.0p6;
287 double y6 = y * SCALE; // Exact.
288
289 DoubleDouble y6_log2_x = fputil::exact_mult(a: y6, b: log2_x.hi);
290 y6_log2_x.lo = fputil::multiply_add(x: y6, y: log2_x.lo, z: y6_log2_x.lo);
291
292 // Check overflow/underflow.
293 double scale = 1.0;
294 bool is_denorm = false;
295
296 // |2^(hi + mid) - exp2_hi_mid| <= ulp(exp2_hi_mid) / 2
297
298 // The fast computation for 2^hi below requires that:
299 // |hi| < 512, or equivalently, |hm| < 512 * 2^6.
300 // This guarantees that the biased exponent:
301 // exp_biased = (hm_i >> 6) + EXP_BIAS = hi + 1023
302 // is strictly within the normal double range:
303 // 1023 - 511 <= exp_biased <= 1023 + 511, or 512 <= exp_biased <= 1534.
304 // Hence, 2^hi is always a normal, non-zero, finite power of 2, and the
305 // multiplication upper * exp2_hi is exact and never underflows or overflows.
306 //
307 // From the edge case checks above:
308 // 2^(-54) / 1074 <= |y| <= 1075 * 2^53, and
309 // |log_2(1 - 2^(-53))| <= |log2(x)| <= 1074.
310 // So their product is bounded by:
311 // 2^-117 < |y * log2(x)| < 2^74.
312 //
313 // The meaningful range for y * log2(x) in double precision is:
314 // -1075 <= y * log2(x) <= 1024.
315 // Any value > 1024 overflows, and any value < -1075 underflows.
316 //
317 // When |y * log2(x)| >= 511, we shift the exponent by an offset S:
318 // y * log2(x) = (y * log2(x) - S) + S
319 // and multiply by scale = 2^S at the end.
320 // To ensure the shifted exponent (y * log2(x) - S) stays within [-511, 511]:
321 // - For positive range [511, 1024]:
322 // 511 - S > -511 ==> S < 1022
323 // 1024 - S < 511 ==> S > 513
324 // - For negative range [-1076, -511]:
325 // -511 + S < 511 ==> S <= 1022
326 // -1076 + S > -511 ==> S >= 565
327 // Combined, the shift must satisfy: 565 <= S <= 1022.
328 // We choose S = 600, which yields:
329 // y * log2(x) - 600 in [-89, 424] for the positive range, and
330 // y * log2(x) + 600 in [-476, 89] for the negative range,
331 // both fitting well within [-511, 511].
332 //
333 // For exponents that are completely out of range:
334 // y * log2(x) > 1025 or y * log2(x) < -1076,
335 // the product can be as large as 2^74, so subtracting or adding 600 * 64
336 // would still overflow a 32-bit int when computing hm_i.
337 // We return early with overflow or underflow in these cases.
338 //
339 // Alternatively, for less branching or in SIMD / vector implementations, one
340 // could clamp y6_log2_x.hi to:
341 // - UPPER_EXP_BOUND (511 * 64) for overflow, which when multiplied by
342 // scale = 2^600 yields 2^1111 and correctly overflows.
343 // - -500 * 64 for underflow, which when multiplied by scale = 2^-600 yields
344 // 2^-1100 and correctly underflows.
345 // However, in this correctly rounded version, such clamping can cause the
346 // Ziv accuracy test to fail on the clamped exponent, unnecessarily
347 // triggering the slow accurate paths for overflow or underflow cases.
348
349 constexpr double UPPER_EXP_BOUND = 511.0 * SCALE;
350 if (LIBC_UNLIKELY(FPBits(y6_log2_x.hi).abs().get_val() >= UPPER_EXP_BOUND)) {
351 if (FPBits(y6_log2_x.hi).sign() == Sign::POS) {
352 if (y6_log2_x.hi > 1025.0 * SCALE)
353 return set_overflow(is_neg);
354 scale = 0x1.0p600;
355 y6_log2_x.hi -= 600.0 * SCALE;
356 } else {
357 if (y6_log2_x.hi <= -1021.0 * SCALE) {
358 if (y6_log2_x.hi < -1076.0 * SCALE)
359 return set_underflow(is_neg);
360 is_denorm = true;
361 } else {
362 scale = 0x1.0p-600;
363 y6_log2_x.hi += 600.0 * SCALE;
364 }
365 }
366 }
367
368 double hm = fputil::nearest_integer(x: y6_log2_x.hi);
369
370 // lo6 = 2^6 * lo.
371 DoubleDouble lo6 = fputil::exact_add(a: y6_log2_x.hi - hm, b: y6_log2_x.lo);
372
373 int hm_i = static_cast<int>(hm);
374 unsigned idx_y = static_cast<unsigned>(hm_i) & 0x3f;
375
376 // 2^hi
377 int hi = hm_i >> 6;
378 double exp2_hi = 1.0;
379 if (LIBC_LIKELY(!is_denorm)) {
380 int64_t exp2_hi_i = static_cast<int64_t>(
381 static_cast<uint64_t>(hi + FPBits::EXP_BIAS) << FPBits::FRACTION_LEN);
382 exp2_hi = FPBits(static_cast<uint64_t>(exp2_hi_i)).get_val();
383 }
384
385 // 2^mid
386 DoubleDouble exp2_mid{.lo: EXP2_MID1[idx_y].mid * sign_d,
387 .hi: EXP2_MID1[idx_y].hi * sign_d};
388
389 // Polynomial expansion for 2^(lo6/64):
390 // 2^(lo6/64) ~ 1 + lo6 * (log(2)/64) + lo6^2 * P(lo6)
391 // The linear term is computed in DoubleDouble, and the degree-3 polynomial
392 // P(lo6) is evaluated with standard double precision.
393 //
394 // hi and lo parts of log(2)/64, generated by Sollya with:
395 // > a = D(log(2)/64);
396 // > b = D(log(2)/64 - a);
397 constexpr DoubleDouble LOG_2_OVER_64 = {.lo: 0x1.abc9e3b39803fp-62,
398 .hi: 0x1.62e42fefa39efp-7};
399
400 // Degree-3 polynomial approximation for (2^(lo6/64) - 1 - lo6*log(2)/64) /
401 // lo6^2: Generated by Sollya with:
402 // > f = (2^(x/64) - 1 - x*log(2)/64) / x^2;
403 // > P = fpminimax(f, 5, [|D...|], [-0.5, 0.5]);
404 // > dirtyinfnorm((f - P) * x^2, [-0.5, 0.5]);
405 // 0x1.1778...p-73
406 constexpr double EXP2_COEFFS[] = {
407 0x1.ebfbdff82c58fp-15, 0x1.c6b08d704a0cp-23, 0x1.3b2ab6fb4d08fp-31,
408 0x1.5d87fe77a735dp-40, 0x1.430d835610044p-49, 0x1.ffe67c38112c3p-59};
409
410 DoubleDouble lo_log_2 = fputil::quick_mult(a: lo6, b: LOG_2_OVER_64);
411 DoubleDouble lo_log_2_p1 = fputil::exact_add(a: 1.0, b: lo_log_2.hi);
412 lo_log_2_p1.lo += lo_log_2.lo;
413
414 double lo6_sq = lo6.hi * lo6.hi;
415 double e0 = fputil::multiply_add(x: lo6.hi, y: EXP2_COEFFS[1], z: EXP2_COEFFS[0]);
416 double e1 = fputil::multiply_add(x: lo6.hi, y: EXP2_COEFFS[3], z: EXP2_COEFFS[2]);
417 double e2 = fputil::multiply_add(x: lo6.hi, y: EXP2_COEFFS[5], z: EXP2_COEFFS[4]);
418
419 double lo6_4 = lo6_sq * lo6_sq;
420 double f0 = fputil::multiply_add(x: lo6_sq, y: e0, z: lo_log_2_p1.lo);
421 double f1 = fputil::multiply_add(x: lo6_sq, y: e2, z: e1);
422
423 lo_log_2_p1.lo = fputil::multiply_add(x: lo6_4, y: f1, z: f0);
424 DoubleDouble r = fputil::quick_mult(a: exp2_mid, b: lo_log_2_p1);
425
426 // Absolute error bound for r:
427 // - Base error from exp2 stage is bounded by 0x1.0p-64.
428 // - Propagated error from log2(x):
429 // |y| * err(log2(x)) * (2^mid * log(2)) <= |y| * 2^-73.5
430 // Dynamic error bound adapting to exponent scale |y|:
431#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
432 double err_r =
433 fputil::multiply_add(FPBits(y).abs().get_val(), 0x1.8p-73, 0x1.0p-64);
434#else // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
435 // Without FMA, intermediate roundings increase log2(x) and exp2 errors.
436 double err_r =
437 fputil::multiply_add(x: FPBits(y).abs().get_val(), y: 0x1.cp-72, z: 0x1.0p-63);
438#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
439
440 if (LIBC_UNLIKELY(is_denorm)) {
441 if (auto r_denorm = ziv_test_denorm(hi, mid: r.hi, lo: r.lo, err: err_r, is_neg);
442 LIBC_LIKELY(r_denorm.has_value())) {
443 double res = r_denorm.value();
444
445 // Since is_denorm is triggered when y * log2(x) <= -1021 * 2^6, the
446 // rounded result can still be a normal number (>= 2^-1022). In that
447 // case, no underflow has occurred, so return directly.
448 if (LIBC_UNLIKELY(FPBits(res).is_normal()))
449 return res;
450
451 // When the rounded result is denormal or zero, underflow and ERANGE
452 // should only be set if the result is inexact. We check for exact
453 // results:
454 // 1. If x = 2^e_x (x_mant == 0), then x^y = 2^(e_x * y) is an exact
455 // power of 2 iff e_x * y is an integer and >= -1074.0 (the smallest
456 // representable power of 2 in double precision).
457 if (LIBC_UNLIKELY(x_mant == 0)) {
458 double ex_y = e_x * y;
459 if (ex_y >= -1074.0 && pow_internal::is_integer(x: ex_y))
460 return res;
461 } else if (LIBC_UNLIKELY(y > 0.0 && y <= 35.0)) {
462 // 2. If x is not a power of 2, exact results in double precision can
463 // only occur for 0 < y <= 35 (Lauter and Lefevre). If it falls on
464 // an exact boundary, convert via DyadicFloat to ensure that
465 // underflow and inexact exceptions are not signaled.
466 uint64_t exact_m = 0;
467 int exact_exp = 0;
468 if (pow_internal::is_exact_rounding_boundary(x, y, exact_m,
469 exact_exp)) {
470 int l = 64 - cpp::countl_zero(value: exact_m);
471 pow_internal::DFloat128 exact_f128(
472 is_neg ? Sign::NEG : Sign::POS, exact_exp + l - 128,
473 pow_internal::MantissaType(exact_m) << (128 - l));
474 exact_f128.normalize();
475 return static_cast<double>(exact_f128);
476 }
477 }
478
479 // Otherwise, the denormal or zero result is inexact, so we signal the
480 // underflow exception and set ERANGE if required.
481 fputil::set_errno_if_required(ERANGE);
482 fputil::raise_underflow_except_if_required<double>();
483 return res;
484 }
485 return pow_accurate(x, y, is_neg, x_e: static_cast<int>(e_x), idx_x, dx);
486 }
487
488 double upper = r.hi + (r.lo + err_r);
489 double lower = r.hi + (r.lo - err_r);
490 if (LIBC_LIKELY(upper == lower)) {
491 double tmp1 = upper * exp2_hi;
492 return tmp1 * scale;
493 }
494
495 return pow_accurate(x, y, is_neg, x_e: static_cast<int>(e_x), idx_x, dx);
496}
497
498} // namespace math
499} // namespace LIBC_NAMESPACE_DECL
500
501#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
502
503#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_H
504