1 //===-- Square root of IEEE 754 floating point 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_H
10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_SQRT_H
11 
12 #include "sqrt_80_bit_long_double.h"
13 #include "src/__support/CPP/Bit.h"
14 #include "src/__support/CPP/TypeTraits.h"
15 #include "src/__support/CPP/UInt128.h"
16 #include "src/__support/FPUtil/FEnvImpl.h"
17 #include "src/__support/FPUtil/FPBits.h"
18 #include "src/__support/FPUtil/PlatformDefs.h"
19 #include "src/__support/FPUtil/builtin_wrappers.h"
20 
21 namespace __llvm_libc {
22 namespace fputil {
23 
24 namespace internal {
25 
26 template <typename T> struct SpecialLongDouble {
27   static constexpr bool VALUE = false;
28 };
29 
30 #if defined(SPECIAL_X86_LONG_DOUBLE)
31 template <> struct SpecialLongDouble<long double> {
32   static constexpr bool VALUE = true;
33 };
34 #endif // SPECIAL_X86_LONG_DOUBLE
35 
36 template <typename T>
37 static inline void normalize(int &exponent,
38                              typename FPBits<T>::UIntType &mantissa) {
39   const int shift = unsafe_clz(mantissa) -
40                     (8 * sizeof(mantissa) - 1 - MantissaWidth<T>::VALUE);
41   exponent -= shift;
42   mantissa <<= shift;
43 }
44 
45 #ifdef LONG_DOUBLE_IS_DOUBLE
46 template <>
47 inline void normalize<long double>(int &exponent, uint64_t &mantissa) {
48   normalize<double>(exponent, mantissa);
49 }
50 #elif !defined(SPECIAL_X86_LONG_DOUBLE)
51 template <>
52 inline void normalize<long double>(int &exponent, UInt128 &mantissa) {
53   const uint64_t hi_bits = static_cast<uint64_t>(mantissa >> 64);
54   const int shift = hi_bits
55                         ? (unsafe_clz(hi_bits) - 15)
56                         : (unsafe_clz(static_cast<uint64_t>(mantissa)) + 49);
57   exponent -= shift;
58   mantissa <<= shift;
59 }
60 #endif
61 
62 } // namespace internal
63 
64 // Correctly rounded IEEE 754 SQRT for all rounding modes.
65 // Shift-and-add algorithm.
66 template <typename T>
67 static inline cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, T>
68 sqrt(T x) {
69 
70   if constexpr (internal::SpecialLongDouble<T>::VALUE) {
71     // Special 80-bit long double.
72     return x86::sqrt(x);
73   } else {
74     // IEEE floating points formats.
75     using UIntType = typename FPBits<T>::UIntType;
76     constexpr UIntType ONE = UIntType(1) << MantissaWidth<T>::VALUE;
77 
78     FPBits<T> bits(x);
79 
80     if (bits.is_inf_or_nan()) {
81       if (bits.get_sign() && (bits.get_mantissa() == 0)) {
82         // sqrt(-Inf) = NaN
83         return FPBits<T>::build_nan(ONE >> 1);
84       } else {
85         // sqrt(NaN) = NaN
86         // sqrt(+Inf) = +Inf
87         return x;
88       }
89     } else if (bits.is_zero()) {
90       // sqrt(+0) = +0
91       // sqrt(-0) = -0
92       return x;
93     } else if (bits.get_sign()) {
94       // sqrt( negative numbers ) = NaN
95       return FPBits<T>::build_nan(ONE >> 1);
96     } else {
97       int x_exp = bits.get_exponent();
98       UIntType x_mant = bits.get_mantissa();
99 
100       // Step 1a: Normalize denormal input and append hidden bit to the mantissa
101       if (bits.get_unbiased_exponent() == 0) {
102         ++x_exp; // let x_exp be the correct exponent of ONE bit.
103         internal::normalize<T>(x_exp, x_mant);
104       } else {
105         x_mant |= ONE;
106       }
107 
108       // Step 1b: Make sure the exponent is even.
109       if (x_exp & 1) {
110         --x_exp;
111         x_mant <<= 1;
112       }
113 
114       // After step 1b, x = 2^(x_exp) * x_mant, where x_exp is even, and
115       // 1 <= x_mant < 4.  So sqrt(x) = 2^(x_exp / 2) * y, with 1 <= y < 2.
116       // Notice that the output of sqrt is always in the normal range.
117       // To perform shift-and-add algorithm to find y, let denote:
118       //   y(n) = 1.y_1 y_2 ... y_n, we can define the nth residue to be:
119       //   r(n) = 2^n ( x_mant - y(n)^2 ).
120       // That leads to the following recurrence formula:
121       //   r(n) = 2*r(n-1) - y_n*[ 2*y(n-1) + 2^(-n-1) ]
122       // with the initial conditions: y(0) = 1, and r(0) = x - 1.
123       // So the nth digit y_n of the mantissa of sqrt(x) can be found by:
124       //   y_n = 1 if 2*r(n-1) >= 2*y(n - 1) + 2^(-n-1)
125       //         0 otherwise.
126       UIntType y = ONE;
127       UIntType r = x_mant - ONE;
128 
129       for (UIntType current_bit = ONE >> 1; current_bit; current_bit >>= 1) {
130         r <<= 1;
131         UIntType tmp = (y << 1) + current_bit; // 2*y(n - 1) + 2^(-n-1)
132         if (r >= tmp) {
133           r -= tmp;
134           y += current_bit;
135         }
136       }
137 
138       // We compute one more iteration in order to round correctly.
139       bool lsb = y & 1; // Least significant bit
140       bool rb = false;  // Round bit
141       r <<= 2;
142       UIntType tmp = (y << 2) + 1;
143       if (r >= tmp) {
144         r -= tmp;
145         rb = true;
146       }
147 
148       // Remove hidden bit and append the exponent field.
149       x_exp = ((x_exp >> 1) + FPBits<T>::EXPONENT_BIAS);
150 
151       y = (y - ONE) | (static_cast<UIntType>(x_exp) << MantissaWidth<T>::VALUE);
152 
153       switch (get_round()) {
154       case FE_TONEAREST:
155         // Round to nearest, ties to even
156         if (rb && (lsb || (r != 0)))
157           ++y;
158         break;
159       case FE_UPWARD:
160         if (rb || (r != 0))
161           ++y;
162         break;
163       }
164 
165       return __llvm_libc::bit_cast<T>(y);
166     }
167   }
168 }
169 
170 } // namespace fputil
171 } // namespace __llvm_libc
172 
173 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_SQRT_H
174