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/// Common utilities and tables for double-precision pow(x, y).
11///
12//===----------------------------------------------------------------------===//
13
14#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POW_UTILS_H
15#define LLVM_LIBC_SRC___SUPPORT_MATH_POW_UTILS_H
16
17#include "hdr/errno_macros.h"
18#include "hdr/fenv_macros.h"
19#include "src/__support/CPP/bit.h"
20#include "src/__support/CPP/optional.h"
21#include "src/__support/FPUtil/FEnvImpl.h"
22#include "src/__support/FPUtil/FPBits.h"
23#include "src/__support/FPUtil/double_double.h"
24#include "src/__support/FPUtil/nearest_integer.h"
25#include "src/__support/FPUtil/rounding_mode.h"
26#include "src/__support/FPUtil/sqrt.h"
27#include "src/__support/common.h"
28#include "src/__support/macros/attributes.h"
29#include "src/__support/macros/config.h"
30#include "src/__support/macros/optimization.h"
31#include "src/__support/macros/properties/cpu_features.h"
32#include "src/__support/math/common_constants.h"
33#include "src/__support/math/exp10.h"
34#include "src/__support/math/exp2.h"
35
36namespace LIBC_NAMESPACE_DECL {
37namespace math {
38namespace pow_internal {
39
40using fputil::DoubleDouble;
41using namespace common_constants_internal;
42
43// Constants for log2(x) range reduction, generated by Sollya with:
44// > for i from 0 to 127 do {
45// r = 2^-8 * ceil( 2^8 * (1 - 2^(-8)) / (1 + i*2^-7) );
46// b = nearestint(log2(r) * 2^41) * 2^-41;
47// c = round(log2(r) - b, D, RN);
48// print("{ ", -c, ", ", -b, " }, ");
49// };
50// This is the same as -log2(RD[i]), with the least significant bits of the
51// high part set to be 2^-41, so that the sum of high parts + e_x is exact in
52// double precision.
53// We also replace the first and the last ones to be 0.
54LIBC_INLINE_VAR constexpr DoubleDouble LOG2_R_DD[128] = {
55 {.lo: 0.0, .hi: 0.0},
56 {.lo: -0x1.19b14945cf6bap-44, .hi: 0x1.72c7ba21p-7},
57 {.lo: -0x1.95539356f93dcp-43, .hi: 0x1.743ee862p-6},
58 {.lo: 0x1.abe0a48f83604p-43, .hi: 0x1.184b8e4c5p-5},
59 {.lo: 0x1.635577970e04p-43, .hi: 0x1.77394c9d9p-5},
60 {.lo: -0x1.401fbaaa67e3cp-45, .hi: 0x1.d6ebd1f2p-5},
61 {.lo: -0x1.5b1799ceaeb51p-43, .hi: 0x1.1bb32a6008p-4},
62 {.lo: 0x1.7c407050799bfp-43, .hi: 0x1.4c560fe688p-4},
63 {.lo: 0x1.da6339da288fcp-43, .hi: 0x1.7d60496cf8p-4},
64 {.lo: 0x1.be4f6f22dbbadp-43, .hi: 0x1.960caf9ab8p-4},
65 {.lo: -0x1.c760bc9b188c4p-45, .hi: 0x1.c7b528b71p-4},
66 {.lo: 0x1.164e932b2d51cp-44, .hi: 0x1.f9c95dc1dp-4},
67 {.lo: 0x1.924ae921f7ecap-45, .hi: 0x1.097e38ce6p-3},
68 {.lo: -0x1.6d25a5b8a19b2p-44, .hi: 0x1.22dadc2ab4p-3},
69 {.lo: 0x1.e50a1644ac794p-43, .hi: 0x1.3c6fb650ccp-3},
70 {.lo: 0x1.f34baa74a7942p-43, .hi: 0x1.494f863b8cp-3},
71 {.lo: -0x1.8f7aac147fdc1p-46, .hi: 0x1.633a8bf438p-3},
72 {.lo: 0x1.f84be19cb9578p-43, .hi: 0x1.7046031c78p-3},
73 {.lo: -0x1.66cccab240e9p-46, .hi: 0x1.8a8980abfcp-3},
74 {.lo: -0x1.3f7a55cd2af4cp-47, .hi: 0x1.97c1cb13c8p-3},
75 {.lo: 0x1.3458cde69308cp-43, .hi: 0x1.b2602497d4p-3},
76 {.lo: -0x1.667f21fa8423fp-44, .hi: 0x1.bfc67a8p-3},
77 {.lo: 0x1.d2fe4574e09b9p-47, .hi: 0x1.dac22d3e44p-3},
78 {.lo: 0x1.367bde40c5e6dp-43, .hi: 0x1.e857d3d36p-3},
79 {.lo: 0x1.d45da26510033p-46, .hi: 0x1.01d9bbcfa6p-2},
80 {.lo: -0x1.7204f55bbf90dp-44, .hi: 0x1.08bce0d96p-2},
81 {.lo: -0x1.d4f1b95e0ff45p-43, .hi: 0x1.169c05364p-2},
82 {.lo: 0x1.c20d74c0211bfp-44, .hi: 0x1.1d982c9d52p-2},
83 {.lo: 0x1.ad89a083e072ap-43, .hi: 0x1.249cd2b13cp-2},
84 {.lo: 0x1.cd0cb4492f1bcp-43, .hi: 0x1.32bfee370ep-2},
85 {.lo: -0x1.2101a9685c779p-47, .hi: 0x1.39de8e155ap-2},
86 {.lo: 0x1.9451cd394fe8dp-43, .hi: 0x1.4106017c3ep-2},
87 {.lo: 0x1.661e393a16b95p-44, .hi: 0x1.4f6fbb2cecp-2},
88 {.lo: -0x1.c6d8d86531d56p-44, .hi: 0x1.56b22e6b58p-2},
89 {.lo: 0x1.c1c885adb21d3p-43, .hi: 0x1.5dfdcf1eeap-2},
90 {.lo: 0x1.3bb5921006679p-45, .hi: 0x1.6552b49986p-2},
91 {.lo: 0x1.1d406db502403p-43, .hi: 0x1.6cb0f6865cp-2},
92 {.lo: 0x1.55a63e278bad5p-43, .hi: 0x1.7b89f02cf2p-2},
93 {.lo: -0x1.66ae2a7ada553p-49, .hi: 0x1.8304d90c12p-2},
94 {.lo: -0x1.66cccab240e9p-45, .hi: 0x1.8a8980abfcp-2},
95 {.lo: -0x1.62404772a151dp-45, .hi: 0x1.921800924ep-2},
96 {.lo: 0x1.ac9bca36fd02ep-44, .hi: 0x1.99b072a96cp-2},
97 {.lo: 0x1.4bc302ffa76fbp-43, .hi: 0x1.a8ff97181p-2},
98 {.lo: 0x1.01fea1ec47c71p-43, .hi: 0x1.b0b67f4f46p-2},
99 {.lo: -0x1.f20203b3186a6p-43, .hi: 0x1.b877c57b1cp-2},
100 {.lo: -0x1.2642415d47384p-45, .hi: 0x1.c043859e3p-2},
101 {.lo: -0x1.bc76a2753b99bp-50, .hi: 0x1.c819dc2d46p-2},
102 {.lo: -0x1.da93ae3a5f451p-43, .hi: 0x1.cffae611aep-2},
103 {.lo: -0x1.50e785694a8c6p-43, .hi: 0x1.d7e6c0abc4p-2},
104 {.lo: 0x1.c56138c894641p-43, .hi: 0x1.dfdd89d586p-2},
105 {.lo: 0x1.5669df6a2b592p-43, .hi: 0x1.e7df5fe538p-2},
106 {.lo: -0x1.ea92d9e0e8ac2p-48, .hi: 0x1.efec61b012p-2},
107 {.lo: 0x1.a0331af2e6feap-43, .hi: 0x1.f804ae8d0cp-2},
108 {.lo: 0x1.9518ce032f41dp-48, .hi: 0x1.0014332bep-1},
109 {.lo: -0x1.b3b3864c60011p-44, .hi: 0x1.042bd4b9a8p-1},
110 {.lo: -0x1.103e8f00d41c8p-45, .hi: 0x1.08494c66b9p-1},
111 {.lo: 0x1.65be75cc3da17p-43, .hi: 0x1.0c6caaf0c5p-1},
112 {.lo: 0x1.3676289cd3dd4p-43, .hi: 0x1.1096015deep-1},
113 {.lo: -0x1.41dfc7d7c3321p-43, .hi: 0x1.14c560fe69p-1},
114 {.lo: 0x1.e0cda8bd74461p-44, .hi: 0x1.18fadb6e2dp-1},
115 {.lo: 0x1.2a606046ad444p-44, .hi: 0x1.1d368296b5p-1},
116 {.lo: 0x1.f9ea977a639cp-43, .hi: 0x1.217868b0c3p-1},
117 {.lo: -0x1.50520a377c7ecp-45, .hi: 0x1.25c0a0463cp-1},
118 {.lo: 0x1.6e3cb71b554e7p-47, .hi: 0x1.2a0f3c3407p-1},
119 {.lo: -0x1.4275f1035e5e8p-48, .hi: 0x1.2e644fac05p-1},
120 {.lo: -0x1.4275f1035e5e8p-48, .hi: 0x1.2e644fac05p-1},
121 {.lo: -0x1.979a5db68721dp-45, .hi: 0x1.32bfee370fp-1},
122 {.lo: 0x1.1ee969a95f529p-43, .hi: 0x1.37222bb707p-1},
123 {.lo: 0x1.bb4b69336b66ep-43, .hi: 0x1.3b8b1c68fap-1},
124 {.lo: 0x1.d5e6a8a4fb059p-45, .hi: 0x1.3ffad4e74fp-1},
125 {.lo: 0x1.3106e404cabb7p-44, .hi: 0x1.44716a2c08p-1},
126 {.lo: 0x1.3106e404cabb7p-44, .hi: 0x1.44716a2c08p-1},
127 {.lo: -0x1.9bcaf1aa4168ap-43, .hi: 0x1.48eef19318p-1},
128 {.lo: 0x1.1646b761c48dep-44, .hi: 0x1.4d7380dcc4p-1},
129 {.lo: 0x1.2f0c0bfe9dbecp-43, .hi: 0x1.51ff2e3021p-1},
130 {.lo: 0x1.29904613e33cp-43, .hi: 0x1.5692101d9bp-1},
131 {.lo: 0x1.1d406db502403p-44, .hi: 0x1.5b2c3da197p-1},
132 {.lo: 0x1.1d406db502403p-44, .hi: 0x1.5b2c3da197p-1},
133 {.lo: -0x1.125d6cbcd1095p-44, .hi: 0x1.5fcdce2728p-1},
134 {.lo: -0x1.bd9b32266d92cp-43, .hi: 0x1.6476d98adap-1},
135 {.lo: 0x1.54243b21709cep-44, .hi: 0x1.6927781d93p-1},
136 {.lo: 0x1.54243b21709cep-44, .hi: 0x1.6927781d93p-1},
137 {.lo: -0x1.ce60916e52e91p-44, .hi: 0x1.6ddfc2a79p-1},
138 {.lo: 0x1.f1f5ae718f241p-43, .hi: 0x1.729fd26b7p-1},
139 {.lo: -0x1.6eb9612e0b4f3p-43, .hi: 0x1.7767c12968p-1},
140 {.lo: -0x1.6eb9612e0b4f3p-43, .hi: 0x1.7767c12968p-1},
141 {.lo: 0x1.fed21f9cb2cc5p-43, .hi: 0x1.7c37a9227ep-1},
142 {.lo: 0x1.7f5dc57266758p-43, .hi: 0x1.810fa51bf6p-1},
143 {.lo: 0x1.7f5dc57266758p-43, .hi: 0x1.810fa51bf6p-1},
144 {.lo: 0x1.5b338360c2ae2p-43, .hi: 0x1.85efd062c6p-1},
145 {.lo: -0x1.96fc8f4b56502p-43, .hi: 0x1.8ad846cf37p-1},
146 {.lo: -0x1.96fc8f4b56502p-43, .hi: 0x1.8ad846cf37p-1},
147 {.lo: -0x1.bdc81c4db3134p-44, .hi: 0x1.8fc924c89bp-1},
148 {.lo: 0x1.36c101ee1344p-43, .hi: 0x1.94c287492cp-1},
149 {.lo: 0x1.36c101ee1344p-43, .hi: 0x1.94c287492cp-1},
150 {.lo: 0x1.e41fa0a62e6aep-44, .hi: 0x1.99c48be206p-1},
151 {.lo: -0x1.d97ee9124773bp-46, .hi: 0x1.9ecf50bf44p-1},
152 {.lo: -0x1.d97ee9124773bp-46, .hi: 0x1.9ecf50bf44p-1},
153 {.lo: -0x1.3f94e00e7d6bcp-46, .hi: 0x1.a3e2f4ac44p-1},
154 {.lo: -0x1.6879fa00b120ap-43, .hi: 0x1.a8ff971811p-1},
155 {.lo: -0x1.6879fa00b120ap-43, .hi: 0x1.a8ff971811p-1},
156 {.lo: 0x1.1659d8e2d7d38p-44, .hi: 0x1.ae255819fp-1},
157 {.lo: 0x1.1e5e0ae0d3f8ap-43, .hi: 0x1.b35458761dp-1},
158 {.lo: 0x1.1e5e0ae0d3f8ap-43, .hi: 0x1.b35458761dp-1},
159 {.lo: 0x1.484a15babcf88p-43, .hi: 0x1.b88cb9a2abp-1},
160 {.lo: 0x1.484a15babcf88p-43, .hi: 0x1.b88cb9a2abp-1},
161 {.lo: 0x1.871a7610e40bdp-45, .hi: 0x1.bdce9dcc96p-1},
162 {.lo: -0x1.2d90e5edaeceep-43, .hi: 0x1.c31a27dd01p-1},
163 {.lo: -0x1.2d90e5edaeceep-43, .hi: 0x1.c31a27dd01p-1},
164 {.lo: -0x1.5dd31d962d373p-43, .hi: 0x1.c86f7b7ea5p-1},
165 {.lo: -0x1.5dd31d962d373p-43, .hi: 0x1.c86f7b7ea5p-1},
166 {.lo: -0x1.9ad57391924a7p-43, .hi: 0x1.cdcebd2374p-1},
167 {.lo: -0x1.3167ccc538261p-44, .hi: 0x1.d338120a6ep-1},
168 {.lo: -0x1.3167ccc538261p-44, .hi: 0x1.d338120a6ep-1},
169 {.lo: 0x1.c7a4ff65ddbc9p-45, .hi: 0x1.d8aba045bp-1},
170 {.lo: 0x1.c7a4ff65ddbc9p-45, .hi: 0x1.d8aba045bp-1},
171 {.lo: -0x1.f9ab3cf74babap-44, .hi: 0x1.de298ec0bbp-1},
172 {.lo: -0x1.f9ab3cf74babap-44, .hi: 0x1.de298ec0bbp-1},
173 {.lo: 0x1.52842c1c1e586p-43, .hi: 0x1.e3b20546f5p-1},
174 {.lo: 0x1.52842c1c1e586p-43, .hi: 0x1.e3b20546f5p-1},
175 {.lo: 0x1.3c6764fc87b4ap-48, .hi: 0x1.e9452c8a71p-1},
176 {.lo: 0x1.3c6764fc87b4ap-48, .hi: 0x1.e9452c8a71p-1},
177 {.lo: -0x1.a0976c0a2827dp-44, .hi: 0x1.eee32e2aedp-1},
178 {.lo: -0x1.a0976c0a2827dp-44, .hi: 0x1.eee32e2aedp-1},
179 {.lo: -0x1.a45314dc4fc42p-43, .hi: 0x1.f48c34bd1fp-1},
180 {.lo: -0x1.a45314dc4fc42p-43, .hi: 0x1.f48c34bd1fp-1},
181 {.lo: 0x1.ef5d00e390ap-44, .hi: 0x1.fa406bd244p-1},
182 {.lo: 0.0, .hi: 1.0},
183};
184
185// Check if x is an odd integer: the lowest set bit must be at the unit
186// position:
187// x_e + lsb == UNIT_EXPONENT.
188LIBC_INLINE bool is_odd_integer(double x) {
189 using FPBits = fputil::FPBits<double>;
190 FPBits xbits(x);
191 uint64_t x_u = xbits.uintval();
192 unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
193 unsigned lsb =
194 static_cast<unsigned>(cpp::countr_zero(value: x_u | FPBits::EXP_MASK));
195 constexpr unsigned UNIT_EXPONENT =
196 static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
197 return (x_e + lsb == UNIT_EXPONENT);
198}
199
200// Check if x is an integer: the lowest set bit must be at or above the unit
201// position:
202// x_e + lsb >= UNIT_EXPONENT.
203LIBC_INLINE bool is_integer(double x) {
204 if (x == 0.0)
205 return true;
206 using FPBits = fputil::FPBits<double>;
207 FPBits xbits(x);
208 uint64_t x_u = xbits.uintval();
209 unsigned x_e = static_cast<unsigned>(xbits.get_biased_exponent());
210 unsigned lsb =
211 static_cast<unsigned>(cpp::countr_zero(value: x_u | FPBits::EXP_MASK));
212 constexpr unsigned UNIT_EXPONENT =
213 static_cast<unsigned>(FPBits::EXP_BIAS + FPBits::FRACTION_LEN);
214 return (x_e + lsb >= UNIT_EXPONENT);
215}
216
217LIBC_INLINE double set_overflow(bool is_neg) {
218 fputil::set_errno_if_required(ERANGE);
219 fputil::raise_overflow_except_if_required<double>();
220 using FPBits = fputil::FPBits<double>;
221#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
222 int rounding = fputil::quick_get_round();
223 if (rounding == FE_TOWARDZERO)
224 return is_neg ? -FPBits::max_normal().get_val()
225 : FPBits::max_normal().get_val();
226 if (rounding == FE_DOWNWARD)
227 return is_neg ? -FPBits::inf().get_val() : FPBits::max_normal().get_val();
228 if (rounding == FE_UPWARD)
229 return is_neg ? -FPBits::max_normal().get_val() : FPBits::inf().get_val();
230#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
231 return is_neg ? -FPBits::inf().get_val() : FPBits::inf().get_val();
232}
233
234LIBC_INLINE double set_underflow(bool is_neg) {
235 fputil::set_errno_if_required(ERANGE);
236 fputil::raise_underflow_except_if_required<double>();
237#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
238 using FPBits = fputil::FPBits<double>;
239 int rounding = fputil::quick_get_round();
240 if (rounding == FE_UPWARD && !is_neg)
241 return FPBits::min_subnormal().get_val();
242 if (rounding == FE_DOWNWARD && is_neg)
243 return -FPBits::min_subnormal().get_val();
244#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
245 return is_neg ? -0.0 : 0.0;
246}
247
248// Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We
249// assume further that
250// 1 <= mid < 2,
251// mid + lo < 2,
252// and |lo| << mid.
253//
254// Notice that, if 0 < x < 2^-1022:
255// double( 2^-1022 + x ) - 2^-1022 = double(x),
256// and if -2^-1022 < x < 0:
257// double(-2^-1022 + x ) - -2^-1022 = double(x).
258//
259// So if we scale |x| up by 2^1022, we can use the "+1 trick" (or "-1 trick"
260// when x < 0):
261// double( 1.0 + 2^1022 * x) - 1.0, for x > 0, and
262// double(-1.0 + 2^1022 * x) - -1.0, for x < 0,
263// to test how x is rounded in the denormal range. By adding +-1.0, the
264// results will have similar rounding points as denormal outputs, aligning the
265// least significant bit of the mantissa with 2^-52 * 2^-1022 = 2^-1074.
266//
267// Finally, subtracting 1023 in the exponent field converts the rounded normal
268// value into the subnormal value.
269template <bool SKIP_ZIV_TEST = false>
270LIBC_INLINE constexpr cpp::optional<double>
271ziv_test_denorm(int hi, double mid, double lo, double err, bool is_neg) {
272 using FPBits = typename fputil::FPBits<double>;
273
274 uint64_t scale_u = static_cast<uint64_t>(hi + 2045) << FPBits::FRACTION_LEN;
275 double exp_scale = cpp::bit_cast<double>(from: scale_u);
276
277 double mid_hi = mid * exp_scale;
278 double lo_scaled = lo * exp_scale;
279
280 double extra_factor = 0.0;
281 uint64_t scale_down = 0x3FE0'0000'0000'0000ULL; // 1022 in the exponent field.
282
283 // Result is denormal if |mid_hi + lo_scaled| < 1.0.
284 // In that case, add +-1.0 to mimic denormal rounding points.
285 if (is_neg) {
286 if ((-1.0 - mid_hi) < lo_scaled) {
287 // Extra rounding step is needed, which adds more rounding errors.
288 extra_factor = -1.0;
289 scale_down = 0x3FF0'0000'0000'0000ULL; // 1023 in the exponent field.
290 err += 0x1.0p-52;
291 }
292 } else {
293 if ((1.0 - mid_hi) > lo_scaled) {
294 // Extra rounding step is needed, which adds more rounding errors.
295 extra_factor = 1.0;
296 scale_down = 0x3FF0'0000'0000'0000ULL; // 1023 in the exponent field.
297 err += 0x1.0p-52;
298 }
299 }
300
301 // By adding +-1.0, the hardware adder performs rounding at the exact same
302 // bit position as denormal rounding.
303 if constexpr (SKIP_ZIV_TEST) {
304 double r = extra_factor + (mid_hi + lo_scaled);
305 return cpp::bit_cast<double>(from: cpp::bit_cast<uint64_t>(from: r) - scale_down);
306 } else {
307 double err_scaled = err * exp_scale;
308
309 double lo_u = lo_scaled + err_scaled;
310 double lo_l = lo_scaled - err_scaled;
311
312 double upper = extra_factor + (mid_hi + lo_u);
313 double lower = extra_factor + (mid_hi + lo_l);
314
315 if (LIBC_LIKELY(upper == lower)) {
316 return cpp::bit_cast<double>(from: cpp::bit_cast<uint64_t>(from: upper) - scale_down);
317 }
318
319 return cpp::nullopt;
320 }
321}
322
323// Upper bound for y = 1075 / |log2(1 - 2^(-53)|, generated by Sollya:
324// > y = round(-1075 / log2(1 - 2^(-53)), D, RU);
325// > y;
326// 0x1.74910d52d3052p62
327// > printdouble(y);
328// 0x43d74910d52d3052
329constexpr uint64_t Y_UPPER_BOUND = 0x43d7'4910'd52d'3052;
330// Lower bound for y = 2^(-54) / 1074), generated by Sollya:
331// > y = round(2^(-54) / 1074), D, RD);
332// > y;
333// 0x1.e829f39aef509p-65
334// > printdouble(y);
335// 0x3bee829f39ae'f509
336constexpr uint64_t Y_LOWER_BOUND = 0x3bee'829f'39ae'f509;
337
338// Fast checks for special inputs:
339// x = 0, +-1, 2^k, 10, +- inf
340// y = 0, +-1, 2, 0.5
341LIBC_ALWAYS_INLINE cpp::optional<double> check_special_inputs(double x,
342 double y) {
343 using FPBits = fputil::FPBits<double>;
344 FPBits xbits(x), ybits(y);
345
346 bool x_sign = xbits.sign() == Sign::NEG;
347 bool y_sign = ybits.sign() == Sign::NEG;
348
349 FPBits x_abs = xbits.abs();
350 FPBits y_abs = ybits.abs();
351
352 uint64_t x_u = xbits.uintval();
353 uint64_t y_u = ybits.uintval();
354 uint64_t y_a = y_abs.uintval();
355
356 if (LIBC_UNLIKELY((x_u & 0x0003'FFFF'FFFF'FFFF) == 0) ||
357 ((y_u & 0x000F'FFFF'FFFF'FFFF) == 0)) {
358 // If x or y is signaling NaN
359 if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
360 fputil::raise_except_if_required(FE_INVALID);
361 return FPBits::quiet_nan().get_val();
362 }
363
364 if (x == 1.0 || y == 0.0)
365 return 1.0;
366
367 if (x == 0.0) {
368 if (y_abs.is_nan())
369 return y;
370 if (y_abs.is_inf())
371 return y_sign ? FPBits::inf().get_val() : 0.0;
372 bool out_is_neg = x_sign && is_odd_integer(x: y);
373 if (y_sign) {
374 // pow(0, negative number) = inf
375 fputil::set_errno_if_required(EDOM);
376 fputil::raise_except_if_required(FE_DIVBYZERO);
377 return FPBits::inf(sign: out_is_neg ? Sign::NEG : Sign::POS).get_val();
378 }
379 // pow(0, positive number) = 0
380 return out_is_neg ? -0.0 : 0.0;
381 }
382
383 if (y == 1.0)
384 return x;
385
386 if (y == 2.0)
387 return x * x;
388
389 if (y == 0.5 && !x_sign)
390 return fputil::sqrt<double>(x);
391
392 // TODO: Add special case y = -0.5 when rsqrt is available.
393
394 if (x == 2.0)
395 return math::exp2(x: y);
396
397 if (x == 10.0)
398 return math::exp10(x: y);
399
400 // For x = 2^(+- 2^n):
401 // x^y = 2^(+- 2^n * y).
402 if (!x_sign && x_abs.is_normal() && xbits.get_mantissa() == 0 &&
403 y_a > Y_LOWER_BOUND && y_a < Y_UPPER_BOUND) {
404 int e_x = xbits.get_exponent();
405 uint32_t abs_e = static_cast<uint32_t>(e_x > 0 ? e_x : -e_x);
406 if (cpp::has_single_bit(value: abs_e)) {
407 double hi = static_cast<double>(e_x) * y;
408 if (hi >= 1024.0)
409 return set_overflow(false);
410 if (hi < -1075.0)
411 return set_underflow(false);
412 if (hi >= -1022.0)
413 return math::exp2(x: hi);
414 }
415 }
416 }
417
418 return cpp::nullopt;
419}
420
421// Filters out extreme input ranges, infinities, NaNs, normalizes denormal
422// inputs, and handles negative bases. Returns a value if the result is
423// determined, or nullopt if regular evaluation should proceed (in which case x,
424// y, e_x, x_mant, is_neg, sign_d may be updated).
425LIBC_ALWAYS_INLINE cpp::optional<double>
426check_exceptional_cases(double &x, double &y, double &e_x, uint64_t &x_mant,
427 bool &is_neg, double &sign_d) {
428 using FPBits = fputil::FPBits<double>;
429
430 FPBits xbits(x), ybits(y);
431 bool x_sign = xbits.sign() == Sign::NEG;
432 bool y_sign = ybits.sign() == Sign::NEG;
433
434 FPBits x_abs = xbits.abs();
435 FPBits y_abs = ybits.abs();
436
437 uint64_t x_u = xbits.uintval();
438 uint64_t y_mant = ybits.get_mantissa();
439 uint64_t x_a = x_abs.uintval();
440 uint64_t y_a = y_abs.uintval();
441
442 if (LIBC_UNLIKELY(y_a <= Y_LOWER_BOUND || y_a >= Y_UPPER_BOUND ||
443 x_u >= FPBits::inf().uintval() ||
444 x_u < FPBits::min_normal().uintval())) {
445 // If x or y is signaling NaN
446 if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
447 fputil::raise_except_if_required(FE_INVALID);
448 return FPBits::quiet_nan().get_val();
449 }
450
451 // |y| > |1075 / log2(1 - 2^-53)|.
452 if (y_a > 0x43d7'4910'd52d'3052) {
453 if (y_a >= 0x7ff0'0000'0000'0000) {
454 // y is inf or nan
455 if (y_mant != 0) {
456 // y is NaN
457 // pow(1, NaN) = 1 is already dealt with.
458 // pow(x, NaN) = NaN
459 return y;
460 }
461
462 // Now y is +-Inf
463 if (x_abs.is_nan()) {
464 // pow(NaN, +-Inf) = NaN
465 return x;
466 }
467
468 if (x_a == 0x3ff0'0000'0000'0000) {
469 // pow(+-1, +-Inf) = 1.0
470 return 1.0;
471 }
472
473 // pow (|x| < 1, -inf) = +inf
474 // pow (|x| < 1, +inf) = 0.0
475 // pow (|x| > 1, -inf) = 0.0
476 // pow (|x| > 1, +inf) = +inf
477 return ((x_a < FPBits::one().uintval()) == y_sign)
478 ? FPBits::inf().get_val()
479 : 0.0;
480 }
481 // x^y will overflow / underflow in double precision. Set y to a
482 // large enough exponent but not too large, so that the computations
483 // won't overflow in double precision.
484 y = y_sign ? -0x1.0p100 : 0x1.0p100;
485 }
486
487 // y is finite and non-zero.
488
489 if (x_a == FPBits::inf().uintval()) {
490 bool out_is_neg = x_sign && is_odd_integer(x: y);
491 Sign out_sign = out_is_neg ? Sign::NEG : Sign::POS;
492 return y_sign ? FPBits::zero(sign: out_sign).get_val()
493 : FPBits::inf(sign: out_sign).get_val();
494 }
495
496 if (x_a > FPBits::inf().uintval()) {
497 // x is NaN.
498 // pow (aNaN, 0) is already taken care above.
499 return x;
500 }
501
502 // Normalize denormal inputs.
503 if (x_a < FPBits::min_normal().uintval()) {
504 FPBits x_norm(x * 0x1.0p64);
505 e_x = static_cast<double>(x_norm.get_exponent()) - 64.0;
506 x_mant = x_norm.get_mantissa();
507 }
508
509 // x is finite and negative, and y is a finite integer.
510 if (x_sign) {
511 if (is_integer(x: y)) {
512 x = -x;
513 if (is_odd_integer(x: y)) {
514 is_neg = true;
515 sign_d = -1.0;
516 }
517 } else {
518 // pow( negative, non-integer ) = NaN
519 fputil::set_errno_if_required(EDOM);
520 fputil::raise_except_if_required(FE_INVALID);
521 return FPBits::quiet_nan().get_val();
522 }
523 }
524
525 if (y_a <= Y_LOWER_BOUND) {
526 volatile double one = 1.0;
527 volatile double eps = ((x_a < FPBits::one().uintval()) == !y_sign)
528 ? -0x1.0p-100
529 : 0x1.0p-100;
530 return one + eps;
531 }
532 }
533
534 return cpp::nullopt;
535}
536
537// Check if x^y is an exact rounding boundary case (a 54-bit dyadic float,
538// which is either an exact 53-bit float or an exact halfway midpoint).
539//
540// Reference:
541// Lauter, C. and Lefevre, V., "Rounding Boundary Cases for Values of the
542// Exponential and Power Functions," IEEE Trans. Comput. 58(8):1063-1074.
543//
544// For a normalized 128-bit dyadic float with MSB at bit 127, the top 54 bits
545// occupy bits [127:74]. For any 54-bit number, bits [73:0] are zero.
546//
547// If x^y is not in F_54:
548// |x^y - o_54(x^y)| / x^y >= 2^-114.
549// With an evaluation error < 2^-115, the distance from bits [73:0] to the
550// nearest 54-bit boundary must satisfy:
551// dist >= (2^-114 - 2^-115) * 2^127 = 2^12.
552// We use a threshold of 2^-114 (dist < 2^13) to capture all boundary
553// candidates:
554//
555// 1. x = 2^e: x^y = 2^(e * y) is in F_54 iff e * y is an integer.
556// 2. x = m * 2^e, y = n * 2^f with m, n odd integers:
557// x^y in F_54 implies:
558// - 0 <= y <= 35, n <= 35, and f >= -5.
559// - e * y is an integer.
560// - When f < 0 (i.e. y = n / 2^(-f)): m^(2^f) is an integer, verified
561// by taking repeated square roots -f times.
562// - The resulting significand (m^(2^f))^n <= 2^54.
563//
564// Returns true and sets exact_m and exact_exp if x^y is an exact boundary.
565LIBC_INLINE bool is_exact_rounding_boundary(double x, double y,
566 uint64_t &exact_m, int &exact_exp) {
567 using FPBits = fputil::FPBits<double>;
568
569 FPBits xbits(x);
570 int x_e = 0;
571 uint64_t x_mant = xbits.get_mantissa();
572 if (LIBC_UNLIKELY(xbits.get_biased_exponent() == 0)) {
573 if (x_mant != 0) {
574 int shift = cpp::countl_zero(value: x_mant) - 11;
575 x_mant = (x_mant << shift) & FPBits::FRACTION_MASK;
576 x_e = -1022 - shift;
577 }
578 } else {
579 x_e = xbits.get_exponent();
580 }
581
582 // Case 1: x = 2^e.
583 if (x_mant == 0) {
584 double e = static_cast<double>(x_e);
585 double ey = e * y;
586 if (is_integer(x: ey)) {
587 exact_m = 1;
588 exact_exp = static_cast<int>(ey);
589 return true;
590 }
591 return false;
592 }
593
594 // Case 2: x is not a power of 2.
595 if (y < 0.0 || y > 35.0)
596 return false;
597
598 // Decompose y = n * 2^f with n odd.
599 FPBits ybits(y);
600 uint64_t y_mant = ybits.get_mantissa() | (1ULL << FPBits::FRACTION_LEN);
601 int y_exp = ybits.get_exponent() - static_cast<int>(FPBits::FRACTION_LEN);
602 int tz_y = cpp::countr_zero(value: y_mant);
603 uint64_t n = y_mant >> tz_y;
604 int f = y_exp + tz_y;
605
606 if (n > 35 || f < -5)
607 return false;
608
609 // Decompose x = m * 2^e with m odd.
610 uint64_t full_x_mant = x_mant | (1ULL << FPBits::FRACTION_LEN);
611 int tz_x = cpp::countr_zero(value: full_x_mant);
612 uint64_t m = full_x_mant >> tz_x;
613 int e = x_e - static_cast<int>(FPBits::FRACTION_LEN) + tz_x;
614
615 if (f < 0) {
616 // Non-integer exponent: y = n / 2^(-f) with -f in {1, ..., 5}.
617 double ey = static_cast<double>(e) * y;
618 if (!is_integer(x: ey))
619 return false;
620
621 // Check if m^(2^f) is an integer by taking repeated square roots -f times.
622 int count = -f;
623 uint64_t cur = m;
624 for (int i = 0; i < count; ++i) {
625 uint64_t s =
626 static_cast<uint64_t>(fputil::sqrt<double>(x: static_cast<double>(cur)));
627 if (s * s != cur)
628 return false;
629 cur = s;
630 }
631
632 // Compute res = cur^n and ensure res <= 2^54.
633 uint64_t res = 1;
634 for (uint64_t i = 0; i < n; ++i) {
635 if (res > (1ULL << 54) / cur)
636 return false;
637 res *= cur;
638 }
639 exact_m = res;
640 exact_exp = static_cast<int>(ey);
641 return true;
642 } else {
643 // Integer exponent: f >= 0, so y is an integer.
644 int y_int = static_cast<int>(y);
645 uint64_t res = 1;
646 for (int i = 0; i < y_int; ++i) {
647 if (res > (1ULL << 54) / m)
648 return false;
649 res *= m;
650 }
651 exact_m = res;
652 exact_exp = e * y_int;
653 return true;
654 }
655}
656
657} // namespace pow_internal
658} // namespace math
659} // namespace LIBC_NAMESPACE_DECL
660
661#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POW_UTILS_H
662