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