1//===-- Implementation header for hypotf16 ----------------------*- C++ -*-===//
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_MATH_HYPOTF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_HYPOTF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "src/__support/FPUtil/FEnvImpl.h"
17#include "src/__support/FPUtil/FPBits.h"
18#include "src/__support/FPUtil/cast.h"
19#include "src/__support/FPUtil/multiply_add.h"
20#include "src/__support/FPUtil/sqrt.h"
21#include "src/__support/common.h"
22#include "src/__support/macros/optimization.h"
23#include "src/__support/macros/properties/types.h"
24
25namespace LIBC_NAMESPACE_DECL {
26
27namespace math {
28
29// For targets where conversion from float to float16 has to be
30// emulated, fputil::hypot<float16> is faster
31LIBC_INLINE float16 hypotf16(float16 x, float16 y) {
32 using FloatBits = fputil::FPBits<float>;
33 using FPBits = fputil::FPBits<float16>;
34
35 FPBits x_abs = FPBits(x).abs();
36 FPBits y_abs = FPBits(y).abs();
37
38 bool x_abs_larger = x_abs.uintval() >= y_abs.uintval();
39
40 FPBits a_bits = x_abs_larger ? x_abs : y_abs;
41 FPBits b_bits = x_abs_larger ? y_abs : x_abs;
42
43 uint16_t a_u = a_bits.uintval();
44 uint16_t b_u = b_bits.uintval();
45
46 // Note: replacing `a_u >= FPBits::EXP_MASK` with `a_bits.is_inf_or_nan()`
47 // generates extra exponent bit masking instructions on x86-64.
48 if (LIBC_UNLIKELY(a_u >= FPBits::EXP_MASK)) {
49 // x or y is inf or nan
50 if (a_bits.is_signaling_nan() || b_bits.is_signaling_nan()) {
51 fputil::raise_except_if_required(FE_INVALID);
52 return FPBits::quiet_nan().get_val();
53 }
54 if (a_bits.is_inf() || b_bits.is_inf())
55 return FPBits::inf().get_val();
56 return a_bits.get_val();
57 }
58
59 float af = fputil::cast<float>(x: a_bits.get_val());
60 float bf = fputil::cast<float>(x: b_bits.get_val());
61
62 // Compiler runtime basic operations for float16 might not be correctly
63 // rounded for all rounding modes.
64 if (LIBC_UNLIKELY(a_u - b_u >=
65 static_cast<uint16_t>((FPBits::FRACTION_LEN + 2)
66 << FPBits::FRACTION_LEN)))
67 return fputil::cast<float16>(x: af + bf);
68
69 // These squares are exact.
70 float a_sq = af * af;
71 float sum_sq = fputil::multiply_add(x: bf, y: bf, z: a_sq);
72
73 FloatBits result(fputil::sqrt<float>(x: sum_sq));
74 uint32_t r_u = result.uintval();
75
76 // If any of the sticky bits of the result are non-zero, except the LSB, then
77 // the rounded result is correct.
78 if (LIBC_UNLIKELY(((r_u + 1) & 0x0000'0FFE) == 0)) {
79 float r_d = result.get_val();
80
81 // Perform rounding correction.
82 float sum_sq_lo = fputil::multiply_add(x: bf, y: bf, z: a_sq - sum_sq);
83 float err = sum_sq_lo - fputil::multiply_add(x: r_d, y: r_d, z: -sum_sq);
84
85 if (err > 0) {
86 r_u |= 1;
87 } else if ((err < 0) && (r_u & 1) == 0) {
88 r_u -= 1;
89 } else if ((r_u & 0x0000'1FFF) == 0) {
90 // The rounded result is exact.
91 fputil::clear_except_if_required(FE_INEXACT);
92 }
93 return fputil::cast<float16>(x: FloatBits(r_u).get_val());
94 }
95
96 return fputil::cast<float16>(x: result.get_val());
97}
98
99} // namespace math
100
101} // namespace LIBC_NAMESPACE_DECL
102
103#endif // LIBC_TYPES_HAS_FLOAT16
104
105#endif // LLVM_LIBC_SRC___SUPPORT_MATH_HYPOTF16_H
106