1//===-- Implementation header for asinbf16 ----------------------*- 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_ASINBF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ASINBF16_H
11
12#include "inv_trigf_utils.h"
13#include "src/__support/FPUtil/FEnvImpl.h"
14#include "src/__support/FPUtil/FPBits.h"
15#include "src/__support/FPUtil/bfloat16.h"
16#include "src/__support/FPUtil/cast.h"
17#include "src/__support/FPUtil/multiply_add.h"
18#include "src/__support/FPUtil/sqrt.h"
19#include "src/__support/macros/optimization.h"
20
21namespace LIBC_NAMESPACE_DECL {
22namespace math {
23
24LIBC_INLINE LIBC_CONSTEXPR bfloat16 asinbf16(bfloat16 x) {
25 // Generated by Sollya using the following command:
26 // > display = hexadecimal;
27 // > round(pi/2, SG, RN);
28 constexpr float PI_2 = 0x1.921fb6p0f;
29
30 using FPBits = fputil::FPBits<bfloat16>;
31 FPBits xbits(x);
32
33 uint16_t x_u = xbits.uintval();
34 uint16_t x_abs = x_u & 0x7fff;
35 float x_sign = (x_u >> 15) ? -1 : 1;
36 float xf = x;
37
38 float xf_abs = (xf < 0 ? -xf : xf);
39 float x_sq = xf_abs * xf_abs;
40
41 // Case 1: |x| <= 0.5
42 if (x_abs <= 0x3F00) { // x_abs <= 0.5
43 // |x| = {0}
44 if (LIBC_UNLIKELY(x_abs == 0))
45 return x; // with sign
46
47 if (LIBC_UNLIKELY(x_abs <= 0x3D00)) {
48#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
49 int rounding = fputil::quick_get_round();
50 if ((xbits.is_pos() && rounding == FE_UPWARD) ||
51 (xbits.is_neg() && rounding == FE_DOWNWARD)) {
52 return fputil::cast<bfloat16>(x: fputil::multiply_add(x: xf, y: 0x1.0p-9f, z: xf));
53 }
54#endif
55 return x;
56 }
57
58 float xp = fputil::cast<float>(x: inv_trigf_utils_internal::asin_eval(xsq: x_sq));
59 float result = xf * (fputil::multiply_add<float>(x: x_sq, y: xp, z: 1.0f));
60 return fputil::cast<bfloat16>(x: result);
61 }
62
63 // Case 2: 0.5 <|x| <= 1
64 // using reduction: asin(x) = pi/2 - 2*asin(sqrt((1-x)/2))
65 if (x_abs <= 0x3F80) { // x_abs <= 1
66 // |x| = {1}
67 if (LIBC_UNLIKELY(x_abs == 0x3F80)) {
68 return fputil::cast<bfloat16>(x: x_sign * PI_2);
69 }
70
71 float t = fputil::multiply_add<float>(x: xf_abs, y: -0.5f, z: 0.5f);
72 float t_sqrt = fputil::sqrt<float>(x: t);
73 float tp = fputil::cast<float>(x: inv_trigf_utils_internal::asin_eval(xsq: t));
74 float asin_sqrt_t = t_sqrt * (fputil::multiply_add<float>(x: t, y: tp, z: 1.0f));
75 float result = fputil::multiply_add<float>(x: -2.0f, y: asin_sqrt_t, z: PI_2);
76 return fputil::cast<bfloat16>(x: x_sign * result);
77 }
78
79 // Case 3: NaN and Inf
80 // NaN
81 if (xbits.is_nan()) {
82 if (xbits.is_signaling_nan()) {
83 fputil::raise_except_if_required(FE_INVALID);
84 return FPBits::quiet_nan().get_val();
85 }
86 return x; // quiet NaN
87 }
88 // |x|>1 & inf
89 fputil::raise_except_if_required(FE_INVALID);
90 fputil::set_errno_if_required(EDOM); // Domain is bounded
91 return FPBits::quiet_nan().get_val();
92}
93
94} // namespace math
95} // namespace LIBC_NAMESPACE_DECL
96
97#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINBF16_H
98