1//===-- Implementation header of sqrtf128 ---------------------------------===//
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#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SQRTF128_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_SQRTF128_H
11
12#include "include/llvm-libc-types/float128.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT128
15
16#include "src/__support/CPP/bit.h"
17#include "src/__support/FPUtil/FEnvImpl.h"
18#include "src/__support/FPUtil/FPBits.h"
19#include "src/__support/FPUtil/rounding_mode.h"
20#include "src/__support/common.h"
21#include "src/__support/macros/optimization.h"
22#include "src/__support/uint128.h"
23
24// Compute sqrtf128 with correct rounding for all rounding modes using integer
25// arithmetic by Alexei Sibidanov (sibid@uvic.ca):
26// https://github.com/sibidanov/llvm-project/tree/as_sqrt_v2
27// https://github.com/sibidanov/llvm-project/tree/as_sqrt_v3
28// TODO: Update the reference once Alexei's implementation is in the CORE-MATH
29// project. https://github.com/llvm/llvm-project/issues/126794
30
31// Let the input be expressed as x = 2^e * m_x,
32// - Step 1: Range reduction
33// Let x_reduced = 2^(e % 2) * m_x,
34// Then sqrt(x) = 2^(e / 2) * sqrt(x_reduced), with
35// 1 <= x_reduced < 4.
36// - Step 2: Polynomial approximation
37// Approximate 1/sqrt(x_reduced) using polynomial approximation with the
38// result errors bounded by:
39// |r0 - 1/sqrt(x_reduced)| < 2^-32.
40// The computations are done in uint64_t.
41// - Step 3: First Newton iteration
42// Let the scaled error defined by:
43// h0 = r0^2 * x_reduced - 1.
44// Then we compute the first Newton iteration:
45// r1 = r0 - r0 * h0 / 2.
46// The result is then bounded by:
47// |r1 - 1 / sqrt(x_reduced)| < 2^-62.
48// - Step 4: Second Newton iteration
49// We calculate the scaled error from Step 3:
50// h1 = r1^2 * x_reduced - 1.
51// Then the second Newton iteration is computed by:
52// r2 = x_reduced * (r1 - r1 * h0 / 2)
53// ~ x_reduced * (1/sqrt(x_reduced)) = sqrt(x_reduced)
54// - Step 5: Perform rounding test and correction if needed.
55// Rounding correction is done by computing the exact rounding errors:
56// x_reduced - r2^2.
57
58namespace LIBC_NAMESPACE_DECL {
59namespace math {
60
61namespace sqrtf128_internal {
62
63template <typename T, typename U = T> LIBC_INLINE constexpr T prod_hi(T, U);
64
65// Get high part of integer multiplications.
66// Use template to prevent implicit conversion.
67template <>
68LIBC_INLINE constexpr uint64_t prod_hi<uint64_t>(uint64_t x, uint64_t y) {
69 return static_cast<uint64_t>(
70 (static_cast<UInt128>(x) * static_cast<UInt128>(y)) >> 64);
71}
72
73// Get high part of unsigned 128x64 bit multiplication.
74template <>
75LIBC_INLINE constexpr UInt128 prod_hi<UInt128, uint64_t>(UInt128 x,
76 uint64_t y) {
77 uint64_t x_lo = static_cast<uint64_t>(x);
78 uint64_t x_hi = static_cast<uint64_t>(x >> 64);
79 UInt128 xyl = static_cast<UInt128>(x_lo) * static_cast<UInt128>(y);
80 UInt128 xyh = static_cast<UInt128>(x_hi) * static_cast<UInt128>(y);
81 return xyh + (xyl >> 64);
82}
83
84// Get high part of signed 64x64 bit multiplication.
85template <>
86LIBC_INLINE constexpr int64_t prod_hi<int64_t>(int64_t x, int64_t y) {
87 return static_cast<int64_t>(
88 (static_cast<Int128>(x) * static_cast<Int128>(y)) >> 64);
89}
90
91// Get high 128-bit part of unsigned 128x128 bit multiplication.
92template <>
93LIBC_INLINE constexpr UInt128 prod_hi<UInt128>(UInt128 x, UInt128 y) {
94 uint64_t x_lo = static_cast<uint64_t>(x);
95 uint64_t x_hi = static_cast<uint64_t>(x >> 64);
96 uint64_t y_lo = static_cast<uint64_t>(y);
97 uint64_t y_hi = static_cast<uint64_t>(y >> 64);
98
99 UInt128 xh_yh = static_cast<UInt128>(x_hi) * static_cast<UInt128>(y_hi);
100 UInt128 xh_yl = static_cast<UInt128>(x_hi) * static_cast<UInt128>(y_lo);
101 UInt128 xl_yh = static_cast<UInt128>(x_lo) * static_cast<UInt128>(y_hi);
102
103 xh_yh += xh_yl >> 64;
104
105 return xh_yh + (xl_yh >> 64);
106}
107
108// Get high 128-bit part of mixed sign 128x128 bit multiplication.
109template <>
110LIBC_INLINE constexpr Int128 prod_hi<Int128, UInt128>(Int128 x, UInt128 y) {
111 UInt128 mask = static_cast<UInt128>(x >> 127);
112 UInt128 negative_part = y & mask;
113 UInt128 prod = prod_hi(x: static_cast<UInt128>(x), y);
114 return static_cast<Int128>(prod - negative_part);
115}
116
117// Newton-Raphson first order step to improve accuracy of the result.
118// For the initial approximation r0 ~ 1/sqrt(x), let
119// h = r0^2 * x - 1
120// be its scaled error. Then the first-order Newton-Raphson iteration is:
121// r1 = r0 - r0 * h / 2
122// which has error bounded by:
123// |r1 - 1/sqrt(x)| < h^2 / 2.
124LIBC_INLINE constexpr uint64_t rsqrt_newton_raphson(uint64_t m, uint64_t r) {
125 uint64_t r2 = prod_hi(x: r, y: r);
126 // h = r0^2*x - 1.
127 int64_t h = static_cast<int64_t>(prod_hi(x: m, y: r2) + r2);
128 // hr = r * h / 2
129 int64_t hr = prod_hi(x: h, y: static_cast<int64_t>(r >> 1));
130 return r - hr;
131}
132
133#ifdef LIBC_MATH_HAS_SMALL_TABLES
134// Degree-12 minimax polynomials for 1/sqrt(x) on [1, 2].
135LIBC_INLINE_VAR constexpr uint32_t RSQRT_COEFFS[12] = {
136 0xb5947a4a, 0x2d651e32, 0x9ad50532, 0x2d28d093, 0x0d8be653, 0x04239014,
137 0x01492449, 0x0066ff7d, 0x001e74a1, 0x000984cc, 0x00049abc, 0x00018340,
138};
139
140LIBC_INLINE constexpr uint64_t rsqrt_approx(uint64_t m) {
141 int64_t x = static_cast<uint64_t>(m) ^ (uint64_t(1) << 63);
142 int64_t x_26 = x >> 2;
143 int64_t z = x >> 31;
144
145 if (LIBC_UNLIKELY(z <= -4294967296))
146 return ~(m >> 1);
147
148 uint64_t x2 = static_cast<uint64_t>(z) * static_cast<uint64_t>(z);
149 uint64_t x2_26 = x2 >> 5;
150 x2 >>= 32;
151 // Calculate the odd part of the polynomial using Horner's method.
152 uint64_t c0 = RSQRT_COEFFS[8] + ((x2 * RSQRT_COEFFS[10]) >> 32);
153 uint64_t c1 = RSQRT_COEFFS[6] + ((x2 * c0) >> 32);
154 uint64_t c2 = RSQRT_COEFFS[4] + ((x2 * c1) >> 32);
155 uint64_t c3 = RSQRT_COEFFS[2] + ((x2 * c2) >> 32);
156 uint64_t c4 = RSQRT_COEFFS[0] + ((x2 * c3) >> 32);
157 uint64_t odd =
158 static_cast<uint64_t>((x >> 34) * static_cast<int64_t>(c4 >> 3)) + x_26;
159 // Calculate the even part of the polynomial using Horner's method.
160 uint64_t d0 = RSQRT_COEFFS[9] + ((x2 * RSQRT_COEFFS[11]) >> 32);
161 uint64_t d1 = RSQRT_COEFFS[7] + ((x2 * d0) >> 32);
162 uint64_t d2 = RSQRT_COEFFS[5] + ((x2 * d1) >> 32);
163 uint64_t d3 = RSQRT_COEFFS[3] + ((x2 * d2) >> 32);
164 uint64_t d4 = RSQRT_COEFFS[1] + ((x2 * d3) >> 32);
165 uint64_t even = 0xd105eb806655d608ul + ((x2 * d4) >> 6) + x2_26;
166
167 uint64_t r = even - odd; // error < 1.5e-10
168 // Newton-Raphson first order step to improve accuracy of the result to almost
169 // 64 bits.
170 return rsqrt_newton_raphson(m, r);
171}
172
173#else
174// Cubic minimax polynomials for 1/sqrt(x) on [1 + k/64, 1 + (k + 1)/64]
175// for k = 0..63.
176LIBC_INLINE_VAR constexpr uint32_t RSQRT_COEFFS[64][4] = {
177 {0xffffffff, 0xfffff780, 0xbff55815, 0x9bb5b6e7},
178 {0xfc0bd889, 0xfa1d6e7d, 0xb8a95a89, 0x938bf8f0},
179 {0xf82ec882, 0xf473bea9, 0xb1bf4705, 0x8bed0079},
180 {0xf467f280, 0xeefff2a1, 0xab309d4a, 0x84cdb431},
181 {0xf0b6848c, 0xe9bf46f4, 0xa4f76232, 0x7e24037b},
182 {0xed19b75e, 0xe4af2628, 0x9f0e1340, 0x77e6ca62},
183 {0xe990cdad, 0xdfcd2521, 0x996f9b96, 0x720db8df},
184 {0xe61b138e, 0xdb16ffde, 0x94174a00, 0x6c913cff},
185 {0xe2b7dddf, 0xd68a967b, 0x8f00c812, 0x676a6f92},
186 {0xdf6689b7, 0xd225ea80, 0x8a281226, 0x62930308},
187 {0xdc267bea, 0xcde71c63, 0x8589702c, 0x5e05343e},
188 {0xd8f7208e, 0xc9cc6948, 0x81216f2e, 0x59bbbcf8},
189 {0xd5d7ea91, 0xc5d428ee, 0x7cecdb76, 0x55b1c7d6},
190 {0xd2c8534e, 0xc1fccbc9, 0x78e8bb45, 0x51e2e592},
191 {0xcfc7da32, 0xbe44d94a, 0x75124a0a, 0x4e4b0369},
192 {0xccd6045f, 0xbaaaee41, 0x7166f40f, 0x4ae66284},
193 {0xc9f25c5c, 0xb72dbb69, 0x6de45288, 0x47b19045},
194 {0xc71c71c7, 0xb3cc040f, 0x6a882804, 0x44a95f5f},
195 {0xc453d90f, 0xb0849cd4, 0x67505d2a, 0x41cae1a0},
196 {0xc1982b2e, 0xad566a85, 0x643afdc8, 0x3f13625c},
197 {0xbee9056f, 0xaa406113, 0x6146361f, 0x3c806169},
198 {0xbc46092e, 0xa7418293, 0x5e70506d, 0x3a0f8e8e},
199 {0xb9aedba5, 0xa458de58, 0x5bb7b2b1, 0x37bec572},
200 {0xb72325b7, 0xa1859022, 0x591adc9a, 0x358c09e2},
201 {0xb4a293c2, 0x9ec6bf52, 0x569865a7, 0x33758476},
202 {0xb22cd56d, 0x9c1b9e36, 0x542efb6a, 0x31797f8a},
203 {0xafc19d86, 0x9983695c, 0x51dd5ffb, 0x2f96647a},
204 {0xad60a1d1, 0x96fd66f7, 0x4fa2687c, 0x2dcab91f},
205 {0xab099ae9, 0x9488e64b, 0x4d7cfbc9, 0x2c151d8a},
206 {0xa8bc441a, 0x92253f20, 0x4b6c1139, 0x2a7449ef},
207 {0xa6785b42, 0x8fd1d14a, 0x496eaf82, 0x28e70cc3},
208 {0xa43da0ae, 0x8d8e042a, 0x4783eba7, 0x276c4900},
209 {0xa20bd701, 0x8b594648, 0x45aae80a, 0x2602f493},
210 {0x9fe2c315, 0x89330ce4, 0x43e2d382, 0x24aa16ec},
211 {0x9dc22be4, 0x871ad399, 0x422ae88c, 0x2360c7af},
212 {0x9ba9da6c, 0x85101c05, 0x40826c88, 0x22262d7b},
213 {0x99999999, 0x83126d70, 0x3ee8af07, 0x20f97cd2},
214 {0x97913630, 0x81215480, 0x3d5d0922, 0x1fd9f714},
215 {0x95907eb8, 0x7f3c62ef, 0x3bdedce0, 0x1ec6e994},
216 {0x93974369, 0x7d632f45, 0x3a6d94a9, 0x1dbfacbb},
217 {0x91a55615, 0x7b955498, 0x3908a2be, 0x1cc3a33b},
218 {0x8fba8a1c, 0x79d2724e, 0x37af80bf, 0x1bd23960},
219 {0x8dd6b456, 0x781a2be4, 0x3661af39, 0x1aeae458},
220 {0x8bf9ab07, 0x766c28ba, 0x351eb539, 0x1a0d21a2},
221 {0x8a2345cc, 0x74c813dd, 0x33e61feb, 0x19387676},
222 {0x88535d90, 0x732d9bdc, 0x32b7823a, 0x186c6f3e},
223 {0x8689cc7e, 0x719c7297, 0x3192747d, 0x17a89f21},
224 {0x84c66df1, 0x70144d19, 0x30769424, 0x16ec9f89},
225 {0x83091e6a, 0x6e94e36c, 0x2f63836f, 0x16380fbf},
226 {0x8151bb87, 0x6d1df079, 0x2e58e925, 0x158a9484},
227 {0x7fa023f1, 0x6baf31de, 0x2d567053, 0x14e3d7ba},
228 {0x7df43758, 0x6a4867d3, 0x2c5bc811, 0x1443880e},
229 {0x7c4dd664, 0x68e95508, 0x2b68a346, 0x13a958ab},
230 {0x7aace2b0, 0x6791be86, 0x2a7cb871, 0x131500ee},
231 {0x79113ebc, 0x66416b95, 0x2997c17a, 0x12863c29},
232 {0x777acde8, 0x64f825a1, 0x28b97b82, 0x11fcc95c},
233 {0x75e9746a, 0x63b5b822, 0x27e1a6b4, 0x11786b03},
234 {0x745d1746, 0x6279f081, 0x2710061d, 0x10f8e6da},
235 {0x72d59c46, 0x61449e06, 0x26445f86, 0x107e05ac},
236 {0x7152e9f4, 0x601591be, 0x257e7b4d, 0x10079327},
237 {0x6fd4e793, 0x5eec9e6b, 0x24be2445, 0x0f955da9},
238 {0x6e5b7d16, 0x5dc9986e, 0x24032795, 0x0f273620},
239 {0x6ce6931d, 0x5cac55b7, 0x234d5496, 0x0ebcefdb},
240 {0x6b7612ec, 0x5b94adb2, 0x229c7cbc, 0x0e56606e},
241};
242
243// Approximate rsqrt with cubic polynomials.
244// The range [1,2] is splitted into 64 equal sub-ranges and the reciprocal
245// square root is approximated by a cubic polynomial by the minimax method in
246// each subrange. The approximation accuracy fits into 32-33 bits and thus it is
247// natural to round coefficients into 32 bit. The constant coefficient can be
248// rounded to 33 bits since the most significant bit is always 1 and implicitly
249// assumed in the table.
250LIBC_INLINE constexpr uint64_t rsqrt_approx(uint64_t m) {
251 // ULP(m) = 2^-64.
252 // Use the top 6 bits as index for looking up polynomial coeffs.
253 uint64_t indx = m >> 58;
254
255 uint64_t c0 = static_cast<uint64_t>(RSQRT_COEFFS[indx][0]);
256 c0 <<= 31; // to 64 bit with the space for the implicit bit
257 c0 |= 1ull << 63; // add implicit bit
258
259 uint64_t c1 = static_cast<uint64_t>(RSQRT_COEFFS[indx][1]);
260 c1 <<= 25; // to 64 bit format
261
262 uint64_t c2 = static_cast<uint64_t>(RSQRT_COEFFS[indx][2]);
263 uint64_t c3 = static_cast<uint64_t>(RSQRT_COEFFS[indx][3]);
264
265 uint64_t d = (m << 6) >> 32; // local coordinate in the subrange [0, 2^32]
266 uint64_t d2 = (d * d) >> 32; // square of the local coordinate
267 uint64_t re = c0 + (d2 * c2 >> 13); // even part of the polynomial (positive)
268 uint64_t ro = d * ((c1 + ((d2 * c3) >> 19)) >> 26) >>
269 6; // odd part of the polynomial (negative)
270 uint64_t r = re - ro; // maximal error < 1.55e-10 and it is less than 2^-32
271 // Newton-Raphson first order step to improve accuracy of the result to almost
272 // 64 bits.
273 r = rsqrt_newton_raphson(m, r);
274 // Adjust in the unlucky case x~1;
275 if (LIBC_UNLIKELY(!r))
276 --r;
277 return r;
278}
279#endif // LIBC_MATH_HAS_SMALL_TABLES
280
281} // namespace sqrtf128_internal
282
283LIBC_INLINE float128 sqrtf128(float128 x) {
284 using namespace sqrtf128_internal;
285 using FPBits = fputil::FPBits<float128>;
286 // Get rounding mode.
287 uint32_t rm = fputil::get_round();
288
289 FPBits xbits(x);
290 UInt128 x_u = xbits.uintval();
291 // Bring leading bit of the mantissa to the highest bit.
292 // ulp(x_frac) = 2^-128.
293 UInt128 x_frac = xbits.get_mantissa() << (FPBits::EXP_LEN + 1);
294
295 int sign_exp = static_cast<int>(x_u >> FPBits::FRACTION_LEN);
296
297 if (LIBC_UNLIKELY(sign_exp == 0 || sign_exp >= 0x7fff)) {
298 // Special cases: NAN, inf, negative numbers
299 if (sign_exp >= 0x7fff) {
300 // x = -0 or x = inf
301 if (xbits.is_zero() || xbits == xbits.inf())
302 return x;
303 // x is nan
304 if (xbits.is_nan()) {
305 // pass through quiet nan
306 if (xbits.is_quiet_nan())
307 return x;
308 // transform signaling nan to quiet and return
309 return xbits.quiet_nan().get_val();
310 }
311 // x < 0 or x = -inf
312 fputil::set_errno_if_required(EDOM);
313 fputil::raise_except_if_required(FE_INVALID);
314 return xbits.quiet_nan().get_val();
315 }
316 // Now x is subnormal or x = +0.
317
318 // x is +0.
319 if (x_frac == 0)
320 return x;
321
322 // Normalize subnormal inputs.
323 sign_exp = -cpp::countl_zero(value: x_frac);
324 int normal_shifts = 1 - sign_exp;
325 x_frac <<= normal_shifts;
326 }
327
328 // For sign_exp = biased exponent of x = real_exponent + 16383,
329 // let f be the real exponent of the output:
330 // f = floor(real_exponent / 2)
331 // Then:
332 // floor((sign_exp + 1) / 2) = f + 8192
333 // Hence, the biased exponent of the final result is:
334 // f + 16383 = floor((sign_exp + 1) / 2) + 8191.
335 // Since the output mantissa will include the hidden bit, we can define the
336 // output exponent part:
337 // e2 = floor((sign_exp + 1) / 2) + 8190
338 unsigned i = static_cast<unsigned>(1 - (sign_exp & 1));
339 uint32_t q2 = (sign_exp + 1) >> 1;
340 // Exponent of the final result
341 uint32_t e2 = q2 + 8190;
342
343 constexpr uint64_t RSQRT_2[2] = {~0ull,
344 0xb504f333f9de6484 /* 2^64/sqrt(2) */};
345
346 // Approximate 1/sqrt(1 + x_frac)
347 // Error: |r_1 - 1/sqrt(x)| < 2^-62.
348 uint64_t r1 = rsqrt_approx(m: static_cast<uint64_t>(x_frac >> 64));
349 // Adjust for the even/odd exponent.
350 uint64_t r2 = prod_hi(x: r1, y: RSQRT_2[i]);
351 unsigned shift = 2 - i;
352
353 // Normalized input:
354 // 1 <= x_reduced < 4
355 UInt128 x_reduced = (x_frac >> shift) | (UInt128(1) << (126 + i));
356 // With r2 ~ 1/sqrt(x) up to 2^-63, we perform another round of Newton-Raphson
357 // iteration:
358 // r3 = r2 - r2 * h / 2,
359 // for h = r2^2 * x - 1.
360 // Then:
361 // sqrt(x) = x * (1 / sqrt(x))
362 // ~ x * r3
363 // = x * (r2 - r2 * h / 2)
364 // = (x * r2) - (x * r2) * h / 2
365 UInt128 sx = prod_hi(x: x_reduced, y: r2);
366 UInt128 h = prod_hi(x: sx, y: r2) << 2;
367 UInt128 ds = static_cast<UInt128>(prod_hi(x: static_cast<Int128>(h), y: sx));
368 UInt128 v = (sx << 1) - ds;
369
370 uint32_t nrst = rm == FE_TONEAREST;
371 // The result lies within (-2,5) of true square root so we now
372 // test that we can correctly round the result taking into account
373 // the rounding mode.
374 // Check the lowest 14 bits (by clearing and sign-extending the top
375 // 32 - 14 = 18 bits).
376 int dd = (static_cast<int>(v) << 18) >> 18;
377
378 if (LIBC_UNLIKELY(dd < 4 && dd >= -8)) { // can round correctly?
379 // m is almost the final result it can be only 1 ulp off so we
380 // just need to test both possibilities. We square it and
381 // compare with the initial argument.
382 UInt128 m = v >> 15;
383 UInt128 m2 = m * m;
384 // The difference of the squared result and the argument
385 Int128 t0 = static_cast<Int128>(m2 - (x_reduced << 98));
386 if (t0 == 0) {
387 // the square root is exact
388 v = m << 15;
389 } else {
390 // Add +-1 ulp to m depend on the sign of the difference. Here
391 // we do not need to square again since (m+1)^2 = m^2 + 2*m +
392 // 1 so just need to add shifted m and 1.
393 Int128 t1 = t0;
394 Int128 sgn = t0 >> 127; // sign of the difference
395 Int128 m_xor_sgn = static_cast<Int128>(m << 1) ^ sgn;
396 t1 -= m_xor_sgn;
397 t1 += Int128(1) + sgn;
398
399 Int128 sgn1 = t1 >> 127;
400 if (LIBC_UNLIKELY(sgn == sgn1)) {
401 t0 = t1;
402 v -= sgn << 15;
403 t1 -= m_xor_sgn;
404 t1 += Int128(1) + sgn;
405 }
406
407 if (t1 == 0) {
408 // 1 ulp offset brings again an exact root
409 v = (m - static_cast<UInt128>((sgn << 1) + 1)) << 15;
410 } else {
411 t1 += t0;
412 Int128 side = t1 >> 127; // select what is closer m or m+-1
413 v &= ~UInt128(0) << 15; // wipe the fractional bits
414 v -= ((sgn & side) | (~sgn & 1)) << (15 + static_cast<int>(side));
415 v |= 1; // add sticky bit since we cannot have an exact mid-point
416 // situation
417 }
418 }
419 }
420
421 unsigned frac = static_cast<unsigned>(v) & 0x7fff; // fractional part
422 unsigned rnd = 0; // round bit
423 if (LIBC_LIKELY(nrst != 0)) {
424 rnd = frac >> 14; // round to nearest tie to even
425 } else if (rm == FE_UPWARD) {
426 rnd = !!frac; // round up
427 } else {
428 rnd = 0; // round down or round to zero
429 }
430
431 v >>= 15; // position mantissa
432 v += rnd; // round
433
434 // Set inexact flag only if square root is inexact
435 // TODO: We will have to raise FE_INEXACT most of the time, but this
436 // operation is very costly, especially in x86-64, since technically, it
437 // needs to synchronize both SSE and x87 flags. Need to investigate
438 // further to see how we can make this performant.
439 // https://github.com/llvm/llvm-project/issues/126753
440
441 // if(frac) fputil::raise_except_if_required(FE_INEXACT);
442
443 v += static_cast<UInt128>(e2) << FPBits::FRACTION_LEN; // place exponent
444 return cpp::bit_cast<float128>(from: v);
445}
446
447} // namespace math
448} // namespace LIBC_NAMESPACE_DECL
449
450#endif // LIBC_TYPES_HAS_FLOAT128
451
452#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SQRTF128_H
453