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 "asin_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 result = xf * asin_internal::asin_eval_float(xsq: x_sq);
59 return fputil::cast<bfloat16>(x: result);
60 }
61
62 // Case 2: 0.5 <|x| <= 1
63 // using reduction: asin(x) = pi/2 - 2*asin(sqrt((1-x)/2))
64 if (x_abs <= 0x3F80) { // x_abs <= 1
65 // |x| = {1}
66 if (LIBC_UNLIKELY(x_abs == 0x3F80)) {
67 return fputil::cast<bfloat16>(x: x_sign * PI_2);
68 }
69
70 float t = fputil::multiply_add<float>(x: xf_abs, y: -0.5f, z: 0.5f);
71 float t_sqrt = fputil::sqrt<float>(x: t);
72 float asin_sqrt_t = t_sqrt * asin_internal::asin_eval_float(xsq: t);
73 float result = fputil::multiply_add<float>(x: -2.0f, y: asin_sqrt_t, z: PI_2);
74 return fputil::cast<bfloat16>(x: x_sign * result);
75 }
76
77 // Case 3: NaN and Inf
78 // NaN
79 if (xbits.is_nan()) {
80 if (xbits.is_signaling_nan()) {
81 fputil::raise_except_if_required(FE_INVALID);
82 return FPBits::quiet_nan().get_val();
83 }
84 return x; // quiet NaN
85 }
86 // |x|>1 & inf
87 fputil::raise_except_if_required(FE_INVALID);
88 fputil::set_errno_if_required(EDOM); // Domain is bounded
89 return FPBits::quiet_nan().get_val();
90}
91
92} // namespace math
93} // namespace LIBC_NAMESPACE_DECL
94
95#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ASINBF16_H
96