1//===-- Implementation header for asinf -------------------------*- 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_ASINHF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF_H
11
12#include "acoshf_utils.h"
13#include "src/__support/FPUtil/FPBits.h"
14#include "src/__support/FPUtil/PolyEval.h"
15#include "src/__support/FPUtil/multiply_add.h"
16#include "src/__support/FPUtil/sqrt.h"
17#include "src/__support/macros/config.h"
18#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
19
20namespace LIBC_NAMESPACE_DECL {
21
22namespace math {
23
24LIBC_INLINE float asinhf(float x) {
25 using namespace acoshf_internal;
26 using FPBits_t = typename fputil::FPBits<float>;
27 FPBits_t xbits(x);
28 uint32_t x_u = xbits.uintval();
29 uint32_t x_abs = xbits.abs().uintval();
30
31 // |x| <= 2^-3
32 if (LIBC_UNLIKELY(x_abs <= 0x3e00'0000U)) {
33 // |x| <= 2^-26
34 if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
35 return static_cast<float>(LIBC_UNLIKELY(x_abs == 0)
36 ? x
37 : (x - 0x1.5555555555555p-3 * x * x * x));
38 }
39
40 // Generated by Sollya with:
41 // > P = fpminimax(asinh(x)/x, [|0, 2, 4, 6, 8, 10, 12|], [|1, D...|],
42 // [0, 2^-3]);
43 // > dirtyinfnorm((asinh(x) - x*P)/asinh(x), [0, 2^-3]);
44 // 0x1.ee29e366e2913deff32ed8fa17f94bfe277a5babbp-62
45 constexpr double COEFFS[] = {
46 -0x1.555555555551ap-3, 0x1.333333330f782p-4, -0x1.6db6dafa7f405p-5,
47 0x1.f1c67120a7cf1p-6, -0x1.6e4b0e52674d3p-6, 0x1.10450cf441118p-6,
48 };
49
50 double x_d = x;
51 double x_sq = x_d * x_d;
52 double c0 = fputil::multiply_add(x: x_sq, y: COEFFS[1], z: COEFFS[0]);
53 double c1 = fputil::multiply_add(x: x_sq, y: COEFFS[3], z: COEFFS[2]);
54 double c2 = fputil::multiply_add(x: x_sq, y: COEFFS[5], z: COEFFS[4]);
55 double x_4 = x_sq * x_sq;
56 double x_3 = x_d * x_sq;
57 double p = fputil::polyeval(x: x_4, a0: c0, a: c1, a: c2);
58 return static_cast<float>(fputil::multiply_add(x: x_3, y: p, z: x_d));
59 }
60
61 constexpr double SIGN[2] = {1.0, -1.0};
62 double x_sign = SIGN[x_u >> 31];
63 double x_a = static_cast<double>(FPBits_t(x_abs).get_val());
64
65#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
66 // Helper functions to set results for exceptional cases.
67 auto round_result_slightly_down = [x_sign](float r) -> float {
68 return fputil::multiply_add(x: static_cast<float>(x_sign), y: r,
69 z: static_cast<float>(x_sign) * (-0x1.0p-24f));
70 };
71 auto round_result_slightly_up = [x_sign](float r) -> float {
72 return fputil::multiply_add(x: static_cast<float>(x_sign), y: r,
73 z: static_cast<float>(x_sign) * 0x1.0p-24f);
74 };
75
76 if (LIBC_UNLIKELY(x_abs >= 0x4b80'0000U)) {
77 // |x| >= 2^24
78 // We can approximate asinh(x) = sign(x) * log(2 * |x|).
79 if (LIBC_UNLIKELY(x_abs >= FPBits_t::inf().uintval())) {
80 if (xbits.is_signaling_nan()) {
81 fputil::raise_except_if_required(FE_INVALID);
82 return FPBits_t::quiet_nan().get_val();
83 }
84
85 return x;
86 }
87
88 // Exceptional cases when x > 2^24.
89 switch (x_abs) {
90 case 0x4bdd65a5: // |x| = 0x1.bacb4ap24f
91 return round_result_slightly_down(0x1.1e0696p4f);
92 case 0x4c803f2c: // |x| = 0x1.007e58p26f
93 return round_result_slightly_down(0x1.2b786cp4f);
94 case 0x4f8ffb03: // |x| = 0x1.1ff606p32f
95 return round_result_slightly_up(0x1.6fdd34p4f);
96 case 0x5c569e88: // |x| = 0x1.ad3d1p57f
97 return round_result_slightly_up(0x1.45c146p5f);
98 case 0x5e68984e: // |x| = 0x1.d1309cp61f
99 return round_result_slightly_up(0x1.5c9442p5f);
100 case 0x62f7a05a: // |x| = 0x1.ef40b4p70f
101 return round_result_slightly_up(0x1.8efc9ap5f);
102 case 0x655890d3: // |x| = 0x1.b121a6p75f
103 return round_result_slightly_down(0x1.a9a3f2p5f);
104 case 0x65de7ca6: // |x| = 0x1.bcf94cp76f
105 return round_result_slightly_up(0x1.af66cp5f);
106 case 0x6eb1a8ec: // |x| = 0x1.6351d8p94f
107 return round_result_slightly_down(0x1.08b512p6f);
108 case 0x76be09de: // |x| = 0x1.7c13bcp110f
109 return round_result_slightly_up(0x1.35569p6f);
110 case 0x7997f30a: // |x| = 0x1.2fe614p116f
111 return round_result_slightly_up(0x1.451436p6f);
112#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
113 case 0x7967ec37: // |x| = 0x1.cfd86ep115f
114 return round_result_slightly_up(0x1.43ff6ep6f);
115 case 0x58719e31: // |x| = 0x1.e33c62p49f
116 return round_result_slightly_down(0x1.1a576cp5f);
117 case 0x71699003: // |x| = 0x1.d32006p99f
118 return round_result_slightly_up(0x1.17aa2p6f);
119#endif // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
120 }
121
122 return static_cast<float>(x_sign * log_eval(x: 2.0 * x_a));
123
124 } else {
125 // Exceptional cases when x < 2^24.
126 if (LIBC_UNLIKELY(x_abs == 0x45abaf26)) {
127 // |x| = 0x1.575e4cp12f
128 return round_result_slightly_down(0x1.29becap3f);
129 }
130 if (LIBC_UNLIKELY(x_abs == 0x49d29048)) {
131 // |x| = 0x1.a5209p20f
132 return round_result_slightly_down(0x1.e1b92p3f);
133 }
134#ifndef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
135 if (LIBC_UNLIKELY(x_abs == 0x45e19b90)) {
136 // |x| = 0x1.c3372p12f
137 return round_result_slightly_down(0x1.327c5cp3f);
138 }
139#endif // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
140 }
141#else
142 if (LIBC_UNLIKELY(x_abs >= FPBits_t::inf().uintval())) {
143 if (xbits.is_signaling_nan()) {
144 fputil::raise_except_if_required(FE_INVALID);
145 return FPBits_t::quiet_nan().get_val();
146 }
147
148 return x;
149 }
150#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
151
152 // asinh(x) = log(x + sqrt(x^2 + 1))
153 return static_cast<float>(
154 x_sign * log_eval(x: x_a + fputil::sqrt<double>(
155 x: fputil::multiply_add(x: x_a, y: x_a, z: 1.0))));
156}
157
158} // namespace math
159
160} // namespace LIBC_NAMESPACE_DECL
161
162#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINHF_H
163