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

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