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/asin_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
39 // For |x| <= 2^-12, 0x1.92p0 < acos(x) < 0x1.93p0: the lower bound is
40 // a bfloat16 value and the upper bound is the midpoint to its successor.
41 // Thus acos(x) rounds like pi/2 in every rounding mode. Return before
42 // squaring x or evaluating the polynomial to avoid intermediate underflow.
43 if (LIBC_UNLIKELY(x_abs <= 0x3980))
44 return fputil::cast<bfloat16>(x: PI_2);
45
46 float xf = x;
47
48 float xf_abs = (xf < 0 ? -xf : xf);
49 float x_sq = xf_abs * xf_abs;
50
51 // case 1: x <= 0.5
52 if (x_abs <= 0x3F00) {
53 float result = xf * asin_internal::asin_eval_float(xsq: x_sq);
54 return fputil::cast<bfloat16>(x: PI_2 - result);
55 }
56
57 // case 2: 0.5< |x|<= 1.0
58 if (x_abs <= 0x3F80) {
59 // |x| = {1}
60 if (x_abs == 0x3F80) {
61 if (sign)
62 return fputil::cast<bfloat16>(x: PI);
63 else
64 return FPBits::zero().get_val();
65 }
66
67 // using reduction for acos:
68 // acos(|x|) = 2*asin(sqrt((1 - |x|)/2)),
69 // and acos(x) = acos(|x|) for x >= 0, pi - acos(|x|) for x < 0
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
74 return fputil::cast<bfloat16>(
75 x: (sign) ? fputil::multiply_add(x: asin_sqrt_t, y: -2.0f, z: PI)
76 : 2.0f * asin_sqrt_t);
77 }
78 // case 3: NaN or |x| > 1
79 // NaN
80 if (xbits.is_nan()) {
81 if (xbits.is_signaling_nan()) {
82 fputil::raise_except_if_required(FE_INVALID);
83 return FPBits::quiet_nan().get_val();
84 }
85 return x; // quiet NaN
86 }
87 // |x| > 1
88 fputil::raise_except_if_required(FE_INVALID);
89 fputil::set_errno_if_required(EDOM); // Domain is bounded
90 return FPBits::quiet_nan().get_val();
91}
92
93} // namespace math
94} // namespace LIBC_NAMESPACE_DECL
95
96#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSBF16_H
97