1//===-- Implementation header for acosf16 -----------------------*- 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_ACOSF16_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF16_H
11
12#include "include/llvm-libc-macros/float16-macros.h"
13
14#ifdef LIBC_TYPES_HAS_FLOAT16
15
16#include "src/__support/FPUtil/FEnvImpl.h"
17#include "src/__support/FPUtil/FPBits.h"
18#include "src/__support/FPUtil/PolyEval.h"
19#include "src/__support/FPUtil/cast.h"
20#include "src/__support/FPUtil/except_value_utils.h"
21#include "src/__support/FPUtil/multiply_add.h"
22#include "src/__support/FPUtil/sqrt.h"
23#include "src/__support/macros/optimization.h"
24
25namespace LIBC_NAMESPACE_DECL {
26
27namespace math {
28
29LIBC_INLINE constexpr float16 acosf16(float16 x) {
30
31 // Generated by Sollya using the following command:
32 // > round(pi/2, SG, RN);
33 // > round(pi, SG, RN);
34 constexpr float PI_OVER_2 = 0x1.921fb6p0f;
35 constexpr float PI = 0x1.921fb6p1f;
36
37#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
38 constexpr size_t N_EXCEPTS = 2;
39
40 constexpr fputil::ExceptValues<float16, N_EXCEPTS> ACOSF16_EXCEPTS{.values: {
41 // (input, RZ output, RU offset, RD offset, RN offset)
42 {.input: 0xacaf, .rnd_towardzero_result: 0x3e93, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 0},
43 {.input: 0xb874, .rnd_towardzero_result: 0x4052, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
44 }};
45#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
46
47 using FPBits = fputil::FPBits<float16>;
48 FPBits xbits(x);
49
50 uint16_t x_u = xbits.uintval();
51 uint16_t x_abs = x_u & 0x7fff;
52 uint16_t x_sign = x_u >> 15;
53
54 // |x| > 0x1p0, |x| > 1, or x is NaN.
55 if (LIBC_UNLIKELY(x_abs > 0x3c00)) {
56 // acosf16(NaN) = NaN
57 if (xbits.is_nan()) {
58 if (xbits.is_signaling_nan()) {
59 fputil::raise_except_if_required(FE_INVALID);
60 return FPBits::quiet_nan().get_val();
61 }
62
63 return x;
64 }
65
66 // 1 < |x| <= +/-inf
67 fputil::raise_except_if_required(FE_INVALID);
68 fputil::set_errno_if_required(EDOM);
69
70 return FPBits::quiet_nan().get_val();
71 }
72
73 float xf = x;
74
75#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
76 // Handle exceptional values
77 if (auto r = ACOSF16_EXCEPTS.lookup(x_bits: x_u); LIBC_UNLIKELY(r.has_value()))
78 return r.value();
79#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
80
81 // |x| == 0x1p0, x is 1 or -1
82 // if x is (-)1, return pi, else
83 // if x is (+)1, return 0
84 if (LIBC_UNLIKELY(x_abs == 0x3c00))
85 return fputil::cast<float16>(x: x_sign ? PI : 0.0f);
86
87 float xsq = xf * xf;
88
89 // |x| <= 0x1p-1, |x| <= 0.5
90 if (x_abs <= 0x3800) {
91 // if x is 0, return pi/2
92 if (LIBC_UNLIKELY(x_abs == 0))
93 return fputil::cast<float16>(x: PI_OVER_2);
94
95 // Note that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x)
96 // Degree-6 minimax polynomial of asin(x) generated by Sollya with:
97 // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
98 float interm =
99 fputil::polyeval(x: xsq, a0: 0x1.000002p0f, a: 0x1.554c2ap-3f, a: 0x1.3541ccp-4f,
100 a: 0x1.43b2d6p-5f, a: 0x1.a0d73ep-5f);
101 return fputil::cast<float16>(x: fputil::multiply_add(x: -xf, y: interm, z: PI_OVER_2));
102 }
103
104 // When |x| > 0.5, assume that 0.5 < |x| <= 1
105 //
106 // Step-by-step range-reduction proof:
107 // 1: Let y = asin(x), such that, x = sin(y)
108 // 2: From complimentary angle identity:
109 // x = sin(y) = cos(pi/2 - y)
110 // 3: Let z = pi/2 - y, such that x = cos(z)
111 // 4: From double angle formula; cos(2A) = 1 - 2 * sin^2(A):
112 // z = 2A, z/2 = A
113 // cos(z) = 1 - 2 * sin^2(z/2)
114 // 5: Make sin(z/2) subject of the formula:
115 // sin(z/2) = sqrt((1 - cos(z))/2)
116 // 6: Recall [3]; x = cos(z). Therefore:
117 // sin(z/2) = sqrt((1 - x)/2)
118 // 7: Let u = (1 - x)/2
119 // 8: Therefore:
120 // asin(sqrt(u)) = z/2
121 // 2 * asin(sqrt(u)) = z
122 // 9: Recall [3]; z = pi/2 - y. Therefore:
123 // y = pi/2 - z
124 // y = pi/2 - 2 * asin(sqrt(u))
125 // 10: Recall [1], y = asin(x). Therefore:
126 // asin(x) = pi/2 - 2 * asin(sqrt(u))
127 // 11: Recall that: acos(x) = pi/2 + asin(-x) = pi/2 - asin(x)
128 // Therefore:
129 // acos(x) = pi/2 - (pi/2 - 2 * asin(sqrt(u)))
130 // acos(x) = 2 * asin(sqrt(u))
131 //
132 // THE RANGE REDUCTION, HOW?
133 // 12: Recall [7], u = (1 - x)/2
134 // 13: Since 0.5 < x <= 1, therefore:
135 // 0 <= u <= 0.25 and 0 <= sqrt(u) <= 0.5
136 //
137 // Hence, we can reuse the same [0, 0.5] domain polynomial approximation for
138 // Step [11] as `sqrt(u)` is in range.
139 // When -1 < x <= -0.5, the identity:
140 // acos(x) = pi - acos(-x)
141 // allows us to compute for the negative x value (lhs)
142 // with a positive x value instead (rhs).
143
144 float xf_abs = (xf < 0 ? -xf : xf);
145 float u = fputil::multiply_add(x: -0.5f, y: xf_abs, z: 0.5f);
146 float sqrt_u = fputil::sqrt<float>(x: u);
147
148 // Degree-6 minimax polynomial of asin(x) generated by Sollya with:
149 // > P = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8|], [|SG...|], [0, 0.5]);
150 float asin_sqrt_u =
151 sqrt_u * fputil::polyeval(x: u, a0: 0x1.000002p0f, a: 0x1.554c2ap-3f,
152 a: 0x1.3541ccp-4f, a: 0x1.43b2d6p-5f, a: 0x1.a0d73ep-5f);
153
154 return fputil::cast<float16>(
155 x: x_sign ? fputil::multiply_add(x: -2.0f, y: asin_sqrt_u, z: PI) : 2 * asin_sqrt_u);
156}
157
158} // namespace math
159
160} // namespace LIBC_NAMESPACE_DECL
161
162#endif // LIBC_TYPES_HAS_FLOAT16
163
164#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF16_H
165