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 #include "src/math/hypotf.h" 9 #include "src/__support/FPUtil/BasicOperations.h" 10 #include "src/__support/FPUtil/FPBits.h" 11 #include "src/__support/FPUtil/sqrt.h" 12 #include "src/__support/common.h" 13 14 namespace LIBC_NAMESPACE { 15 16 LLVM_LIBC_FUNCTION(float, hypotf, (float x, float y)) { 17 using DoubleBits = fputil::FPBits<double>; 18 using FPBits = fputil::FPBits<float>; 19 20 FPBits x_bits(x), y_bits(y); 21 22 uint16_t x_exp = x_bits.get_biased_exponent(); 23 uint16_t y_exp = y_bits.get_biased_exponent(); 24 uint16_t exp_diff = (x_exp > y_exp) ? (x_exp - y_exp) : (y_exp - x_exp); 25 26 if (exp_diff >= FPBits::FRACTION_LEN + 2) { 27 return fputil::abs(x) + fputil::abs(y); 28 } 29 30 double xd = static_cast<double>(x); 31 double yd = static_cast<double>(y); 32 33 // These squares are exact. 34 double x_sq = xd * xd; 35 double y_sq = yd * yd; 36 37 // Compute the sum of squares. 38 double sum_sq = x_sq + y_sq; 39 40 // Compute the rounding error with Fast2Sum algorithm: 41 // x_sq + y_sq = sum_sq - err 42 double err = (x_sq >= y_sq) ? (sum_sq - x_sq) - y_sq : (sum_sq - y_sq) - x_sq; 43 44 // Take sqrt in double precision. 45 DoubleBits result(fputil::sqrt(sum_sq)); 46 47 if (!DoubleBits(sum_sq).is_inf_or_nan()) { 48 // Correct rounding. 49 double r_sq = result.get_val() * result.get_val(); 50 double diff = sum_sq - r_sq; 51 constexpr uint64_t MASK = 0x0000'0000'3FFF'FFFFULL; 52 uint64_t lrs = result.uintval() & MASK; 53 54 if (lrs == 0x0000'0000'1000'0000ULL && err < diff) { 55 result.set_uintval(result.uintval() | 1ULL); 56 } else if (lrs == 0x0000'0000'3000'0000ULL && err > diff) { 57 result.set_uintval(result.uintval() - 1ULL); 58 } 59 } else { 60 FPBits bits_x(x), bits_y(y); 61 if (bits_x.is_inf_or_nan() || bits_y.is_inf_or_nan()) { 62 if (bits_x.is_inf() || bits_y.is_inf()) 63 return FPBits::inf().get_val(); 64 if (bits_x.is_nan()) 65 return x; 66 return y; 67 } 68 } 69 70 return static_cast<float>(result.get_val()); 71 } 72 73 } // namespace LIBC_NAMESPACE 74