1//===-- Implementation header for asinhf16 ----------------------*- 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_ASINHF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "acoshf_utils.h"
17#include "src/__support/FPUtil/FEnvImpl.h"
18#include "src/__support/FPUtil/FPBits.h"
19#include "src/__support/FPUtil/PolyEval.h"
20#include "src/__support/FPUtil/cast.h"
21#include "src/__support/FPUtil/except_value_utils.h"
22#include "src/__support/FPUtil/multiply_add.h"
23#include "src/__support/FPUtil/rounding_mode.h"
24#include "src/__support/FPUtil/sqrt.h"
25#include "src/__support/macros/config.h"
26#include "src/__support/macros/optimization.h"
27
28namespace LIBC_NAMESPACE_DECL {
29
30namespace math {
31
32LIBC_INLINE constexpr float16 asinhf16(float16 x) {
33
34#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
35 constexpr size_t N_EXCEPTS = 8;
36
37 constexpr fputil::ExceptValues<float16, N_EXCEPTS> ASINHF16_EXCEPTS{.values: {
38 // (input, RZ output, RU offset, RD offset, RN offset)
39
40 // x = 0x1.da4p-2, asinhf16(x) = 0x1.ca8p-2 (RZ)
41 {.input: 0x3769, .rnd_towardzero_result: 0x372a, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
42 // x = 0x1.d6cp-1, asinhf16(x) = 0x1.a58p-1 (RZ)
43 {.input: 0x3b5b, .rnd_towardzero_result: 0x3a96, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
44 // x = 0x1.c7cp+3, asinhf16(x) = 0x1.accp+1 (RZ)
45 {.input: 0x4b1f, .rnd_towardzero_result: 0x42b3, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
46 // x = 0x1.26cp+4, asinhf16(x) = 0x1.cd8p+1 (RZ)
47 {.input: 0x4c9b, .rnd_towardzero_result: 0x4336, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
48 // x = -0x1.da4p-2, asinhf16(x) = -0x1.ca8p-2 (RZ)
49 {.input: 0xb769, .rnd_towardzero_result: 0xb72a, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 1},
50 // x = -0x1.d6cp-1, asinhf16(x) = -0x1.a58p-1 (RZ)
51 {.input: 0xbb5b, .rnd_towardzero_result: 0xba96, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
52 // x = -0x1.c7cp+3, asinhf16(x) = -0x1.accp+1 (RZ)
53 {.input: 0xcb1f, .rnd_towardzero_result: 0xc2b3, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 0},
54 // x = -0x1.26cp+4, asinhf16(x) = -0x1.cd8p+1 (RZ)
55 {.input: 0xcc9b, .rnd_towardzero_result: 0xc336, .rnd_upward_offset: 0, .rnd_downward_offset: 1, .rnd_tonearest_offset: 1},
56 }};
57#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
58
59 using namespace acoshf_internal;
60 using FPBits = fputil::FPBits<float16>;
61 FPBits xbits(x);
62
63 uint16_t x_u = xbits.uintval();
64 uint16_t x_abs = x_u & 0x7fff;
65
66 if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
67 if (xbits.is_signaling_nan()) {
68 fputil::raise_except_if_required(FE_INVALID);
69 return FPBits::quiet_nan().get_val();
70 }
71
72 return x;
73 }
74
75#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
76 // Handle exceptional values
77 if (auto r = ASINHF16_EXCEPTS.lookup(x_bits: x_u); LIBC_UNLIKELY(r.has_value()))
78 return r.value();
79#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
80
81 float xf = x;
82 const float SIGN[2] = {1.0f, -1.0f};
83 float x_sign = SIGN[x_u >> 15];
84
85 // |x| <= 0.25
86 if (LIBC_UNLIKELY(x_abs <= 0x3400)) {
87 // when |x| < 0x1.718p-5, asinhf16(x) = x. Adjust by 1 ULP for certain
88 // rounding types.
89 if (LIBC_UNLIKELY(x_abs < 0x29c6)) {
90#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
91 int rounding = fputil::quick_get_round();
92 if ((rounding == FE_UPWARD || rounding == FE_TOWARDZERO) && xf < 0)
93 return fputil::cast<float16>(x: xf + 0x1p-24f);
94 if ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) && xf > 0)
95 return fputil::cast<float16>(x: xf - 0x1p-24f);
96#endif
97 return fputil::cast<float16>(x: xf);
98 }
99
100 float x_sq = xf * xf;
101 // Generated by Sollya with:
102 // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 2^-2]);
103 // The last coefficient 0x1.bd114ep-6f has been changed to 0x1.bd114ep-5f
104 // for better accuracy.
105 float p = fputil::polyeval(x: x_sq, a0: 1.0f, a: -0x1.555552p-3f, a: 0x1.332f6ap-4f,
106 a: -0x1.6c53dep-5f, a: 0x1.bd114ep-5f);
107
108 return fputil::cast<float16>(x: xf * p);
109 }
110
111 // General case: asinh(x) = ln(x + sqrt(x^2 + 1))
112 float sqrt_term = fputil::sqrt<float>(x: fputil::multiply_add(x: xf, y: xf, z: 1.0f));
113 return fputil::cast<float16>(
114 x: x_sign * log_eval(x: fputil::multiply_add(x: xf, y: x_sign, z: sqrt_term)));
115}
116
117} // namespace math
118
119} // namespace LIBC_NAMESPACE_DECL
120
121#endif // LIBC_TYPES_HAS_FLOAT16
122
123#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF16_H
124