Warning: This file is not a C or C++ file. It does not have highlighting.

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
22namespace LIBC_NAMESPACE_DECL {
23namespace fputil {
24
25namespace internal {
26
27template <typename T>
28LIBC_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 static_cast<T>((T(1) << shift_length));
34}
35
36} // namespace internal
37
38template <typename T> struct DoubleLength;
39
40template <> struct DoubleLength<uint16_t> {
41 using Type = uint32_t;
42};
43
44template <> struct DoubleLength<uint32_t> {
45 using Type = uint64_t;
46};
47
48template <> 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//
106template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
107LIBC_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 = static_cast<StorageType>((r << 1)) +
211 ((tail_bits & current_bit) ? 1 : 0);
212 StorageType tmp = static_cast<StorageType>((y_new << 1)) +
213 current_bit; // 2*y_new(n - 1) + 2^(-n)
214 if (r >= tmp) {
215 r -= tmp;
216 y_new += current_bit;
217 }
218 }
219
220 bool round_bit = y_new & StorageType(1);
221 bool lsb = y_new & StorageType(2);
222
223 if (y_new >= ONE) {
224 y_new -= ONE;
225
226 if (out_exp == 0) {
227 out_exp = 1;
228 }
229 }
230
231 y_new >>= 1;
232
233 // Round to the nearest, tie to even.
234 int round_mode = quick_get_round();
235 switch (round_mode) {
236 case FE_TONEAREST:
237 // Round to nearest, ties to even
238 if (round_bit && (lsb || sticky_bits || (r != 0)))
239 ++y_new;
240 break;
241 case FE_UPWARD:
242 if (round_bit || sticky_bits || (r != 0))
243 ++y_new;
244 break;
245 }
246
247 if (y_new >= (ONE >> 1)) {
248 y_new -= ONE >> 1;
249 ++out_exp;
250 if (out_exp >= FPBits_t::MAX_BIASED_EXPONENT) {
251 if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
252 return FPBits_t::inf().get_val();
253 return FPBits_t::max_normal().get_val();
254 }
255 }
256
257 y_new |= static_cast<StorageType>(out_exp) << FPBits_t::FRACTION_LEN;
258
259 if (!(round_bit || sticky_bits || (r != 0)))
260 fputil::clear_except_if_required(FE_INEXACT);
261
262 return cpp::bit_cast<T>(y_new);
263}
264
265} // namespace fputil
266} // namespace LIBC_NAMESPACE_DECL
267
268#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_HYPOT_H
269

Warning: This file is not a C or C++ file. It does not have highlighting.

source code of libc/src/__support/FPUtil/Hypot.h