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
14namespace LIBC_NAMESPACE {
15
16LLVM_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(x: 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(x: 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

source code of libc/src/math/generic/hypotf.cpp