1 //===-- Square root of x86 long double numbers ------------------*- C++ -*-===//
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_FPUTIL_GENERIC_SQRT_80_BIT_LONG_DOUBLE_H
10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_SQRT_80_BIT_LONG_DOUBLE_H
11
12 #include "src/__support/CPP/UInt128.h"
13 #include "src/__support/FPUtil/FEnvImpl.h"
14 #include "src/__support/FPUtil/FPBits.h"
15 #include "src/__support/FPUtil/PlatformDefs.h"
16 #include "src/__support/FPUtil/builtin_wrappers.h"
17
18 namespace __llvm_libc {
19 namespace fputil {
20 namespace x86 {
21
normalize(int & exponent,UInt128 & mantissa)22 inline void normalize(int &exponent, UInt128 &mantissa) {
23 const int shift =
24 unsafe_clz(static_cast<uint64_t>(mantissa)) -
25 (8 * sizeof(uint64_t) - 1 - MantissaWidth<long double>::VALUE);
26 exponent -= shift;
27 mantissa <<= shift;
28 }
29
30 // if constexpr statement in sqrt.h still requires x86::sqrt to be declared
31 // even when it's not used.
32 static inline long double sqrt(long double x);
33
34 // Correctly rounded SQRT for all rounding modes.
35 // Shift-and-add algorithm.
36 #if defined(SPECIAL_X86_LONG_DOUBLE)
sqrt(long double x)37 static inline long double sqrt(long double x) {
38 using UIntType = typename FPBits<long double>::UIntType;
39 constexpr UIntType ONE = UIntType(1)
40 << int(MantissaWidth<long double>::VALUE);
41
42 FPBits<long double> bits(x);
43
44 if (bits.is_inf_or_nan()) {
45 if (bits.get_sign() && (bits.get_mantissa() == 0)) {
46 // sqrt(-Inf) = NaN
47 return FPBits<long double>::build_nan(ONE >> 1);
48 } else {
49 // sqrt(NaN) = NaN
50 // sqrt(+Inf) = +Inf
51 return x;
52 }
53 } else if (bits.is_zero()) {
54 // sqrt(+0) = +0
55 // sqrt(-0) = -0
56 return x;
57 } else if (bits.get_sign()) {
58 // sqrt( negative numbers ) = NaN
59 return FPBits<long double>::build_nan(ONE >> 1);
60 } else {
61 int x_exp = bits.get_exponent();
62 UIntType x_mant = bits.get_mantissa();
63
64 // Step 1a: Normalize denormal input
65 if (bits.get_implicit_bit()) {
66 x_mant |= ONE;
67 } else if (bits.get_unbiased_exponent() == 0) {
68 normalize(x_exp, x_mant);
69 }
70
71 // Step 1b: Make sure the exponent is even.
72 if (x_exp & 1) {
73 --x_exp;
74 x_mant <<= 1;
75 }
76
77 // After step 1b, x = 2^(x_exp) * x_mant, where x_exp is even, and
78 // 1 <= x_mant < 4. So sqrt(x) = 2^(x_exp / 2) * y, with 1 <= y < 2.
79 // Notice that the output of sqrt is always in the normal range.
80 // To perform shift-and-add algorithm to find y, let denote:
81 // y(n) = 1.y_1 y_2 ... y_n, we can define the nth residue to be:
82 // r(n) = 2^n ( x_mant - y(n)^2 ).
83 // That leads to the following recurrence formula:
84 // r(n) = 2*r(n-1) - y_n*[ 2*y(n-1) + 2^(-n-1) ]
85 // with the initial conditions: y(0) = 1, and r(0) = x - 1.
86 // So the nth digit y_n of the mantissa of sqrt(x) can be found by:
87 // y_n = 1 if 2*r(n-1) >= 2*y(n - 1) + 2^(-n-1)
88 // 0 otherwise.
89 UIntType y = ONE;
90 UIntType r = x_mant - ONE;
91
92 for (UIntType current_bit = ONE >> 1; current_bit; current_bit >>= 1) {
93 r <<= 1;
94 UIntType tmp = (y << 1) + current_bit; // 2*y(n - 1) + 2^(-n-1)
95 if (r >= tmp) {
96 r -= tmp;
97 y += current_bit;
98 }
99 }
100
101 // We compute one more iteration in order to round correctly.
102 bool lsb = y & 1; // Least significant bit
103 bool rb = false; // Round bit
104 r <<= 2;
105 UIntType tmp = (y << 2) + 1;
106 if (r >= tmp) {
107 r -= tmp;
108 rb = true;
109 }
110
111 // Append the exponent field.
112 x_exp = ((x_exp >> 1) + FPBits<long double>::EXPONENT_BIAS);
113 y |= (static_cast<UIntType>(x_exp)
114 << (MantissaWidth<long double>::VALUE + 1));
115
116 switch (get_round()) {
117 case FE_TONEAREST:
118 // Round to nearest, ties to even
119 if (rb && (lsb || (r != 0)))
120 ++y;
121 break;
122 case FE_UPWARD:
123 if (rb || (r != 0))
124 ++y;
125 break;
126 }
127
128 // Extract output
129 FPBits<long double> out(0.0L);
130 out.set_unbiased_exponent(x_exp);
131 out.set_implicit_bit(1);
132 out.set_mantissa((y & (ONE - 1)));
133
134 return out;
135 }
136 }
137 #endif // SPECIAL_X86_LONG_DOUBLE
138
139 } // namespace x86
140 } // namespace fputil
141 } // namespace __llvm_libc
142
143 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_SQRT_80_BIT_LONG_DOUBLE_H
144