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 evaluation implementation for powf(x, y).
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
16
17#include "src/__support/CPP/bit.h"
18#include "src/__support/CPP/optional.h"
19#include "src/__support/FPUtil/FEnvImpl.h"
20#include "src/__support/FPUtil/FPBits.h"
21#include "src/__support/FPUtil/double_double.h"
22#include "src/__support/FPUtil/multiply_add.h"
23#include "src/__support/FPUtil/nearest_integer.h"
24#include "src/__support/common.h"
25#include "src/__support/macros/config.h"
26#include "src/__support/macros/optimization.h"
27#include "src/__support/macros/properties/cpu_features.h"
28#include "src/__support/math/common_constants.h"
29#include "src/__support/math/exp_constants.h"
30#include "src/__support/math/powf_utils.h"
31
32namespace LIBC_NAMESPACE_DECL {
33namespace math {
34namespace double_eval {
35
36namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;
37
38#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
39// Accurate evaluation in DoubleDouble precision when Ziv's test fails.
40LIBC_INLINE double powf_accurate(float x, float y, double e_x, uint64_t sign) {
41 using FloatBits = fputil::FPBits<float>;
42 using DoubleBits = fputil::FPBits<double>;
43 using fputil::DoubleDouble;
44
45 uint32_t x_u = FloatBits(x).uintval();
46 uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
47 int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 7));
48 // Append hidden bit.
49 float m_x = cpp::bit_cast<float>(from: x_mant | 0x3f80'0000);
50 double dx = fputil::multiply_add(
51 x: static_cast<double>(m_x),
52 y: static_cast<double>(common_constants_internal::R[idx_x]), z: -1.0);
53
54 // Degree-11 minimax polynomial for log2(1 + x)/x on [-2^-8, 2^-7] in
55 // double-double, generated by Sollya with:
56 // > display = hexadecimal;
57 // > prec = 256;
58 // > P = fpminimax(log2(1 + x)/x, 11, [|DD...|], [-2^(-8), 2^(-7)]);
59 // > for i from 0 to 11 do {
60 // c = coeff(P, i);
61 // hi = round(c, D, RN);
62 // lo = round(c - hi, D, RN);
63 // print("{", lo, ", ", hi, "},");
64 // };
65 // > dirtyinfnorm((log2(1 + x) - x * P) / x, [-2^(-8), 2^(-7)]);
66 // 0x1.c598...p-104 < 2^-103.
67 constexpr DoubleDouble LOG2_COEFFS[] = {
68 {.lo: 0x1.777d0ffda0d34p-56, .hi: 0x1.71547652b82fep0},
69 {.lo: -0x1.777d0ffd88ca8p-57, .hi: -0x1.71547652b82fep-1},
70 {.lo: 0x1.d27f052d9d41dp-56, .hi: 0x1.ec709dc3a03fdp-2},
71 {.lo: -0x1.777f259c6c0ap-58, .hi: -0x1.71547652b82fep-2},
72 {.lo: 0x1.e5c055599afdbp-56, .hi: 0x1.2776c50ef9bfep-2},
73 {.lo: 0x1.e923056eef156p-58, .hi: -0x1.ec709dc3a03fdp-3},
74 {.lo: 0x1.1f24e26b5ad0dp-57, .hi: 0x1.a61762a7add76p-3},
75 {.lo: 0x1.16d9b5013b55ep-58, .hi: -0x1.71547652bf304p-3},
76 {.lo: -0x1.1551148c3d1d8p-57, .hi: 0x1.484b13f00425ep-3},
77 {.lo: -0x1.a1cf327c6703p-60, .hi: -0x1.2776ca937f335p-3},
78 {.lo: -0x1.d24f6254c2e6ap-57, .hi: 0x1.0c9211efe6d7ep-3},
79 {.lo: 0x1.a81cff2946d7cp-59, .hi: -0x1.e1f45baaf7f1dp-4},
80 };
81
82 // Evaluate P(dx) ~ log2(1 + dx) / dx using Horner's scheme.
83 DoubleDouble p = LOG2_COEFFS[11];
84 for (int i = 10; i >= 0; --i)
85 p = fputil::add(a: LOG2_COEFFS[i], b: fputil::quick_mult(a: dx, b: p));
86
87 // log2(1 + dx) = dx * P(dx) in DoubleDouble.
88 DoubleDouble log2_1p = fputil::quick_mult(a: dx, b: p);
89
90 // -log2(r) is represented in TripleDouble as {lo, mid, hi}.
91 // We decompose log2(x) = e_x - log2(r) + log2(1 + dx) into:
92 // log2(x) = (e_x + hi) + (log2_1p + mid + lo)
93 // = (e_x + hi) + log2_tail.
94 // The lower part log2_tail accumulates log2_1p, mid, and lo in DoubleDouble:
95 DoubleDouble log2_tail =
96 fputil::add<false>(a: log2_1p, b: powf_internal::LOG2_R_TD[idx_x].mid);
97 log2_tail = fputil::add<false>(a: log2_tail, b: powf_internal::LOG2_R_TD[idx_x].lo);
98
99 // Range reduction for 2^(y * log2(x)):
100 // Scale by 2^6 = 64 to use the 64-entry EXP2_MID1 table:
101 // 64 * y * log2(x) = y6 * (e_x + hi) + y6 * log2_tail.
102 //
103 // Since e_x + hi has at most 8 + 20 = 28 bits of precision and y has 24 bits,
104 // the product y6 * (e_x + hi) has at most 28 + 24 = 52 bits, making y6_hi
105 // exact in double precision.
106 constexpr double SCALE = 0x1.0p6;
107 double y6 = static_cast<double>(y) * SCALE;
108 DoubleDouble prod = fputil::quick_mult(a: y6, b: log2_tail);
109
110 double e_x_hi = e_x + powf_internal::LOG2_R_TD[idx_x].hi;
111 DoubleDouble y6_hi = fputil::exact_mult(a: y6, b: e_x_hi);
112
113 // hm = round(64 * y * log2(x)) is the nearest integer index.
114 double hm_d = y6_hi.hi + prod.hi;
115 double hm = fputil::nearest_integer(x: hm_d);
116
117 // lo6 = 64 * y * log2(x) - hm = (y6_hi - hm) + prod:
118 DoubleDouble lo6_hi = fputil::exact_add<false>(a: y6_hi.hi, b: -hm);
119 lo6_hi.lo += y6_hi.lo;
120
121 DoubleDouble lo6 = fputil::add<false>(a: prod, b: lo6_hi);
122
123 // lo = lo6 / 64, with |lo| <= 2^-7.
124 DoubleDouble lo = fputil::quick_mult(a: 0x1.0p-6, b: lo6);
125
126 int64_t hm_i = static_cast<int64_t>(hm);
127 int idx_y = static_cast<int>(hm_i & 0x3f);
128 int hi = static_cast<int>(hm_i >> 6);
129 int64_t exp_hi_i = static_cast<int64_t>(static_cast<uint64_t>(hi)
130 << DoubleBits::FRACTION_LEN);
131
132 // Degree-8 minimax polynomial for (2^x - 1)/x on [-2^-7, 2^-7] in
133 // double-double, generated by Sollya with:
134 // > display = hexadecimal;
135 // > prec = 256;
136 // > P = fpminimax((2^x - 1)/x, 8, [|DD...|], [-2^(-7), 2^(-7)]);
137 // > for i from 0 to 8 do {
138 // c = coeff(P, i);
139 // hi = round(c, D, RN);
140 // lo = round(c - hi, D, RN);
141 // print("{", lo, ", ", hi, "},");
142 // };
143 // > dirtyinfnorm((2^x - (1 + x * P)) / (2^x), [-2^(-7), 2^(-7)]);
144 // 0x1.e62d...p-106 < 2^-105.
145 constexpr DoubleDouble EXP2_COEFFS_DD[] = {
146 {.lo: 0x1.abc9e3b39803ep-56, .hi: 0x1.62e42fefa39efp-1},
147 {.lo: -0x1.5e43a54066432p-57, .hi: 0x1.ebfbdff82c58fp-3},
148 {.lo: -0x1.d331626e9fe04p-59, .hi: 0x1.c6b08d704a0cp-5},
149 {.lo: 0x1.4f492022eeb74p-62, .hi: 0x1.3b2ab6fba4e77p-7},
150 {.lo: 0x1.035733b00abdp-66, .hi: 0x1.5d87fe78a6731p-10},
151 {.lo: 0x1.9e4f81ec951dep-68, .hi: 0x1.430912f86c123p-13},
152 {.lo: 0x1.7be627c299948p-74, .hi: 0x1.ffcbfc588b997p-17},
153 {.lo: -0x1.c8a006d5b39acp-74, .hi: 0x1.62c03345a64f6p-20},
154 {.lo: -0x1.fd14495e4b96p-81, .hi: 0x1.b5253eef8c0a9p-24},
155 };
156
157 DoubleDouble exp2_p = EXP2_COEFFS_DD[8];
158 for (int i = 7; i >= 0; --i)
159 exp2_p = fputil::add(a: EXP2_COEFFS_DD[i], b: fputil::quick_mult(a: lo, b: exp2_p));
160
161 DoubleDouble p_lo = fputil::quick_mult(a: lo, b: exp2_p);
162 DoubleDouble exp2_lo = fputil::exact_add(a: 1.0, b: p_lo.hi);
163 exp2_lo.lo += p_lo.lo;
164
165 DoubleDouble exp2_hi_mid_dd;
166 exp2_hi_mid_dd.hi = cpp::bit_cast<double>(
167 from: exp_hi_i + cpp::bit_cast<int64_t>(from: EXP2_MID1[idx_y].hi) + sign);
168 exp2_hi_mid_dd.lo =
169 (idx_y != 0)
170 ? cpp::bit_cast<double>(
171 from: exp_hi_i + cpp::bit_cast<int64_t>(from: EXP2_MID1[idx_y].mid) + sign)
172 : 0.0;
173
174 DoubleDouble rr = fputil::quick_mult(a: exp2_hi_mid_dd, b: exp2_lo);
175
176 DoubleDouble r = fputil::exact_add(a: rr.hi, b: rr.lo);
177
178 // Round to odd to avoid double rounding when converting to float.
179 uint64_t r_bits = cpp::bit_cast<uint64_t>(from: r.hi);
180 if (LIBC_UNLIKELY(((r_bits & 0x0fff'ffffULL) == 0) && (r.lo != 0.0))) {
181 if (DoubleBits(r.hi).sign() == DoubleBits(r.lo).sign()) {
182 ++r_bits;
183 } else if ((r_bits & DoubleBits::FRACTION_MASK) > 0) {
184 --r_bits;
185 }
186 }
187
188 return cpp::bit_cast<double>(from: r_bits);
189}
190#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
191
192#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \
193 !defined(LIBC_MATH_HAS_SMALL_TABLES)
194// Check if x^y is an exact rounding boundary.
195// When x^y = exact_m * 2^exact_exp with exact_m <= 2^25:
196// - exact_m has at most 26 bits, so it fits in double precision without
197// rounding.
198// - For single precision, exact_exp is in [-200, 127], well within the normal
199// range [-1022, 1023] of double precision.
200// We compute exact_d = exact_m * 2^exact_exp in double precision and cast to
201// float for the final rounded result.
202LIBC_INLINE cpp::optional<float>
203check_exact_boundary(float x, float y, uint64_t sign,
204 Sign out_sign = Sign::POS) {
205 uint32_t exact_m = 0;
206 int exact_exp = 0;
207 if (LIBC_UNLIKELY(powf_internal::is_exact_rounding_boundary(x, y, exact_m,
208 exact_exp))) {
209 // Number of bits in exact_m: 2^(l - 1) <= exact_m < 2^l.
210 int l = 32 - cpp::countl_zero(value: exact_m);
211 // Unbiased exponent of exact_m * 2^exact_exp.
212 int unbiased_exp = exact_exp + l - 1;
213
214 if (LIBC_UNLIKELY(unbiased_exp > 127))
215 return powf_internal::set_overflow(out_sign);
216
217 // Below the minimum subnormal float (2^-149).
218 if (LIBC_UNLIKELY(unbiased_exp < -149)) {
219 fputil::set_errno_if_required(ERANGE);
220 fputil::raise_underflow_except_if_required<float>();
221 }
222
223 // Scale factor = sign * 2^exact_exp.
224 using DoubleBits = fputil::FPBits<double>;
225 double scale = cpp::bit_cast<double>(
226 from: (static_cast<uint64_t>(exact_exp + DoubleBits::EXP_BIAS)
227 << DoubleBits::FRACTION_LEN) |
228 sign);
229 double exact_d = static_cast<double>(exact_m) * scale;
230 return static_cast<float>(exact_d);
231 }
232 return cpp::nullopt;
233}
234#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS && !LIBC_MATH_HAS_SMALL_TABLES
235
236// Overview of powf(x, y) = x^y computation in double precision:
237//
238// Let x = 2^(e_x) * m_x > 0. Then:
239// x^y = 2^( y * log2(x) )
240// = 2^( y * ( e_x + log2(m_x) ) )
241// = 2^( k + f )
242// = 2^k * 2^f,
243// where:
244// k = round(y * log2(x)),
245// f = y * log2(x) - k.
246//
247// In particular, k is an integer, and |f| <= 0.5.
248// For the final result to fit in single precision:
249// -150 <= y * log2(x) <= 128,
250// and anything outside that range overflows or underflows.
251//
252// Fast pass:
253// 1. Range reduction for log2(m_x) using a 32-entry lookup table:
254// dx = r * m_x - 1 in [-2^-6, 2^-5],
255// log2(m_x) = log2(1 + dx) - log2(r),
256// where -log2(r) is retrieved from LOG2_R_32[idx_x], and log2(1 + dx)/dx is
257// approximated by a degree-6 polynomial in double precision with Estrin's
258// scheme.
259// 2. Exponent reduction:
260// z = y * log2(x),
261// k = round(z),
262// f = z - k in [-0.5, 0.5].
263// 3. 2^f is evaluated using a degree-7 polynomial.
264// Absorbing +-ERR into the constant term 1.0 +- ERR of the final FMA allows
265// computing the upper and lower bounds for Ziv's rounding test in parallel.
266// 4. If upper == lower in single precision, return upper.
267// 5. If Ziv's test fails, fall back to check_exact_boundary and powf_accurate
268// (DoubleDouble / 128-entry table).
269//
270// Accurate pass:
271// 1. Range reduction using a 128-entry table: -log2(r) = {lo, mid, hi}.
272// log2(1 + dx) is evaluated using a degree-11 polynomial in DoubleDouble.
273// 2. Decompose: log2(x) = (e_x + hi) + log2_tail, where:
274// log2_tail = log2(1 + dx) + mid + lo.
275// Because e_x + hi has <= 28 bits of precision and y has 24 bits,
276// y * (e_x + hi) is exact in double precision.
277// 3. Exponent reduction with 64-entry EXP2_MID1 table and degree-8 polynomial
278// in DoubleDouble.
279LIBC_INLINE float powf(float x, float y) {
280 using FloatBits = fputil::FPBits<float>;
281 using DoubleBits = fputil::FPBits<double>;
282 using fputil::DoubleDouble;
283 using namespace common_constants_internal;
284
285 FloatBits xbits(x), ybits(y);
286 uint32_t x_u = xbits.uintval();
287 uint32_t y_u = ybits.uintval();
288 uint32_t y_a = ybits.abs().uintval();
289
290 // Quick filter for special inputs:
291 // - x: 0, +-1, 2^k, +-Inf, NaN.
292 // - y: 0, +-1, +-2, +-0.5, +-Inf, NaN.
293 if (LIBC_UNLIKELY((x_u & 0x001F'FFFF) == 0 || (y_u & 0x007F'FFFF) == 0)) {
294 if (auto r = powf_internal::check_special_inputs(x, y);
295 LIBC_UNLIKELY(r.has_value()))
296 return r.value();
297 }
298
299#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \
300 !defined(LIBC_MATH_HAS_SMALL_TABLES)
301 float orig_x_abs = xbits.abs().get_val();
302 float orig_y = y;
303#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS && !LIBC_MATH_HAS_SMALL_TABLES
304
305 int ex = -FloatBits::EXP_BIAS;
306 uint64_t sign = 0;
307
308 // Check for exceptional cases:
309 // - |y| <= 2^-40 or |y| >= 2^40.
310 // - x < 0, x is subnormal, zero, Inf, or NaN.
311 if (LIBC_UNLIKELY(y_a <= powf_internal::Y_LOWER_BOUND ||
312 y_a >= powf_internal::Y_UPPER_BOUND ||
313 x_u >= FloatBits::inf().uintval() ||
314 x_u < FloatBits::min_normal().uintval())) {
315 if (auto r = powf_internal::check_exceptional_cases(x, y, ex, sign);
316 LIBC_UNLIKELY(r.has_value()))
317 return r.value();
318 }
319
320 Sign out_sign = (sign == 0) ? Sign::POS : Sign::NEG;
321
322 // x^y = 2^( y * log2(x) ) = 2^( y * ( e_x + log2(m_x) ) )
323 // Compute log2(x) = e_x + log2(m_x)
324 x_u = FloatBits(x).uintval();
325
326 // Extract exponent field of x.
327 ex += (x_u >> FloatBits::FRACTION_LEN);
328 double e_x = static_cast<double>(ex);
329 uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
330 // Top 5 bits of mantissa for 32-entry table lookup.
331 int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 5));
332 // Embed mantissa into double format: m_x in [1.0, 2.0).
333 uint64_t m_x_u =
334 (static_cast<uint64_t>(x_mant) << 29) | 0x3ff0'0000'0000'0000ULL;
335 double m_x = cpp::bit_cast<double>(from: m_x_u);
336
337 // Range reduction: dx = m_x * r - 1.0 in [-2^-6, 2^-5].
338 // Since r has <= 7 bits and m_x has 24 bits, m_x * r - 1.0 is exact in
339 // double.
340 double dx =
341 fputil::multiply_add(x: m_x, y: powf_internal::R_32_D[idx_x], z: -1.0); // Exact
342
343 // Degree-6 minimax polynomial approximation for log2(1 + dx)/dx on
344 // [-2^-6, 2^-5]. Generated by Sollya with:
345 // > display = hexadecimal;
346 // > p6 = fpminimax(log2(1 + x)/x, 6, [|D...|], [-2^-6, 2^-5]);
347 // > for i from 0 to 6 do print(round(coeff(p6, i), D, RN), ",");
348 // > dirtyinfnorm((log2(1 + x) - x * p6) / log2(1 + x), [-2^-6, 2^-5]);
349 // 0x1.050a...p-47
350 // > dirtyinfnorm(log2(1 + x) - x * p6, [-2^-6, 2^-5]);
351 // 0x1.712b...p-52
352 constexpr double COEFFS[] = {
353 0x1.71547652b831fp0, -0x1.71547652b2ff5p-1, 0x1.ec709dbd02463p-2,
354 -0x1.7154787f73adep-2, 0x1.2777b7f57b29fp-2, -0x1.ec55ea50eaba7p-3,
355 0x1.92c78e7d5a931p-3,
356 };
357
358 // Evaluate P(dx) ~ log2(1 + dx)/dx with Estrin's scheme:
359 double dx2 = dx * dx;
360 double lp0 = fputil::multiply_add(x: dx, y: COEFFS[1], z: COEFFS[0]);
361 double lp1 = fputil::multiply_add(x: dx, y: COEFFS[3], z: COEFFS[2]);
362 double lp2 = fputil::multiply_add(x: dx, y: COEFFS[5], z: COEFFS[4]);
363
364 double dx4 = dx2 * dx2;
365 double lq0 = fputil::multiply_add(x: dx2, y: lp1, z: lp0);
366 double lq1 = fputil::multiply_add(x: dx2, y: COEFFS[6], z: lp2);
367
368 double p = fputil::multiply_add(x: dx4, y: lq1, z: lq0);
369
370 // log2(x) = e_x - log2(r) + dx * P(dx).
371 double log2_r_ex = powf_internal::LOG2_R_32[idx_x] + e_x;
372 double s = fputil::multiply_add(x: dx, y: p, z: log2_r_ex);
373
374 double y_d = static_cast<double>(y);
375 double z = y_d * s;
376
377 // y * log2(x) = k + f, with k integer and |f| <= 0.5.
378 double kd = fputil::nearest_integer(x: z);
379 int64_t k = static_cast<int64_t>(kd);
380 double f = fputil::multiply_add(x: y_d, y: s, z: -kd);
381
382 // Degree-7 polynomial approximation P(f) ~ (2^f - 1)/f on [-0.5, 0.5]
383 // Generated by Sollya with:
384 // > display = hexadecimal;
385 // > p = fpminimax((2^x - 1)/x, 7, [|D...|], [-0.5, 0.5]);
386 // > for i from 0 to 7 do print(round(coeff(p, i), D, RN), ",");
387 // > dirtyinfnorm(2^x - (1 + x * p), [-0.5, 0.5]);
388 // 0x1.0496...p-39
389 // > dirtyinfnorm((2^x - (1 + x * p)) / 2^x, [-0.5, 0.5]);
390 // 0x1.0496...p-39
391 constexpr double EXP2_COEFFS[] = {
392 0x1.62e42fef9cdf7p-1, 0x1.ebfbdff85f49p-3, 0x1.c6b08da69c963p-5,
393 0x1.3b2ab6a0e1172p-7, 0x1.5d87762dfc134p-10, 0x1.43099b4e829a9p-13,
394 0x1.00c080f7699bep-16, 0x1.62c0108a065d8p-20,
395 };
396
397 // Evaluate P(f) ~ (2^f - 1)/f on [-0.5, 0.5] with Estrin's scheme:
398 double f2 = f * f;
399 double f4 = f2 * f2;
400
401 double p0 = fputil::multiply_add(x: f, y: EXP2_COEFFS[1], z: EXP2_COEFFS[0]);
402 double p1 = fputil::multiply_add(x: f, y: EXP2_COEFFS[3], z: EXP2_COEFFS[2]);
403 double p2 = fputil::multiply_add(x: f, y: EXP2_COEFFS[5], z: EXP2_COEFFS[4]);
404 double p3 = fputil::multiply_add(x: f, y: EXP2_COEFFS[7], z: EXP2_COEFFS[6]);
405
406 double q0 = fputil::multiply_add(x: f2, y: p1, z: p0);
407 double q1 = fputil::multiply_add(x: f2, y: p3, z: p2);
408
409 double poly = fputil::multiply_add(x: f4, y: q1, z: q0);
410
411 // Normal range: -125 <= k <= 128.
412 // (k + 125) as unsigned <= 253 covers all normal outputs.
413 uint64_t k_u = static_cast<uint64_t>(k + 125);
414
415 if (LIBC_LIKELY(k_u <= 253)) {
416 // Scale by 2^k.
417 uint64_t exp2_k_i =
418 (static_cast<uint64_t>(static_cast<int>(k) + DoubleBits::EXP_BIAS)
419 << DoubleBits::FRACTION_LEN) |
420 sign;
421 double exp2_k = cpp::bit_cast<double>(from: exp2_k_i);
422
423#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
424 double pp = fputil::multiply_add(f, poly, 1.0);
425 double r_d = pp * exp2_k;
426 return static_cast<float>(r_d);
427#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
428#ifdef LIBC_TARGET_CPU_HAS_FMA
429 // Absorbing +-ERR into 1.0 computes both bounds using two parallel FMAs.
430 constexpr double ERR = 0x1.5p-39;
431 double pp_hi = fputil::multiply_add(f, poly, 1.0 + ERR);
432 double pp_lo = fputil::multiply_add(f, poly, 1.0 - ERR);
433 double r_d = pp_hi * exp2_k;
434#else // !LIBC_TARGET_CPU_HAS_FMA
435 // Without FMA, compute unbiased pp = f * poly + 1.0 directly.
436 double pp = fputil::multiply_add(x: f, y: poly, z: 1.0);
437 double r_d = pp * exp2_k;
438#endif // LIBC_TARGET_CPU_HAS_FMA
439
440 float res = static_cast<float>(r_d);
441
442 // Since 2^k is an exact power of 2, testing if the bounds round to the same
443 // float is equivalent to testing the fully scaled bounds.
444#ifdef LIBC_TARGET_CPU_HAS_FMA
445 float upper = static_cast<float>(pp_hi);
446 float lower = static_cast<float>(pp_lo);
447#else // !LIBC_TARGET_CPU_HAS_FMA
448 constexpr double ERR = 0x1.5p-39;
449 float upper = static_cast<float>(pp + ERR);
450 float lower = static_cast<float>(pp - ERR);
451#endif // LIBC_TARGET_CPU_HAS_FMA
452
453 if (LIBC_LIKELY(upper == lower))
454 return res;
455
456 // Accurate fallback when Ziv's test fails:
457#ifndef LIBC_MATH_HAS_SMALL_TABLES
458 if (auto r = check_exact_boundary(x: orig_x_abs, y: orig_y, sign, out_sign);
459 LIBC_UNLIKELY(r.has_value()))
460 return r.value();
461#endif // !LIBC_MATH_HAS_SMALL_TABLES
462
463 double r_dd = powf_accurate(x, y, e_x, sign);
464 res = static_cast<float>(r_dd);
465 if (LIBC_UNLIKELY(FloatBits(res).is_inf()))
466 return powf_internal::set_overflow(out_sign);
467 return res;
468#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
469 }
470
471#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
472 // Exceptional path: k > 128 (overflow), k < -155 (underflow), or denormal.
473 if (k > 128)
474 return powf_internal::set_overflow(out_sign);
475
476 if (k < -155)
477 return powf_internal::set_underflow(out_sign);
478
479 // Denormal and underflow path for -155 <= k <= -126:
480#ifndef LIBC_MATH_HAS_SMALL_TABLES
481 if (auto r = check_exact_boundary(x: orig_x_abs, y: orig_y, sign, out_sign);
482 LIBC_UNLIKELY(r.has_value()))
483 return r.value();
484#endif // !LIBC_MATH_HAS_SMALL_TABLES
485
486 int64_t exp2_k_i = (static_cast<uint64_t>(k + DoubleBits::EXP_BIAS)
487 << DoubleBits::FRACTION_LEN) |
488 sign;
489 double exp2_k = cpp::bit_cast<double>(from: exp2_k_i);
490
491 // Add +-2^(-126 + (53 - 24)) = +-2^-97 to mimic single precision denormal
492 // rounding in double precision.
493 // For -155 <= k <= -126, the least significant bit of the double precision
494 // sum (r_d + denorm_bias) is aligned at 2^(-97 - 52) = 2^-149, matching the
495 // least significant bit of single precision denormals 2^(-126 - 23) = 2^-149.
496 // This keeps the values in the double precision normal range while performing
497 // rounding according to the current rounding mode.
498 double denorm_bias = cpp::bit_cast<double>(from: 0x39e0'0000'0000'0000ULL | sign);
499
500 constexpr double ERR = 0x1.5p-39;
501
502#ifdef LIBC_TARGET_CPU_HAS_FMA
503 // Absorbing +-ERR into 1.0 computes both bounds using two parallel FMAs.
504 double pp_hi = fputil::multiply_add(f, poly, 1.0 + ERR);
505 double pp_lo = fputil::multiply_add(f, poly, 1.0 - ERR);
506 double u_hi = fputil::multiply_add(pp_hi, exp2_k, denorm_bias);
507 double u_lo = fputil::multiply_add(pp_lo, exp2_k, denorm_bias);
508#else // !LIBC_TARGET_CPU_HAS_FMA
509 // Without FMA, compute unbiased pp = f * poly + 1.0 directly.
510 double pp = fputil::multiply_add(x: f, y: poly, z: 1.0);
511 double r_d = pp * exp2_k;
512 double err = ERR * r_d;
513 double u_hi = (r_d + err) + denorm_bias;
514 double u_lo = (r_d - err) + denorm_bias;
515#endif // LIBC_TARGET_CPU_HAS_FMA
516
517 if (LIBC_LIKELY(u_hi == u_lo)) {
518 if (LIBC_UNLIKELY(u_hi == denorm_bias))
519 return powf_internal::set_underflow(out_sign);
520
521 float res = static_cast<float>(u_hi - denorm_bias);
522 if (LIBC_UNLIKELY(FloatBits(res).is_normal())) {
523#ifdef LIBC_TARGET_CPU_HAS_FMA
524 double pp = fputil::multiply_add(f, poly, 1.0);
525#endif // LIBC_TARGET_CPU_HAS_FMA
526 if (static_cast<float>(pp) < 1.0f)
527 fputil::raise_underflow_except_if_required<float>();
528 return res;
529 }
530
531 fputil::set_errno_if_required(ERANGE);
532 fputil::raise_underflow_except_if_required<float>();
533 return res;
534 }
535
536 // Ziv's test failed for denormal input, fall back to accurate pass.
537 double r_dd = powf_accurate(x, y, e_x, sign);
538 float res = static_cast<float>(r_dd);
539 if (LIBC_UNLIKELY(FloatBits(res).is_normal())) {
540 if (static_cast<float>(r_dd * 0x1.0p126) < 1.0f)
541 fputil::raise_underflow_except_if_required<float>();
542 return res;
543 }
544 fputil::set_errno_if_required(ERANGE);
545 fputil::raise_underflow_except_if_required<float>();
546 return res;
547#else // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
548 if (k > 128)
549 return powf_internal::set_overflow(out_sign);
550 return powf_internal::set_underflow(out_sign);
551#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
552}
553
554} // namespace double_eval
555} // namespace math
556} // namespace LIBC_NAMESPACE_DECL
557
558#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
559