1//===-- Implementation header for acosbf16 ----------------------*- 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_ACOSBF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSBF16_H
11
12#include "src/__support/FPUtil/FEnvImpl.h"
13#include "src/__support/FPUtil/FPBits.h"
14#include "src/__support/FPUtil/bfloat16.h"
15#include "src/__support/FPUtil/cast.h"
16#include "src/__support/FPUtil/multiply_add.h"
17#include "src/__support/FPUtil/sqrt.h"
18#include "src/__support/macros/optimization.h"
19#include "src/__support/math/inv_trigf_utils.h"
20
21namespace LIBC_NAMESPACE_DECL {
22namespace math {
23
24LIBC_INLINE bfloat16 acosbf16(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 // > round(pi, SG, RN);
30 constexpr float PI = 0x1.921fb6p1f;
31
32 using FPBits = fputil::FPBits<bfloat16>;
33 FPBits xbits(x);
34
35 uint16_t x_u = xbits.uintval();
36 uint16_t x_abs = x_u & 0x7fff;
37 bool sign = (x_u >> 15);
38 float xf = x;
39
40 float xf_abs = (xf < 0 ? -xf : xf);
41 float x_sq = xf_abs * xf_abs;
42
43 // case 1: x <= 0.5
44 if (x_abs <= 0x3F00) {
45 // |x| = {0}
46 if (LIBC_UNLIKELY(x_abs == 0))
47 return fputil::cast<bfloat16>(x: PI_2);
48
49 float xp = fputil::cast<float>(x: inv_trigf_utils_internal::asin_eval(xsq: x_sq));
50 float result = xf * fputil::multiply_add(x: x_sq, y: xp, z: 1.0f);
51 return fputil::cast<bfloat16>(x: PI_2 - result);
52 }
53
54 // case 2: 0.5< |x|<= 1.0
55 if (x_abs <= 0x3F80) {
56 // |x| = {1}
57 if (x_abs == 0x3F80) {
58 if (sign)
59 return fputil::cast<bfloat16>(x: PI);
60 else
61 return FPBits::zero().get_val();
62 }
63
64 // using reduction for acos:
65 // acos(|x|) = 2*asin(sqrt((1 - |x|)/2)),
66 // and acos(x) = acos(|x|) for x >= 0, pi - acos(|x|) for x < 0
67 float t = fputil::multiply_add<float>(x: xf_abs, y: -0.5f, z: 0.5f);
68 float t_sqrt = fputil::sqrt<float>(x: t);
69 // TODO: Use bfloat16 version for inv_trigf_utils_internals after they are
70 // available Tracking issue :
71 // https://github.com/llvm/llvm-project/issues/202079
72 float tp = fputil::cast<float>(x: inv_trigf_utils_internal::asin_eval(xsq: t));
73 float asin_sqrt_t = t_sqrt * (fputil::multiply_add(x: t, y: tp, z: 1.0f));
74
75 return fputil::cast<bfloat16>(
76 x: (sign) ? fputil::multiply_add(x: asin_sqrt_t, y: -2.0f, z: PI)
77 : 2.0f * asin_sqrt_t);
78 }
79 // case 3: NaN or |x| > 1
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
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_ACOSBF16_H
98