1 //===-- Implementation of hypotf function ---------------------------------===//
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_HYPOT_H
10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_HYPOT_H
11
12 #include "BasicOperations.h"
13 #include "FEnvImpl.h"
14 #include "FPBits.h"
15 #include "rounding_mode.h"
16 #include "src/__support/CPP/bit.h"
17 #include "src/__support/CPP/type_traits.h"
18 #include "src/__support/common.h"
19 #include "src/__support/macros/config.h"
20 #include "src/__support/uint128.h"
21
22 namespace LIBC_NAMESPACE_DECL {
23 namespace fputil {
24
25 namespace internal {
26
27 template <typename T>
find_leading_one(T mant,int & shift_length)28 LIBC_INLINE T find_leading_one(T mant, int &shift_length) {
29 shift_length = 0;
30 if (mant > 0) {
31 shift_length = (sizeof(mant) * 8) - 1 - cpp::countl_zero(mant);
32 }
33 return T(1) << shift_length;
34 }
35
36 } // namespace internal
37
38 template <typename T> struct DoubleLength;
39
40 template <> struct DoubleLength<uint16_t> {
41 using Type = uint32_t;
42 };
43
44 template <> struct DoubleLength<uint32_t> {
45 using Type = uint64_t;
46 };
47
48 template <> struct DoubleLength<uint64_t> {
49 using Type = UInt128;
50 };
51
52 // Correctly rounded IEEE 754 HYPOT(x, y) with round to nearest, ties to even.
53 //
54 // Algorithm:
55 // - Let a = max(|x|, |y|), b = min(|x|, |y|), then we have that:
56 // a <= sqrt(a^2 + b^2) <= min(a + b, a*sqrt(2))
57 // 1. So if b < eps(a)/2, then HYPOT(x, y) = a.
58 //
59 // - Moreover, the exponent part of HYPOT(x, y) is either the same or 1 more
60 // than the exponent part of a.
61 //
62 // 2. For the remaining cases, we will use the digit-by-digit (shift-and-add)
63 // algorithm to compute SQRT(Z):
64 //
65 // - For Y = y0.y1...yn... = SQRT(Z),
66 // let Y(n) = y0.y1...yn be the first n fractional digits of Y.
67 //
68 // - The nth scaled residual R(n) is defined to be:
69 // R(n) = 2^n * (Z - Y(n)^2)
70 //
71 // - Since Y(n) = Y(n - 1) + yn * 2^(-n), the scaled residual
72 // satisfies the following recurrence formula:
73 // R(n) = 2*R(n - 1) - yn*(2*Y(n - 1) + 2^(-n)),
74 // with the initial conditions:
75 // Y(0) = y0, and R(0) = Z - y0.
76 //
77 // - So the nth fractional digit of Y = SQRT(Z) can be decided by:
78 // yn = 1 if 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n),
79 // 0 otherwise.
80 //
81 // 3. Precision analysis:
82 //
83 // - Notice that in the decision function:
84 // 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n),
85 // the right hand side only uses up to the 2^(-n)-bit, and both sides are
86 // non-negative, so R(n - 1) can be truncated at the 2^(-(n + 1))-bit, so
87 // that 2*R(n - 1) is corrected up to the 2^(-n)-bit.
88 //
89 // - Thus, in order to round SQRT(a^2 + b^2) correctly up to n-fractional
90 // bits, we need to perform the summation (a^2 + b^2) correctly up to (2n +
91 // 2)-fractional bits, and the remaining bits are sticky bits (i.e. we only
92 // care if they are 0 or > 0), and the comparisons, additions/subtractions
93 // can be done in n-fractional bits precision.
94 //
95 // - For single precision (float), we can use uint64_t to store the sum a^2 +
96 // b^2 exact up to (2n + 2)-fractional bits.
97 //
98 // - Then we can feed this sum into the digit-by-digit algorithm for SQRT(Z)
99 // described above.
100 //
101 //
102 // Special cases:
103 // - HYPOT(x, y) is +Inf if x or y is +Inf or -Inf; else
104 // - HYPOT(x, y) is NaN if x or y is NaN.
105 //
106 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
107 LIBC_INLINE T hypot(T x, T y) {
108 using FPBits_t = FPBits<T>;
109 using StorageType = typename FPBits<T>::StorageType;
110 using DStorageType = typename DoubleLength<StorageType>::Type;
111
112 FPBits_t x_abs = FPBits_t(x).abs();
113 FPBits_t y_abs = FPBits_t(y).abs();
114
115 bool x_abs_larger = x_abs.uintval() >= y_abs.uintval();
116
117 FPBits_t a_bits = x_abs_larger ? x_abs : y_abs;
118 FPBits_t b_bits = x_abs_larger ? y_abs : x_abs;
119
120 if (LIBC_UNLIKELY(a_bits.is_inf_or_nan())) {
121 if (x_abs.is_signaling_nan() || y_abs.is_signaling_nan()) {
122 fputil::raise_except_if_required(FE_INVALID);
123 return FPBits_t::quiet_nan().get_val();
124 }
125 if (x_abs.is_inf() || y_abs.is_inf())
126 return FPBits_t::inf().get_val();
127 if (x_abs.is_nan())
128 return x;
129 // y is nan
130 return y;
131 }
132
133 uint16_t a_exp = a_bits.get_biased_exponent();
134 uint16_t b_exp = b_bits.get_biased_exponent();
135
136 if ((a_exp - b_exp >= FPBits_t::FRACTION_LEN + 2) || (x == 0) || (y == 0))
137 return x_abs.get_val() + y_abs.get_val();
138
139 uint64_t out_exp = a_exp;
140 StorageType a_mant = a_bits.get_mantissa();
141 StorageType b_mant = b_bits.get_mantissa();
142 DStorageType a_mant_sq, b_mant_sq;
143 bool sticky_bits;
144
145 // Add an extra bit to simplify the final rounding bit computation.
146 constexpr StorageType ONE = StorageType(1) << (FPBits_t::FRACTION_LEN + 1);
147
148 a_mant <<= 1;
149 b_mant <<= 1;
150
151 StorageType leading_one;
152 int y_mant_width;
153 if (a_exp != 0) {
154 leading_one = ONE;
155 a_mant |= ONE;
156 y_mant_width = FPBits_t::FRACTION_LEN + 1;
157 } else {
158 leading_one = internal::find_leading_one(a_mant, y_mant_width);
159 a_exp = 1;
160 }
161
162 if (b_exp != 0)
163 b_mant |= ONE;
164 else
165 b_exp = 1;
166
167 a_mant_sq = static_cast<DStorageType>(a_mant) * a_mant;
168 b_mant_sq = static_cast<DStorageType>(b_mant) * b_mant;
169
170 // At this point, a_exp >= b_exp > a_exp - 25, so in order to line up aSqMant
171 // and bSqMant, we need to shift bSqMant to the right by (a_exp - b_exp) bits.
172 // But before that, remember to store the losing bits to sticky.
173 // The shift length is for a^2 and b^2, so it's double of the exponent
174 // difference between a and b.
175 uint16_t shift_length = static_cast<uint16_t>(2 * (a_exp - b_exp));
176 sticky_bits =
177 ((b_mant_sq & ((DStorageType(1) << shift_length) - DStorageType(1))) !=
178 DStorageType(0));
179 b_mant_sq >>= shift_length;
180
181 DStorageType sum = a_mant_sq + b_mant_sq;
182 if (sum >= (DStorageType(1) << (2 * y_mant_width + 2))) {
183 // a^2 + b^2 >= 4* leading_one^2, so we will need an extra bit to the left.
184 if (leading_one == ONE) {
185 // For normal result, we discard the last 2 bits of the sum and increase
186 // the exponent.
187 sticky_bits = sticky_bits || ((sum & 0x3U) != 0);
188 sum >>= 2;
189 ++out_exp;
190 if (out_exp >= FPBits_t::MAX_BIASED_EXPONENT) {
191 if (int round_mode = quick_get_round();
192 round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
193 return FPBits_t::inf().get_val();
194 return FPBits_t::max_normal().get_val();
195 }
196 } else {
197 // For denormal result, we simply move the leading bit of the result to
198 // the left by 1.
199 leading_one <<= 1;
200 ++y_mant_width;
201 }
202 }
203
204 StorageType y_new = leading_one;
205 StorageType r = static_cast<StorageType>(sum >> y_mant_width) - leading_one;
206 StorageType tail_bits = static_cast<StorageType>(sum) & (leading_one - 1);
207
208 for (StorageType current_bit = leading_one >> 1; current_bit;
209 current_bit >>= 1) {
210 r = (r << 1) + ((tail_bits & current_bit) ? 1 : 0);
211 StorageType tmp = (y_new << 1) + current_bit; // 2*y_new(n - 1) + 2^(-n)
212 if (r >= tmp) {
213 r -= tmp;
214 y_new += current_bit;
215 }
216 }
217
218 bool round_bit = y_new & StorageType(1);
219 bool lsb = y_new & StorageType(2);
220
221 if (y_new >= ONE) {
222 y_new -= ONE;
223
224 if (out_exp == 0) {
225 out_exp = 1;
226 }
227 }
228
229 y_new >>= 1;
230
231 // Round to the nearest, tie to even.
232 int round_mode = quick_get_round();
233 switch (round_mode) {
234 case FE_TONEAREST:
235 // Round to nearest, ties to even
236 if (round_bit && (lsb || sticky_bits || (r != 0)))
237 ++y_new;
238 break;
239 case FE_UPWARD:
240 if (round_bit || sticky_bits || (r != 0))
241 ++y_new;
242 break;
243 }
244
245 if (y_new >= (ONE >> 1)) {
246 y_new -= ONE >> 1;
247 ++out_exp;
248 if (out_exp >= FPBits_t::MAX_BIASED_EXPONENT) {
249 if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
250 return FPBits_t::inf().get_val();
251 return FPBits_t::max_normal().get_val();
252 }
253 }
254
255 y_new |= static_cast<StorageType>(out_exp) << FPBits_t::FRACTION_LEN;
256
257 if (!(round_bit || sticky_bits || (r != 0)))
258 fputil::clear_except_if_required(FE_INEXACT);
259
260 return cpp::bit_cast<T>(y_new);
261 }
262
263 } // namespace fputil
264 } // namespace LIBC_NAMESPACE_DECL
265
266 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_HYPOT_H
267