1//===-- Implementation header for acosf -------------------------*- 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_ACOSF_H
10#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF_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/except_value_utils.h"
16#include "src/__support/FPUtil/multiply_add.h"
17#include "src/__support/FPUtil/sqrt.h"
18#include "src/__support/macros/config.h"
19#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
20
21namespace LIBC_NAMESPACE_DECL {
22
23namespace math {
24
25namespace acosf_internal {
26
27#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
28
29LIBC_INLINE_VAR constexpr size_t N_EXCEPTS = 4;
30
31// Exceptional values when |x| <= 0.5
32LIBC_INLINE_VAR constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS =
33 {.values: {
34 // (inputs, RZ output, RU offset, RD offset, RN offset)
35 // x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
36 {.input: 0x328885a3, .rnd_towardzero_result: 0x3fc90fda, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
37 // x = -0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
38 {.input: 0xb28885a3, .rnd_towardzero_result: 0x3fc90fda, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
39 // x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ)
40 {.input: 0x39826222, .rnd_towardzero_result: 0x3fc907b4, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
41 // x = -0x1.04c444p-12, acosf(x) = 0x1.923p0 (RZ)
42 {.input: 0xb9826222, .rnd_towardzero_result: 0x3fc91800, .rnd_upward_offset: 1, .rnd_downward_offset: 0, .rnd_tonearest_offset: 1},
43 }};
44
45#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
46
47} // namespace acosf_internal
48
49LIBC_INLINE float acosf(float x) {
50 constexpr double M_MATH_PI = 0x1.921fb54442d18p+1;
51 constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
52
53 using namespace acosf_internal;
54 using namespace asin_internal;
55 using FPBits = typename fputil::FPBits<float>;
56
57 FPBits xbits(x);
58 uint32_t x_uint = xbits.uintval();
59 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
60 uint32_t x_sign = x_uint >> 31;
61
62 // |x| <= 0.5
63 if (LIBC_UNLIKELY(x_abs <= 0x3f00'0000U)) {
64 // |x| < 0x1p-10
65 if (LIBC_UNLIKELY(x_abs < 0x3a80'0000U)) {
66 // When |x| < 2^-10, we use the following approximation:
67 // acos(x) = pi/2 - asin(x)
68 // ~ pi/2 - x - x^3 / 6
69
70#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
71 // Check for exceptional values
72 if (auto r = ACOSF_EXCEPTS.lookup(x_bits: x_uint); LIBC_UNLIKELY(r.has_value()))
73 return r.value();
74#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
75
76 double xd = static_cast<double>(x);
77 return static_cast<float>(fputil::multiply_add(
78 x: -0x1.5555555555555p-3 * xd, y: xd * xd, z: M_MATH_PI_2 - xd));
79 }
80
81 // For |x| <= 0.5, we approximate acosf(x) by:
82 // acos(x) = pi/2 - asin(x) = pi/2 - x * P(x^2)
83 // Where P(X^2) = Q(X) is a degree-24 minimax even polynomial approximating
84 // asin(x)/x on [0, 0.5] generated by Sollya with:
85 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20,
86 // 22, 24|], [|1, D...|], [0, 0.5]);
87 double xd = static_cast<double>(x);
88 double xsq = xd * xd;
89 double x3 = xd * xsq;
90 double r = asin_eval(xsq);
91 return static_cast<float>(fputil::multiply_add(x: -x3, y: r, z: M_MATH_PI_2 - xd));
92 }
93
94 // |x| >= 1, return 0, 2pi, or NaNs.
95 if (LIBC_UNLIKELY(x_abs >= 0x3f80'0000U)) {
96 if (x_abs == 0x3f80'0000U)
97 return x_sign ? /* x == -1.0f */ fputil::round_result_slightly_down(
98 value_rn: 0x1.921fb6p+1f)
99 : /* x == 1.0f */ 0.0f;
100
101 if (xbits.is_signaling_nan()) {
102 fputil::raise_except_if_required(FE_INVALID);
103 return FPBits::quiet_nan().get_val();
104 }
105
106 // |x| <= +/-inf
107 if (x_abs <= 0x7f80'0000U) {
108 fputil::set_errno_if_required(EDOM);
109 fputil::raise_except_if_required(FE_INVALID);
110 }
111
112 return x + FPBits::quiet_nan().get_val();
113 }
114
115 // When 0.5 < |x| < 1, we perform range reduction as follow:
116 //
117 // Assume further that 0.5 < x <= 1, and let:
118 // y = acos(x)
119 // We use the double angle formula:
120 // x = cos(y) = 1 - 2 sin^2(y/2)
121 // So:
122 // sin(y/2) = sqrt( (1 - x)/2 )
123 // And hence:
124 // y = 2 * asin( sqrt( (1 - x)/2 ) )
125 // Let u = (1 - x)/2, then
126 // acos(x) = 2 * asin( sqrt(u) )
127 // Moreover, since 0.5 < x <= 1,
128 // 0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
129 // And hence we can reuse the same polynomial approximation of asin(x) when
130 // |x| <= 0.5:
131 // acos(x) ~ 2 * sqrt(u) * P(u).
132 //
133 // When -1 < x <= -0.5, we use the identity:
134 // acos(x) = pi - acos(-x)
135 // which is reduced to the postive case.
136
137 xbits.set_sign(Sign::POS);
138 double xd = static_cast<double>(xbits.get_val());
139 double u = fputil::multiply_add(x: -0.5, y: xd, z: 0.5);
140 double cv = 2 * fputil::sqrt<double>(x: u);
141
142 double r3 = asin_eval(xsq: u);
143 double r = fputil::multiply_add(x: cv * u, y: r3, z: cv);
144 return static_cast<float>(x_sign ? M_MATH_PI - r : r);
145}
146
147} // namespace math
148
149} // namespace LIBC_NAMESPACE_DECL
150
151#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSF_H
152