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> 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